From: Aoife Fogarty (aoife.fogarty_at_ens.fr)
Date: Mon Mar 07 2011 - 03:43:14 CST
Hello all,
I need to output the forces on certain atoms at each timestep in my
simulation. As a test case, I took an equilibrated box of SPC/E water,
did NVT at 300 K with PME and printed the forces on every H atom at
every timestep using loadtotalforces. (Being SPC/E water, these forces
should be purely electrostatic)
The relevant part of the tcl script is below. The entire tcl script,
output force files, and other input and output files are at
http://www.chimie.ens.fr/theorie/namd-files.php
with the output force files in the format: timestep, magnitude of
force on H, x component of force, y component, z component
The output forces (magnitude and x y and z components) are quite
noisy, making them impossible to use afterwards. This doesn't seem
normal to me. I've tested with a couple of other MD programs (YASP,
DL_POLY) and the forces should change smoothly with time.
To try and track down where this noise was coming from, I did an NVT
run of one "water" molecule in a simulation cell by itself, but with
the O and H charges set to zero in the psf file. The box length is
larger than twice the cutoff for van der Waals. There should be no
forces at all on the water molecule in this case. And indeed, all
energies except kinetic are zero.
However, loadtotalforces still outputs a force on the two hydrogen
atoms! It is about 1.5 kcal mol-1 A-1 in magnitude, and extremely
noisy with time. I suppose this is where the noise in the "proper"
simulations with normal density water and correct parameters comes
from. However, does anyone know what the source of this mysterious
force is? It seems rather strange to me. Has anyone else noticed this
before?
Best wishes,
Aoife Fogarty
relevant part of tcl script:
proc calcforces {} {
global numatoms force_freq numwatmol watmols file
set frame [getstep]
# load the forces
if {int(fmod([expr $frame-1],$force_freq))==0 && ([expr $frame-1]!=0)} {
loadtotalforces p
# get force on H atoms
for {set i 1} {$i <= $numwatmol} {incr i 1} {
if {[array exists p]} {
set h $watmols($i)
set j [expr $h + 1]
set k [expr $h + 2]
set fileHef1 [open "force_on_H$j.dat" "a"]
set fileHef2 [open "force_on_H$k.dat" "a"]
set f1 $p($j)
set f2 $p($k)
set f1x [lindex $f1 0]
set f1y [lindex $f1 1]
set f1z [lindex $f1 2]
set f2x [lindex $f2 0]
set f2y [lindex $f2 1]
set f2z [lindex $f2 2]
set forcmag1 [expr {sqrt($f1x*$f1x + $f1y*$f1y + $f1z*$f1z)} ]
set forcmag2 [expr {sqrt($f2x*$f2x + $f2y*$f2y + $f2z*$f2z)} ]
puts $fileHef1 [list $frame $forcmag1 $f1x $f1y $f1z]
puts $fileHef2 [list $frame $forcmag2 $f2x $f2y $f2z]
close $fileHef1
close $fileHef2
}
}
}
}
----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:56:43 CST