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