Use tclforces to find systematic forces on fixed atoms?

From: Bill K (bill_at_mercury.chem.pitt.edu)
Date: Tue Oct 07 2008 - 16:20:44 CDT

Hello everyone,
     I've encountered an interesting (and yet undiscussed?) problem using
the tclforces interface. I'm trying to keep a couple of atoms fixed
(in different areas of my simulation box - they're not interacting
with each other) and measure the forces exerted upon them by the
mobile atoms of the system. I used a fixedatomfile and set my target
atoms' beta values to 1 in the fixed atom pdb to immobilize them. My
tcl script works fine, by the way, as I've tested it on mobile atoms.
 The problem is that calcforces proc prints out values of 0.0 0.0 0.0
for the force on these fixed atoms.
     When I do this calculation, the fixed atoms must be completely
restrained, instead of simply applying a harmonic restraint, since
I'm calculating a force autocorrelation function. If anyone has
tackled this problem but just didn't post it, or has a good
suggestion, I would greatly appreciate your time and consideration.
I'm thinking that calcforces is printing out the forces after the
fixed atom forces have already been set to zero (which they must be
in order for their displacement to be zero). I'll paste my tcl
script below just in case you're curious. Thank you very much for
your time.
    -Bill

set atm1 250
set atm2 4560

addatom $atm1
addatam $atm2

set Tclfreq 10
set t 0

set outfilename forces1_tcl1.out
open $outfilename w

proc calcforces {} {

global atm1 atm2 Tclfreq
global t outfilename

loadtotalforces forc

if {$t>1} {
if {[array exists forc]} {
print "Forces success"
set f1 $forc($atm1)
set f2 $forc($atm2)

print "$f1 $f2"

set outfile [open $outfilename a]
set time [expr $t/1]
puts $outfile "$time $f1 $f2"
close $outfile
}
} else {
print "Forces Failure (after waiting for 1 timestep)"
}
incr t
return
}

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:49:56 CST