Project

General

Profile

Bug #1348

OPLS/AA not compatible with TIP5P

Added by Grzegorz Wieczorek almost 4 years ago. Updated about 2 years ago.

Status:
Closed
Priority:
Normal
Category:
documentation
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

It is not possible to run successful energy minimization or MD simulation of protein containing positively charged dimensionless hydrogens (like in OPLS/AA) in TIP5P water model. During EM, these hydrogens collapse with lone electron pairs of the water, penetrating the repulsive LJ space of water's oxygen, which leads to building up of huge forces (orders of magnitude up to 1.0e+25). MD simulations with OPLS/AA and TIP5P (of systems that were previously minimized and equilibrated with any working model, replaced later to TIP5P) stop prematurely with segmentation fault preceded by warnings about violaton of listed nonbonded interactions (pointing dimensionless hydrogens as cause of the problem), failures of settling water molecules, and if constraints on bond lengths are used - LINCS warnings (for charged aminoacids on surface of the protein). Repeatability is 100%, tried with several proteins and peptides, with and without constraints and restraints. I discussed this with Chris Neale on gmx-users (title "OPLS/AA + TIP5P, anybody?"). Here I'm attaching the simplest inputs for reproducing the phenomena. It is not gromacs bug per se, but I think that the manual or on-line help should not suggest TIP5P as compatible with OPLSAA, as it does.

em00.mdp (530 Bytes) em00.mdp Grzegorz Wieczorek, 10/01/2013 03:03 PM
peptide.pdb (13.6 KB) peptide.pdb Grzegorz Wieczorek, 10/01/2013 03:03 PM
do (279 Bytes) do run this to reproduce Grzegorz Wieczorek, 10/01/2013 03:03 PM
do (286 Bytes) do run this to reproduce Grzegorz Wieczorek, 10/02/2013 08:19 AM
em00.mdp (542 Bytes) em00.mdp Grzegorz Wieczorek, 10/02/2013 08:19 AM
peptide.pdb (13.7 KB) peptide.pdb Grzegorz Wieczorek, 10/02/2013 08:19 AM
view.gro (756 KB) view.gro Grzegorz Wieczorek, 10/02/2013 08:19 AM
tip5p.itp (1.73 KB) tip5p.itp Rossen Apostolov, 06/12/2014 05:23 PM

Associated revisions

Revision 88e39531 (diff)
Added by Rossen Apostolov about 3 years ago

Added tip5p to the list of watermodels in some FF.

Added support for flexible tip5p to OPLSAA and parameters
were taken from the Amber/Charmm tip5p.itp files.

There are, however, various problems with that water model.
This commit doesn't resolve them so the bug report below
will need to be followed up at later stage.

watermodels.dat now includes a note about these issues
and refers users to the redmine report.

Refs #1348.

Change-Id: I3adaae21f35efbb4240716e91f1bb5d256c8853e

History

#1 Updated by Justin Lemkul almost 4 years ago

I have tried to follow the thread on gmx-users, but I have had nothing particularly interesting to add. Some thoughts, though:

1. It seems as if pure water runs stably, if I remember one of the posts correctly. Is this the case?
2. Can you reproduce the problem with other force fields, or just OPLS-AA?
3. Can you point to the place in the manual (PDF or online) where it is suggested that OPLS is compatible with TIP5P? There are plenty of options offered with each force field due to a generic selection mechanism, but that doesn't indicate any suggestion or endorsement beyond the one recommended water model for each force field.

#2 Updated by Grzegorz Wieczorek almost 4 years ago

Ad 1) Yes - the pure water runs stably.
Ad 2) By now I havent reproduced it for other ff, but I'm starting to check amber family
Ad 3) I apologize - there is nothing like this mentioned in the documentation. However, TIP5P is part of the OPLS/AA in gromacs - which is standard actualy. As Chris Neale suggested, maybe it would be worth to warn a user at least, for example when she/he chooses OPLS/AA and TIP5P in pdb2gmx?

#3 Updated by Justin Lemkul almost 4 years ago

Grzegorz Wieczorek wrote:

Ad 1) Yes - the pure water runs stably.

Good.

Ad 2) By now I havent reproduced it for other ff, but I'm starting to check amber family

Good, thanks.

Ad 3) I apologize - there is nothing like this mentioned in the documentation. However, TIP5P is part of the OPLS/AA in gromacs - which is standard actualy. As Chris Neale suggested, maybe it would be worth to warn a user at least, for example when she/he chooses OPLS/AA and TIP5P in pdb2gmx?

We avoid trying to out-think the user when it comes to combinations. This was an issue some time ago and we decided against any such warnings. That would just introduce a complex decision structure with pdb2gmx getting worse than it currently is ;) It's better to get to the root of the problem, and probably a good test is to see if the issue is (1) with TIP5P and any force field, (2) the specific implementation in OPLS-AA, or (3) handling of virtual sites in general. Does TIP4P work? It uses different OPLS atom types but the principle is the same, a massless, charge-bearing particle with no LJ terms. I know people have used TIP4P and OPLS quite often, so demonstrating that it still works would be good evidence that something is specifically wrong with TIP5P.

#4 Updated by Grzegorz Wieczorek almost 4 years ago

Justin Lemkul wrote:

It's better to get to the root of the problem, and probably a good test is to see if the issue is (1) with TIP5P and any force field,

I have ran several simulations with amber94 and TIP5P already, the error did not occur. In all amber ffs, proteins have hydrogens without LJ in serine and threonine. Parameters of arginine in OPLS and serine in amber94 (charges, sigmas, epsilons, bond lengths) do not really convince me that OPLS is much more prone to this error, so I will run simulations for other ambers as well. Charmm has been used with TIP5P by many, nobody complained - it does not have charged hydrogens without LJ.

 the specific implementation in OPLS-AA,

It could be done in Tinker, but I have never used it...

 handling of virtual sites in general. Does TIP4P work?

TIP4P works perfectly well for the same cases. The virtual site is buried deeply in the oxygen there is no way that any other atom can approach it.

#5 Updated by Grzegorz Wieczorek almost 4 years ago

I ran energy minimization for all forcefield that include tip5p (oplsaa, all ambers and charmm27). To my surprise, in ambers and charmm27 tip5p can be treated as flexible. Thus - 2 tests for each ff: flexible and rigid. The system consisted of a peptide, water and ions (0.15 mol). Results for opls are known already.
With all ambers and charmm27 simulations went smoothly with rigid water.
I could not minimize the system if the water was flexible for any of the above forcefields. In case of all ambers, the OH hydrogen of serine was collapsing with the LP of the water.
Now, in the case of charmm27 the problem was different. It was showing itself not in a proximity of the protein. Bonds of at least one water molecule were stretched with giant force during EM. This includes distances between oxygen and ... lone electron pairs (virtual sites)... I don't know how to interpret this.
4 files attached on 2nd of October in the morning are needed to reproduce. After simulation finishes, take a look at bond lengths of the water molecule with biggest force on it. view.gro is confout.gro that I got, after trjconv -pbc mol -ur compact. In my case its 701SOL that is badly deformed.
Conclusion for now is that TIP5P should be rigid, as it was from the beginning of this model. How the flexibility terms, taken probably from TIP4P, ended up in TIP5P? Stretching the distance between atom and virtual site growing out of it is a completely different story.

#6 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1348.
Uploader: Rossen Apostolov ()
Change-Id: I3adaae21f35efbb4240716e91f1bb5d256c8853e
Gerrit URL: https://gerrit.gromacs.org/3587

#7 Updated by Rossen Apostolov over 3 years ago

I fixed the model in 4-6, it will be merged in 5-0 also. In the meantime, can you check with the attached itp file whether it works properly?

#8 Updated by Mark Abraham over 3 years ago

Rossen Apostolov wrote:

I fixed the model in 4-6, it will be merged in 5-0 also. In the meantime, can you check with the attached itp file whether it works properly?

That change makes it possible to use a flexible TIP5P with OPLS/AA, but the issues here seem to be that the rigid TIP5P does not work with OPLS/AA, and that flexible TIP5P did not work with EM for any force field. I can find no reason for that, except perhaps that oplsaa.ff/tip5p.itp does not have a mass column and charmm27.ff and amber99sb.ff do. All three force fields have the right masses for their respective atom types. tip5p.gro seems sane.

#9 Updated by Rossen Apostolov over 3 years ago

tip5p doesn't have masses but neither does tip3p and tip4p or any of the others. we can keep the issue open

#10 Updated by Rossen Apostolov over 3 years ago

  • Status changed from Fix uploaded to In Progress

#11 Updated by Rossen Apostolov about 3 years ago

  • Status changed from In Progress to Blocked, need info

88e3953148d396001df1154a04fa39efadb19c49 is merged but the problem need more thorough investigation.

#12 Updated by Berk Hess about 2 years ago

Flexible TIP5P sounds dangerous to me, since the lone pair positions are very sensitive to position of the masses, so it wouldn't surprise mail if that causes problems.
Normal rigid TIP5P should work fine, I don't see how that could lead to issues.

#13 Updated by Erik Lindahl about 2 years ago

  • Status changed from Blocked, need info to Rejected

OK - finally had time to investigate this. Unfortunately (well, depending on your p-o-v...) this is not a Gromacs issue, but a problem in the OPLS force field.

A pure TIP5p simulation will work fine in OPLS too, but combining it with a protein is fragile. The problem is that OPLS does not use any Lennard-Jones interactions on the hydrogens e.g. in Arginine, and the LJ radius of the nitrogen in the NH2 groups is also smaller compared to Charmm27 (Amber99 also uses hydrogens with a small LJ radius in Arginine). This means the lone pairs of TIP5p (which only has LJ on the oxygen) can get extremely close to the NH2 hydrogens in Arginine, which leads to the gigantic forces and subsequent crash.

The only way to get around this is to change the force field parameters. However, apart from the fact that we don't want to alter data in a force field away from their published values, OPLS has explicitly been parametrized using TIP4p water, and one obviously doesn't want to ruin that case.

For now I guess the obvious solution is to simply stay away from TIP5p with OPLS (since the OPLS protein force field is NOT optimized for it, and is likely to be unstable in more programs if they use the same parameters).

If you still really want to use TIP5p with OPLS, I would advice to alter the NH2 hydrogen atomtype just-so-slightly by copying Amber. If you set the OPLS_301 nonbonded data to

 opls_301   H3   1   1.00800     0.460       A    1.06908e-01  6.56888e-02

then the files in this issue appear to work fine. It's a very small change compared to the electrostatic interactions for this atom, so I can't imagine it will have any bad side effects, even if you would formally be simulating something slightly different than official OPLS-AA.

#14 Updated by Erik Lindahl about 2 years ago

  • Status changed from Rejected to Closed

#15 Updated by Chris Neale about 2 years ago

Perhaps we should add a note or a warning to pdb2gmx for this combination?

However, the authors of the previously noted 2007 Biophysical Journal article titled "Structural and Hydration Properties of the Partially Unfolded States of the Prion Protein" ( http://www.sciencedirect.com/science/article/pii/S0006349507713863 ) have confirmed that their gromacs simulations of this force field combination were stable in 50 ns simulations starting from structures that were pre-equilibrated with a different force field (personal communication), so perhaps some equilibration procedures using different initial force field combinations do allow opls/tip5p simulations, or perhaps their arginine residues were simply protected from solvent by salt bridges, etc. I agree with the general gromacs philosophy of not out-thinking the user, but this seems like a sufficiently problematic force field combination to warrant some type of user notification.

Also available in: Atom PDF