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 (
http://ambermd.org/tutorials/advanced/tutorial22/section1.htm), 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 restraints.
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 (
http://archive.ambermd.org/201210/0430.html). 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 error
NSTEP ENERGY RMS GMAX NAME NUMBER
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:
Min1.in
Min2.in
Min3.in
&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,
restraintmask=':1-376',
restraint_wt = 10.0,
&end
/
&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 /
Heat1.in
Heat2.in
Equil.in
Heating System NVT 0.5 ns
&cntrl
ntx=1, irest=0, nmropt=1,
ntb=1,
cut=10.0,
ntpr=500, ntwr=500, ntwx=500, iwrap=1,
ntc=2, ntf=2, ntb=1, cut=10.0, nsnb=20,
nstlim=250000,dt=0.002,
nscm=500, ntt=1,
temp0=0.0, tempi=0.0, tautp=0.5
ntr=1,
restraintmask=':1-376',
restraint_wt=10.0,
&end
&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
heat
&cntrl
irest=1,ntx=5,
cut=10.0, ntb=2, ntp=1, taup=1.0,
nstlim=250000,dt=0.002,
ntc=2,ntf=2,
ntpr=500, ntwx=500,
ntt=3, gamma_ln=2.0,
temp0=300.0,iwrap=1,
ntr=1, restraintmask=':1-376',
restraint_wt=5.0,
/
/
equil NPT 5 ns
&cntrl
irest=1,ntx=5,
nstlim=2500000,dt=0.002,
cut=10.0, ntb=2, ntp=1, taup=2.0,
ntc=2,ntf=2,ig=-1,
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.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Nov 19 2015 - 03:00:03 PST
