From: Jim Parker (jim.parker_at_prismsciences.com)
Date: Wed Dec 02 2009 - 00:33:37 CST
Hello,
I'm trying to generate a time variable force (harmonic perturbation)
for use in NAMD. I thought a straight-forward approach would be to
compute the perturbation using the "getstep" function in Tcl. Then
addforce to an atom-id with size =Amplitude*sin(2*pi*freq*t), 0,0
However, executing the code fails with FATAL ERROR: expected
floating-point number. (actual error below)
The problem seems to be that the TCL interpreter is not substituting a
value for "getstep" into the argument for sin() and thus the argument is
not a floating-point number. Can anyone help point me in the right
direction?
I did see that several people mentioned float-point errors in Tcl. I
think this a syntax error, but maybe related to that??
Cheers,
--Jim Parker
UTSA Physics and Astronomy
Error and code follows:
*****************************************
Reason: FATAL ERROR: expected floating-point number but got
"0.1*sin(2*3.1416*1000*getstep/1000000*2)"
while executing
"vecscale $xvec $compF"
(procedure "calcforces" line 10)
invoked from within
"calcforces"
Charm++ fatal error:
FATAL ERROR: expected floating-point number but got
"0.1*sin(2*3.1416*1000*getstep/1000000*2)"
while executing
"vecscale $xvec $compF"
(procedure "calcforces" line 10)
invoked from within
"calcforces"
Aborted
*****************************************
TCL snipet follows
tclForcesScript {
# save timestep
set t [getstep] ; #t is in 2fs steps
print "The current time is $t"
# The IDs of the atoms defining the oscillation
set aidleft 1
set aidright 2
addatom $aidleft
addatom $aidright
set PI 3.1416
# real frequency = RFfreq * 10^9 Hz
set RFfreq 1000 ; #1 THz
# RFampl in angstroms
set RFampl 0.1
# print "Freq is $RFfreq and amplitude is $RFampl"
proc calcforces {} {
global aidleft aidright PI RFfreq RFampl t
loadcoords p
# Calculate the harmonic force component
# note time = t * time/step
set compF "$RFampl*sin(2*$PI*$RFfreq*$t/1000000*2)"
# The force to be applied on each atom along x-axis (basevector)
set basevec {1. 0. 0.}
addforce $aidleft [vecscale $basevec $compF]
addforce $aidright [vecscale $basevec -1.0*$compF]
}
}
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:33 CST