tclforce script in FEP: global variable

From: Sebastian Stolzenberg (s.stolzenberg_at_gmail.com)
Date: Tue Nov 25 2008 - 01:27:09 CST

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