Re: tclforce script in FEP: global variable

From: Jerome Henin (jhenin_at_cmm.chem.upenn.edu)
Date: Tue Nov 25 2008 - 08:57:24 CST

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