Extracting force on subset of atoms from tclForces

From: Nate Walkins (walkinsnate_at_gmail.com)
Date: Tue Dec 07 2021 - 04:07:11 CST

I recently ran a long trajectory for an all-atom system and now wish to
compute the force acting only on a subset of the atoms. I set about doing
this with tclForces by reading in a dcd files and then calling a custom
calcforces script. I read that the force can only be used computed after a
step of MD has been run as it computes the forces on the "stale"
coordinates so I thought to get around this by running each step twice
(inefficient, but I care about the forces and not speed in this case) and
passing a flag variable the second time. Eventually this is just going to
be a loop, but for now I just copied and pasted it twice to see if it was
possible to continue in this manner.

The relevant part of the namd configuration files are below. I am getting
the error message that "FATAL ERROR: can't read "p": no such variable" and
if commenting out that line then "FATAL ERROR: can't read "f": no such
variable", however both seem to be defined above.

tclForces on
tclForcesScript forces.tcl

coorfile open dcd short.dcd
coorfile read
set flag 0
set ts 0
firstTimestep $ts
run 0
set flag 1
run 0

coorfile read
set flag 0
set ts [expr $ts+1]
firstTimestep $ts
run 0
set flag 1
run 0
The accompanying tclForces script is:
set inStream [open cg-beads.ind r]
set atoms {}
foreach line [split [read $inStream] \n] {
    set ind [string trim [string range $line 0 3]]
    lappend atoms $ind
close $inStream
puts $atoms

proc calcforces {} {
    global atoms flag
    set coord_step [getstep]
    print $flag $coord_step
    if {$flag > 0} {
        loadforces f
        loadcoords p
        print $p
        print $f
cg-beads.ind is just a text file of the atom indices that I wish to get the
forces from.

This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:12 CST