Re: [AMBER] Need some advice on calculating diffusion constant

From: Brian Radak <>
Date: Mon, 16 Dec 2013 09:05:15 -0500

On Sun, Dec 15, 2013 at 2:05 AM, Alan <> wrote:

> Hi Brian,
> I'm using AMBER to do some research regarding the diffusion constant of
> water. And I saw your post on the mailing list,
>, discussing which ensemble to
> use. It helps a lot, though there is still something I don't quite
> understand.
> The reason not to include any kind of heat bath in MD is to avoid the
> velocity scaling effect on the diffusion constant,

so It's best if NVE ensemble is used.

Yes. Nearly all thermostatted simulations have completely fictitious
trajectories, in the sense that they do not abide (on the microscopic
level) any physical laws that I am familiar with. NVE dynamics, on the
other hand, clearly abide Newtonian mechanics (whether or not that is the
"correct" set of physical laws to abide is of course another question

> And to make the un-ergodic NVE ensemble to be a little more ergodic, an
> NPT MD needs to be ran so different initial status can be extracted from
> it. Am I right?

What makes you think NVE trajectories are non-ergodic (or rather, any less
ergodic than any other kind of trajectory)? Perhaps this is not correct,
but I have always though of NVT and NpT simulations as simply convenience
tools that remove the task of manually selecting a high Boltzmann weight
starting point.

> But How can I decide the weight of each initial status?

They are automatically drawn with the appropriate weight by a properly
converged NVT simulation. Thus a uniform weighted average across samples
will correspond to the correct average.

> Or if I just run a single MD from the most probable initial status as you
> also suggested, how do I pick it out?
That's the exact problem, there is no way to know as far as I am aware. An
energy near the mode of the energy distribution would likely be the "right"
one, but to know that you have to run an NVT simulation!

I have ran some tests with NVE ensemble, pure water box and with protein
> with tip4p water. Here is my input file for the water box MD:
> production
> &cntrl
> imin=0,irest=1,ntx=5,
> nstlim=1000000,dt=0.002,
> ntc=2,ntf=2,
> cut=8.0,ntt=0,ntb=1,
> ntpr=1000,ntwx=5,ntwv=5,ntwr=1000,
> temp0=300.0
> /
> However, the diffusion constant calculated from the trajectories is always
> in the approximation of 4, which is very high compared with the theoretical
> value of tip4p water 3.29. This doesn't happen with NPT, so I think it must
> be a problem with NVE ensemble.
> I see that you have a 2fs time step and the relatively low default SHAKE
tolerance, does the energy for the NVE run remain constant (within ~0.1
kcal/mol)? What does your initial structure and set of velocities look
like? The temp0 value here is ignored by AMBER.

Another concern, especially if this value is older, is to check that your
model is *actually* the same. Did the original work use a finite cutoff or
an Ewald scheme? How big of a cutoff?

I'd strongly suggest looking for a paper in the literature that reports
this sort of calculation. My favorite has always been the TIP4P-Ew paper,
as it is extremely thorough and succinctly describes the full simulation
protocol. They also performed several simulations under different

Horn, *et al.* J. Chem. Phys. 2004, 120, 9665.

> Thank you.
> Alan
> -------------------------------
> Using Opera's mail client:

================================ Current Address =======================
 Brian Radak                                             :     BioMaPS
Institute for Quantitative Biology
 PhD candidate - York Research Group       :     Rutgers, The State
University of New Jersey
 University of Minnesota - Twin Cities         :     Center for Integrative
Proteomics Room 308
 Graduate Program in Chemical Physics     :     174 Frelinghuysen Road,
 Department of Chemistry                          :     Piscataway, NJ
08854-8066                                 :
Sorry for the multiple e-mail addresses, just use the institute appropriate
AMBER mailing list
Received on Mon Dec 16 2013 - 06:30:02 PST
Custom Search