From: Aoife Fogarty (aoife.fogarty_at_ens.fr)
Date: Tue Mar 15 2011 - 13:15:04 CDT
Hi all,
Hugh, thank you very much for the pointer. Yes, indeed, the noise in
the forces turned out to be a combination of random forces from
Langevin dynamics, and not doing a full electrostatics evaluation
every timestep. (The other MD programs I mentioned use a Berendsen
thermostat.)
However, even in NVE loadtotalforces still outputs a (not noisy) force
on the two hydrogens and the oxygen, less than 1 kcal mol-1, in my
test case of a "water" molecule alone in a box with all charges set to
zero, as I described before. The total force on the molecule is zero,
therefore these forces are obviously intramolecular forces, even
though the molecule is rigid. (I'm using "rigidBonds all", and indeed,
when I check the output trajectory file, the molecule really is rigid.)
In my NVE simulation of a box of SPC/E water, again, the total force
on each molecule corresponds to that calculated by another MD program
(DL_POLY, where water molecules are implemented as rigid bodies) but
the individual forces on the atoms do not. Again, this suggests there
are intramolecular forces, even though the molecule is rigid. These
are present whether I use SETTLE or SHAKE to constrain the water
molecules.
It was my understanding that with SHAKE, the velocities and positions
at t+dt are first calculated for the unconstrained system (using
velocities, positions and forces at time t), and then adjusted to
satisfy each constraint. So, when the forces at time t were
calculated, it was from a configuration where all the bond lengths and
angles were at their equilibrium (constrained) values, i.e. bond and
angle forces should have been zero. So where are these intramolecular
forces in rigid water coming from? Does anyone know if the forces I
just described are not actually the forces loadtotalforces outputs?
If not, does anyone know of a way to implement genuine rigid bodies in
NAMD, rather than using SHAKE or SETTLE? (Note- I'm not talking about
fixing atoms, but about rigid molecules whose centre of mass undergoes
dynamics in the usual way but which have fixed internal coordinates.)
Thanks!
(And thanks again, Hugh.)
Aoife Fogarty
Thanks againQuoting Hugh Heldenbrand <helde010_at_umn.edu>:
> 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 message was sent using IMP, the Internet Messaging Program.
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2012 - 23:19:56 CST