Re: [AMBER] Distorted and Failed Phosphate Simulations

From: <>
Date: Thu, 19 Nov 2015 14:23:03 +0100 (CET)

Hi Sigurd,
do you have some non-zero van der Waals radius on hydrogen atom? In this
case, you cannot use standard HO atom type with zero VdW radius as there is
attractive (albeit scaled as it is 1-4 interaction) electrostatic
interaction between hydrogen of OH group and another oxygen bound to
phosphorus, while there is no VdW interaction that would protect the
hydrogen from coming closer and closer to that oxygen. In other cases, the
hydrogen is "protected" by VdW radius of the OH oxygen, but in this case the
oxygen-oxygen interaction is 1-3 interaction that is excluded, so it does
not count.

hope this helps,


Pavel Banáš
Department of Physical Chemistry, 
Palacky University Olomouc 
Czech Republic 
---------- Původní zpráva ----------
Od: Sigurd Friis Truelsen <>
Komu: <>
Datum: 19. 11. 2015 11:51:33
Předmět: [AMBER] Distorted and Failed Phosphate Simulations
"Hello Amber mailing list.
I have problems running regular MD simulations with inorganic phosphate HPO4
(2-), H3PO4(-) and H3PO4 . The objective is to set up simulations where the 
anion is in complex with a phosphate binding protein, and to study their 
interaction. I follow the approach of the accelerated MD tutorial found here
(, although I 
have made some minor changes especially in the simulation times.
1. Minimize the water around the protein with a restraint of 10 kcal/mol.
2. Relax the water with an NPT for 0.5 ns, restraining the protein and 
phosphate with 10 kcal/mol.
3. The water and protein are minimized together, keeping the protein 
restrained with 10 kcal/mol.
4. The system is then heated from 0 K to 300 in 0.5 ns in an NVT ensemble 
with 5 kcal/mol restraint on the protein and phosphate.
5. The system is then equilibrated at 300 K in 0.5 ns in an NPT ensemble 
with 5 kcal/mol restrain.
6. The system is equilibrated at 300 K in 5 ns in an NPT ensemble with no 
The system is made with explicit water in an octahedral box of size 12 Å, 
counter ions are added until the system is neutral, which is only one single
Na+ ion. The cutoff is 10 Å, and I use the force field ff14SB for the 
protein and GAFF for the anion. I have set up the HPO4 ion with the 
antechamber tool, and found no error messages in the .frcmod file. Finally 
the simulations are run with pmemd.
I can perfectly run simulations with the protein by itself and also the 
protein phosphate complex when applying QM on the anion. However, when I 
remove the QM treatment on the phosphate it is broken or distorted during 
the NPT relaxation in the second simulation step. Initially I used a 2 fs 
time step with SHAKE enabled, during the second step (NPT for 0.5 ns) where 
I get the error
Coordinate resetting cannot be accomplished,
deviation is too large
iter_cnt, my_bond_idx, i and j are : 1 1 5327 5330
The atoms 5327-5330 are the oxygens of the phosphate anion. By visualizing 
the trajectory the phosphate anion does not attain the tetrahedral shape 
that it should, and over a few hundred steps the hydrogen begins to 
circulate with increasing speed until the simulation stops.
I searched for solutions in the Amber community, I found that reducing the 
time step to 0.001 or 0.0005 and removing shake (setting ntc = ntf = 1) may 
solve the problem ( Setting dt 
to 0.001 , the simulation would advance through the second simulation step 
with a distorted phosphate, and failed in the third simulation step with the
307 -2.3981E+07 3.4735E+10 7.7191E+12 H 5330
BOND = 14499.1961 ANGLE = 765.0943 DIHED = 4071.9067
VDWAALS = 25468.8996 EEL = -266120.4952 HBOND = 0.0000
1-4 VDW = 1108.5824 1-4 EEL = ************* RESTRAINT = 0.0000
Eventually the BOND is also written out as ************* , and the 
simulation stops.
I tested the 0.0005 time step for phosphate without the protein, since it 
should be stable by itself. I set up the same kind of simulations in a water
box with counter ions to make it charge neutral. However these simulations 
are (perhaps unsurprisingly) not stable, for any of the 0.5 fs, 1 fs or 2 fs
time step. The HPO4(2-) is distorted already during the first minimization.
I also tested the influence of the charge of the anion, and set simulations 
with the H2PO4(-) and H3PO4 together with the protein. However the same 
thing occurs, as the anions suddenly are distorted and the simulation stops 
either during the minimization step 3, or due to an SHAKE error during the 
heating or equilibration.
It seems that there is an issue with the force field compatibility with the 
phosphate. I checked the gaff force field (gaff.dat) and as far as I 
understand, the phosphate should be parametrized: ”Phosphate with four 
connected atoms, such as O=P(OH)3”. So I figured at least the H3PO4 anion 
should work, but maybe I am misunderstanding something.
Perhaps I can do something to stabilize the system? Any suggestions of what 
can be tried to fix this problem would be greatly appreciated.
Here are the in-files for the shake enabled simulation:
&cntrl imin=1, ntmin=1, nmropt=0, ntb=1, cut=10.0, drms=0.1 maxcyc=10000, 
ncyc=5000, ntx=1, irest=0, ntpr=100, ntwr=100, ntwx=100, iwrap=0, ntf=1, 
nsnb=20, ntr=1, restraintmask="!:WAT", restraint_wt = 10.0, &end /
&cntrl ntx = 1, irest = 0, nmropt = 0,
ntb = 2, ntp = 1 , taup = 1.0,
cut = 10.0,
ntrx = 1,
ntxo = 1,
ntpr = 1, ntwx = 1, ntwr = 1,
nscm = 2500,
ntf = 2, ntc = 2,
dt = 0.002,
nsnb = 20,
nstlim = 250000,
iwrap = 1,
t = 0.0,
temp0 = 300.0,
tempi = 200.0,
tautp = 0.5,
ntt = 1,
tol = 0.00001,
ntr = 1,
restraint_wt = 10.0,
&cntrl imin=1, ntmin=1, nmropt=0 ntb=1, cut=10.0, drms=0.1 maxcyc=2000, ncyc
=1500, ntx=1, irest=0, ntpr=1, ntwr=1, ntwx=1, iwrap=0, ntf=1, ntb=1, cut=
10.0, nsnb=20, restraintmask="!:WAT", restraint_wt = 10.0, &end /
Heating System NVT 0.5 ns
ntx=1, irest=0, nmropt=1,
ntpr=500, ntwr=500, ntwx=500, iwrap=1,
ntc=2, ntf=2, ntb=1, cut=10.0, nsnb=20,
nscm=500, ntt=1,
temp0=0.0, tempi=0.0, tautp=0.5
&wt type='REST', istep1 = 0, istep2=0, value1=1.0, value2=1.0, &end
&wt type='TEMP0', istep1=0, istep2=250000, value1=0.0, value2=300, &end
&wt type='END'
equil NTP 0.5 ns
cut=10.0, ntb=2, ntp=1, taup=1.0,
ntpr=500, ntwx=500,
ntt=3, gamma_ln=2.0,
ntr=1, restraintmask=':1-376',
equil NPT 5 ns
cut=10.0, ntb=2, ntp=1, taup=2.0,
ntpr=1000, ntwx=1000,
ntt=3, gamma_ln=2.0,
temp0 = 300.0,
Med venlig hilsen / best regards
Sigurd Friis Truelsen
PhD student
The Biomimetics Group
DTU Environment
AMBER mailing list"
AMBER mailing list
Received on Thu Nov 19 2015 - 05:30:05 PST
Custom Search