Hello Sigurd,
I didn't carefully look at all of your input parameters, but from what I
see I would also say that it is a force field problem. You could look up
the original Gaff paper and see what they say about modeling phosphate
molecules. I have done a little bit of this as well for methylated
phosphates, and used the p5 atom type that you mention as my phosphorous,
with o, oh and os as my oxygen atoms. However, in each of my molecules I
ended up doing some parameterization. My guess is that you will need to do
some as well. If I have some time, I will try to look into this more and
will respond again if I find something interesting.
Bests,
Karl
On Thu, Nov 19, 2015 at 11:48 AM, Sigurd Friis Truelsen <sigut.env.dtu.dk>
wrote:
>
> 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
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Nov 19 2015 - 03:30:04 PST
