Hello,
A co-worker and I have been struggling to implement a simple molecular
dynamics package in order to train a new type of classical water models.
We have implemented SETTLE, a very fast real-space non-bonded force
computation, support for many different charge and parameter sets (as
well as different potentials), and two working thermostats (Nose-Hoover
and Berendsen). There is also a Nose-Hoover thermo-barostat in place
which works to properly regulate the system volume in response to the
user-supplied external pressure (1 bar) in the ideal gas case when the
virials have no role.
However, when I try to include a significant virial for something like
liquid water, the system blows up on me. I'm trying to compare to AMBER
to see what's going wrong, but getting information out of a code that only
does PME-based non-bonded forces for periodic systems and applying it to a
code that currently works on periodic systems but does not yet include
Ewald is painful.
Still, I've made some progress with comparing my results to AMBER. I
have some short questions that hopefully will keep me on the right track.
1.) What are the formats of the information in the "-o mdout" file? They
seem pretty self-explanatory, except for:
PRESSURE: kilopascals?
All energy units: kcal? (the VDWAALS term certainly is given in kcal)
VIRIAL: kcal?
2.) I tried a case where I created a 216-water system in an appropriate
periodic box. Using a 9A cutoff, I noted that the VDWAALS energy I
obtained differed from AMBER's by about 1% (I'm not ready to compare
electrostatics in periodic systems yet). I tried a similar, 1728-water
system with an 18A cutoff and found that the energies were similar to
within 0.15%. What could be going on here? Is there also an Ewald
treatment of long-range LJ interactions that I'm missing? Please note
that as of right now my code only implements real-space interactions.
3.) I assumed that the energies matched up sufficiently and compared the
virials. To do this, I stripped the system of all charges so that only
VDWAALS interactions were taking place. In two different cases
(different sizes of water boxes, with different numbers of particles
charge-free water particles), I found that my virial computation was
almost exactly twice what AMBER's was: twice that value to within 0.25%.
I've tried looking through the sander code, but can anyone state
succinctly what all is going into the "VIRIAL = ###" number reported in
the "-o mdout" file? I compute the atomic virial using the equation given
in the DL-POLY 3 manual, sum (-r_ij dot f_ij). I also compute a
constraint virial, also by the formula given in the DL-POLY 3 manual, but
it seems to be very small compared to the atomic virial.
Thanks for your help! This is the last project I intend to do as a grad
student. If it succeeds, I'm hoping it will really push the envelope of
conventional, fixed-charge MD simulations.
Dave
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Jan 31 2007 - 06:07:47 PST
