From: Jacqueline Schmidt (jacqueline_schmidt_1983_at_yahoo.com)
Date: Wed Dec 26 2012 - 12:04:23 CST
Dear NAMD users,
I trying to add an smd force on a group of atoms (a peptide). I have wrote my calcforces routine as below. The problem is the if statement (to get initial center of mass of my groups of atoms) does not work! I just need them for the force calculations, can anyone please help me! why the if { $t==0} (which must occur at the beginning of the simulation) does not work?
set grp2 {}
for {set j 1 } { $j <=107 } {incr j } {
lappend grp2 $j
}
set a2 [addgroup $grp2]
# set the output frequency, initialize the time counter
set Tclfreq 50
set t 0
# force constant (kcal/mol/A^2)
set k 10
# pulling velocity (A/timestep)
set v 0.0001
set outfilename klpgwsg_smd_tcl.out
open $outfilename w
proc calcforces {} {
global Tclfreq t k v a2 c2x c2y c2z outfilename
loadcoords coordinate
if { $t==0 } {
set r20 $coordinate($a2)
set r20z [lindex $r20 2]
#print
}
set r2 $coordinate($a2)
set r2x [lindex $r2 0]
set r2y [lindex $r2 1]
set r2z [lindex $r2 2]
set f2x 0
set f2y 0
set f2z [expr $k*($v*$t-$r2z+$r20z)]
lappend f2 $f2x $f2y $f2z
addforce $a2 $f2
set foo [expr $t % $Tclfreq]
if { $foo == 0 } {
set outfile [open $outfilename a]
set time [expr $t*2/1000.0]
puts $outfile "$time $r2z $f2z"
close $outfile
}
incr t
return
}
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:22:23 CST