Re: Re: tclforce script in FEP: global variable

From: Sebastian Stolzenberg (s.stolzenberg_at_gmail.com)
Date: Tue Nov 25 2008 - 18:21:12 CST

Dear Jerome,

good news: my problem mentioned below is not a bug and can be fixed easily:

Someone in our lab pointed out possible conflicts between the
"tclforcescript" and "fep.tcl".
In particular, my error message

variable "at" is not defined

arises, whenever I call "runFEPlist()" or "runFEP()" in "fep.tcl", in
which my $at variable is not defined. I now believe that these
procedures directly call "calcforces()" (because they contain a "run"
command) and thus, $at cannot be found.

Simply adding

global at k r0

to the beginning of the definitions for "runFEPlist()" and "runFEP()" fixes any issues and tclforces and FEP-NAMD seems to work just fine - so far.

Does my reasoning make sense?

Cheers,
Sebastian

Jerome Henin wrote:
> Dear Sebastian,
> Thanks for reporting this. It confirms an issue that has been brought
> to the attention of namd developpers recently, and is not solved yet.
> When the error occurs, does namd exit gracefully or with a segfault?
> Jerome
>
> On Tue, Nov 25, 2008 at 2:27 AM, Sebastian Stolzenberg
> <s.stolzenberg_at_gmail.com> wrote:
>
>> Dear Jerome,
>>
>> is it possible that tclforce scripts have to be set up in NAMD-FEP slightly
>> differently than in plain MD?
>>
>> In the attached script, I inserted a tclforcescript part that worked in MD,
>> but now doesn't in FEP.
>> In particular, in the procedure calcforces(), it complains that the global
>> variable "at" is not defined.
>> So the first
>> print "at"
>> print $at
>> in the script works, at the second, it crashes.
>>
>> Thank you,
>> Best,
>> Sebastian
>>
>>
>>>> and the real story is that we are
>>>> using your FEP NAMD implementation on single-residue mutated proteins to
>>>> get
>>>> an end structure for further MD simulations. Sure enough, MD without FEP
>>>> in
>>>> principle should lead to the same results, but it seems that people in
>>>> my
>>>> lab have experienced faster MD convergence using your FEP
>>>> implementations
>>>> prior to MD production runs.
>>>>
>> Okay, that makes sense. In this case, you don't have to worry about a
>> nicely converged free energy change. I would still mutate the ion at
>>
>>
>>
>>
>> source parameter.tcl
>>
>> set input0dir ../../prep_FEP/build/transporters/$transporter/input
>> set inputdir ../../prep_FEP/build/transporters/$transporter/output
>> set outputdir ../output/$transporter
>>
>> #############################################################
>> ## ADJUSTABLE PARAMETERS ##
>> #############################################################
>>
>> set temperature 310
>> set inputname $input0dir/$transporter
>> set outputname $transporter\_fep
>> set firststep 0
>>
>> structure $inputdir/$transporter\_fep.psf
>> coordinates $inputdir/$transporter\_fep.fep
>> #bincoordinates $inputname.restart.coor
>> #binvelocities $inputname.restart.vel
>> #extendedSystem $inputname.restart.xsc
>> extendedSystem $inputname.restart.xsc
>> restartname $outputdir/$outputname.restart
>> restartsave no
>> #firsttimestep $firststep
>>
>>
>> #############################################################
>> ## SIMULATION PARAMETERS ##
>> #############################################################
>>
>> # Input
>> paraTypeCharmm on
>> parameters $common/toppar/par_all27_prot_lipid.inp
>> temperature $temperature
>>
>> # Force-Field Parameters
>> exclude scaled1-4
>> 1-4scaling 1.0
>> cutoff 12.
>> switching on
>> switchdist 10.
>> pairlistdist 13.5
>> margin 1.5
>>
>> # Integrator Parameters
>> timestep 1.0 ;# 1fs/step
>> rigidBonds all ;# needed for 2fs steps
>> rigidTolerance 0.00000001
>> nonbondedFreq 2
>> fullElectFrequency 4
>> stepspercycle 20
>>
>>
>> # Constant Temperature Control
>> langevin on ;# do langevin dynamics
>> langevinDamping 1 ;# damping coefficient (gamma)
>> langevinTemp $temperature
>> langevinHydrogen off ;# don't couple langevin bath to hydrogens
>>
>> LangevinPiston on
>> LangevinPistonTarget 1.01325
>> LangevinPistonPeriod 200
>> LangevinPistonDecay 50
>> LangevinPistonTemp $temperature
>>
>>
>> wrapAll on
>>
>> # PME (for full-system periodic electrostatics)
>> PME yes
>> PMEGridSizeX 108
>> PMEGridSizeY 108
>> PMEGridSizeZ 108
>>
>>
>> # Constant Pressure Control (variable volume)
>> useGroupPressure yes ;# needed for rigidBonds
>> useFlexibleCell yes
>> useConstantArea no
>>
>> # Output
>> outputName $outputdir/$outputname
>>
>> restartfreq 5000
>> dcdfreq 5000
>> xstFreq 5000
>> outputEnergies 5000
>> outputPressure 5000
>>
>> binaryoutput yes
>> binaryrestart yes
>> outputTiming 5000
>>
>> #constraints on
>> #consRef $inputdir/$transporter\_fep.restraints.pdb
>> #consKFile $inputdir/$transporter\_fep.restraints.pdb
>> #consKCol B
>>
>>
>> tclforces on
>>
>> tclforcesScript {
>> set at {77334 77337}
>>
>>
>> #********
>> print "at"
>> print $at
>> #********
>>
>>
>>
>> for {set i 0} {$i <= 1} {incr i} {
>> addatom [lindex $at $i]
>> }
>>
>> set k {10}
>> set r0 {0.000}
>>
>> proc calcforces {} {
>> global at k r0
>> print "TCL calcforces()"
>>
>>
>> #********
>> print "at"
>> print $at
>> #********
>>
>>
>>
>> loadcoords c
>> set r01 [vecsub $c([lindex $at 1]) $c([lindex $at 0]) ]
>> set rl01 [veclength $r01 ]
>> set n01 [vecnorm $r01]
>> set force1 [expr [lindex $k 0]*($rl01-[lindex $r0 0])]
>>
>> set f0 [vecscale $force1 $n01 ]
>> set f1 [vecscale [expr 0 - $force1] $n01 ]
>>
>> addforce [lindex $at 0] $f0
>> addforce [lindex $at 1] $f1
>> }
>>
>> proc veclength {v} {
>> set retval 0
>> foreach term $v {
>> set retval [expr $retval + $term * $term ]
>> }
>> return [expr sqrt($retval)]
>> }
>> proc vecnorm {v} {
>> set sum 0
>> foreach term $v {
>> set sum [expr $sum + $term*$term]
>> }
>> set sum [expr sqrt($sum)]
>> set retval {}
>> foreach term $v {
>> lappend retval [expr $term / $sum]
>> }
>> return $retval
>> }
>> }
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> #############################################################
>> ## EXECUTION SCRIPT ##
>> #############################################################
>>
>> source fep.tcl
>>
>> fep on
>> fepFile $inputdir/$transporter\_fep.fep
>> fepCol B
>> fepOutFile $outputdir/$outputname.fepout
>> fepOutFreq 100
>> fepEquilSteps 5000
>>
>> set numSteps 100000
>>
>> runFEPlist {.0 .0 } 1000000
>>
>> #5 intervals
>> runFEPlist {.000001 .00001 .0001 .001 .003 .005 .007 .01 } $numSteps
>> #45 intervals
>> runFEP .01 .13 .005 $numSteps
>> runFEP .13 .39 .01 $numSteps
>> runFEP .39 .99 .015 $numSteps
>> #15 intervals
>> runFEP .99 .999 .003 $numSteps
>> runFEP .999 .9999 .0003 $numSteps
>> runFEP .9999 .99999 .00003 $numSteps
>> runFEP .99999 .999999 .000003 $numSteps
>> runFEP .999999 .9999999 .0000003 $numSteps
>> runFEP .9999999 1 .0000001 1000000
>>
>>
>>
>>
>
>

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:50:09 CST