[gmx-users] OPLS/AA + TIP5P, anybody?

[gigo]

Dear Chris,
I am really grateful for your help. This is what I did, with additional 
LJ terms on LP1 and LP2 of tip5p:
- 5000 steps of steepest descent with positions restraints on protein 
and flexible water (flexibility like in tip4p),
- 5000 steps of steep, no restraints, flexible water,
- 5000 steps of cg, flexible water,
- 100000 steps of MD with posres and constrained bond lengths, very weak 
temp coupling (v-rescale), starting without generating of initial 
velocities, so its heating very slowly to 300K, no pressure coupling
- 100000 steps of MD with posres, no constraints, v-rescale, nvt
- 100000 steps of MD, no posres, nvt,
- 100000 steps of MD, v-rescale, Berendsen pressure coupling,
- 100000 steps of MD, v-rescale, Parrinello-Rahman pressure coupling.

Output of this chain, after removing LJ from LP, became input for 4 
simulations just like the last one from the above chain, with or without 
posres and constraints turned on.
Results:
1) Posres off, constraints off :
"step 0"
"WARNING: Listed nonbonded interaction between particles 1157 and 1163" 
<== Why does it say "particles" and not "atoms"? Nevermind, its lysine 
on the surface of the protein, one of atoms is dimensionless hydrogen.
(...)
"step 63: Water molecule starting at atom 7923 can not be settled."
(...)
"Segmentation fault"

2) Posres off, constraints on :
Warning like above + LINCS warnings, all for charged aminoacids on 
surface of the protein
Segfault at step 63

3) Posres on, constraints off, had to add refcoord_scaling = COM :
"WARNING: Listed nonbonded interaction between particles 3075 and 3079" 
<== arginine on the surface, dimensionless hydrogen

4) Posres on, constraints on, refcoord_scaling = COM :
Same warnings, several other positively charged dimensionless hydrogens 
listed, waters could not be settled, segfault.

I tried to run few other peptides and proteins with tip5p, 100% of 
crashes.

> Also, tip5p has been used successfully with the charmm foce field:
> http://pubs.acs.org/doi/abs/10.1021/ct700053u

Yes, there are no dimensionless charged hydrogens there except of these 
on ... tip3p, tip4p and tip5p water.

> (...)
> http://pubs.acs.org/doi/abs/10.1021/ct300180w
> (...)
> Are you using halogens with oplsaa-x?

No, I use OPLSAA distributed with Gromacs.


> Standard oplsaa (non-x) and tip5p seem to be a fine combination:
> http://www.sciencedirect.com/science/article/pii/S0006349507713863
> 
> you might want to contact the authors of that paper to see if they
> ever had such problems.

Thank you, I will ask them. I am starting to be sure, though, that it is 
not possible to run any simulations with tip5p and proteins (containing 
arginine for example) without some tricks, strengthening/constraining 
bonds at least. It happened during the minimization, that the distance 
between HHxy and NHx in arginine grew up to 1.48A, while HHxy was 
landing on LP of the water. Repeatability is 100%.
Regards,
Grzegorz

> 
> Chris.
> 
> 
> Dear Chris,
> Thank you for your reply. I defined a new virtual atomtype (type D) 
> with
> LJ sigma 1.72A ( 2*0.7+1.72=3.12, which is sigma of the oxygen) and
> played a bit with epsilon till the new LJ repulsion was able to prevent
> the build up of a huge force on the oxygen while interacting with
> dimensionless charged hydrogens. Also, I copied the flexibility
> parameters from tip4p to see if it helps in minimization before I turn
> it into rigid water - it seems that it does. I was able to minimize the
> system with such water. Also, I minimized the system with tip4p and
> replaced it with tip5p with a script. I tried to minimize the system
> afterwards with true tip5p, which did not work. My question is, besides
> the correctness of water model, why do you think it is safe to remove
> the LJ on lone electron pairs in MD? Will it not collapse like in 
> energy
> minimization?
> Best Regards,
> Grzegorz