Re: [AMBER] Distorted and Failed Phosphate Simulations

From: Jason Swails <>
Date: Tue, 24 Nov 2015 08:17:54 -0500

On Tue, Nov 24, 2015 at 5:08 AM, Sigurd Friis Truelsen <>

> Hi again,
> Thank you all for your valuable input. I have made some tests based on
> your suggestions, where I increased the angle force constants as suggested
> by Karl and David. I am not particularly concerned about the phosphate
> conformations during the MD simulation – as I will use QMMM in this case,
> but I would like to avoid QM in order to increase my sampling. So I am
> looking for solutions to do this. Also I am mostly interested in the
> HPO4(2-) anion.
> When increasing the angle force constants (with 15 kcal/mol rad2) I get
> the following results:
> The H2PO4(-) in the protein complex completes without any errors, and it
> looks fine in the trajectory. However it is not the case for HPO4(2-)
> (which you also pointed out is less stable). I tested for 1 and 2 fs time
> steps.

​What about trying larger force constants for HPO4(2-)?

The 1 fs time step is able to complete the minimization, heating and
> equilibration, but during the heating the phosphate is distorted as it
> interacts with a Na+ counterions which is in close proximity. The
> hydrogen-free oxygens orient toward the Na+ ion and the OH group points
> away. To me this seems like a charge interaction issue. I think it is
> necessary to have charge neutrality in the system, but only to get a stable
> simulation. When I build the system in leap, the counterions are placed
> just next to the phosphate, and I don't know if I should move them and stop
> including them?

​Charge neutrality is formally required by PME (if you think about it, it's
pretty clear that if this is *not* the case, then the electrostatic energy
of an infinitely periodic system would be infinity due to the divergent
nature of the infinite harmonic series).​ However, in practice a
neutralizing plasma is applied uniformly to the charge grid to prevent this
from happening in PME. That said, omitting positively charged ions because
they interact with a negatively charged species in a way you don't want is
poor justification for omitting them in this case (at least in my opinion).

The 2 fs time step is suddenly distorted during the equilibration, where
> the hydrogen is attracted towards one of the other oxygens and the
> simulation stops due to a shake error. This could be a van der Waals
> interaction problem that Pavel suggested, as the hydrogen is unprotected. I
> have used the gaff.dat hydroxyl group, which should have a VDW radius of
> 1.008 if I understand the parameter file correctly.

​You are looking at the mass. The radius is near the bottom of gaff.dat:

  ho 0.0000 0.0000 OPLS Jorgensen,

The well-depth is 0, which means it does not have any L-J interactions at
all (regardless of what Rmin/2 is set to).

> If this is the standard HO atom with insufficient VDW radius, which HO
> atom could I use, or would it require a full parametrization?

​No, you just need *some* L-J term to add some repulsive force. For
instance, this is what GLYCAM_06j.dat has for its hydroxyl H atom type:

  Ho 0.2000 0.0300 M.B. Tessier 2011

Give those a try and see if it works.


Jason M. Swails
Rutgers University
Postdoctoral Researcher
AMBER mailing list
Received on Tue Nov 24 2015 - 05:30:06 PST
Custom Search