Re: [AMBER] Drifting molecules despite strong restraints

From: Jason Swails <>
Date: Fri, 20 Mar 2015 11:23:15 -0400

On Fri, Mar 20, 2015 at 10:18 AM, Sanmeet Chahal <>

> Hello all,
> Thank-you for the reply Jason. I am interested in the non-bonded only
> dynamics so I think ntf = 7 is the right choice.

​I'm curious what application this has. Note that you are, in effect,
inventing your own force field that nobody has used before and that has
never been validated. You are *not* using ff14SB or ff99SB or whatever
leaprc you sourced in tleap (or anything that even looks like it)... it is
something entirely of your own creation.

In fact, this is a very weird potential. You are ignoring all bonds and
angles between 1-2 and 1-3 atom pairs, but you are ALSO ignoring their
electrostatic and van der Waals interactions... so those atoms have no
impact on one another except through many-body effects (i.e., two
"non-interacting" particles can "interact" through a mediating atom). This
is due to the fact that these atom pairs are in each others' exclusion
lists... I *think*. I'm really not entirely sure that my description fully
details what is being done.

Note that these ntf options are *intended* to be used alongside
constraints. If you are constraining a certain degree of freedom, there's
no reason to calculate a spring force between them (since the net force
will be 0 thanks to Newton's 3rd law). If you are *not* constraining those
degrees of freedom, then skipping their calculation is... weird (and wrong,
from perspective of the resulting ensemble compared to the parametrization
of the force field).

> Also I have changed to
> doing NVT heating over 20 ps which maintains the atoms in their current
> position as you said. However at the end of this process the water density
> is too low (around 0.63 g/cc) so I think that means I should do an NPT
> equilibration to bring the water density to its equilibrium value, but this
> again produces the drifting issue.

​That's why I said do a short heating, then switch to NPT. Barostats can
behave strangely when the system is at very low temperatures, so it's best
to turn them on when the temperature is closer to the target temperature.

All the .in, inpcrd, prmtop and .out (equil-file.out shows the water
> density being equilibrated):
> Here are the steps I take: 1) Minimize, 2) Heat at NVT for 20 ps, 3)
> Equilibrate at NPT
> Here are the commands I used:
> sander -O -i -o min-file.out -p file.prmtop -c file.inpcrd \
> -r min-file.rst -inf min-file.mdinfo -ref file.inpcrd
> sander -O -i -o heat-file.out -p file.prmtop -c min-file.rst -r
> heat-file.rst \
> -x heat-file.mdcrd -inf heat-file.mdinfo -ref file.inpcrd
> sander -O -i e <> -o equil-file.out -p file.prmtop
> -c
> heat-file.rst -r equil-file.rst \
> -x equil-file.mdcrd -inf equil-file.mdinfo -ref
> file.inpcrd
> Is there a way I can equilibrate the water density while keeping the
> restrained molecules in the same position?

​No. Think about what you're asking it to do. The way pressure scaling is
implemented is that the center of mass of every molecule is translated
either in toward the center of the box (when compressing, so everything
gets a *little* closer together), or you are spreading the molecules apart
by scaling them away from the center of the box (when expanding, so
everything gets a *little* further apart). The restraints are implemented
so as to respect the barostat. If you didn't scale the restraints
alongside the actual coordinates when you did the pressure scaling, then
you would introduce an impulse from the restraint forces every time the box
was scaled. In turn, this impulse would affect whatever barostat you were
using and introduce pressure artifacts that really, in my opinion, should
not exist.


Jason M. Swails
Rutgers University
Postdoctoral Researcher
AMBER mailing list
Received on Fri Mar 20 2015 - 08:30:04 PDT
Custom Search