From: Hugh Heldenbrand (helde010_at_umn.edu)
Date: Mon Mar 07 2011 - 07:10:09 CST
If you are using a thermostat that should be putting random forces on
the atoms, shouldn't it?
-Hugh Heldenbrand
Graduate Student
University of Minnesota
On 03/07/2011 03:43 AM, Aoife Fogarty wrote:
> 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 : Mon Dec 31 2012 - 23:19:53 CST