Project

General

Profile

Bug #908

PME with non-linear virtual sites gives incorrect virial and pressure

Added by Berk Hess over 7 years ago. Updated over 6 years ago.

Status:
Closed
Priority:
High
Assignee:
Category:
mdrun
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

The contribution of the PME mesh potential to the virial is determined separately and differently from all other virial contributions, which are determined directly from the particle coordinates and forces. With virtual sites forces get redistributed, which can affect the virial. Linear virtual sites (as in TIP4P) do not alter the virial, but non-linear virtual sites (as in TIP5P or when removing hydrogen vibrations) do. This effect was not corrected for in all Gromacs versions up to 4.5.5. The effect is small and isotropic for unordered systems. But it can lead to artifacts in ordered systems, such as all-atom lipid bilayer with hydrogens replaced with virtual sites (e.g. 20% increase of the area per lipid).

data.tar.gz (1.6 MB) data.tar.gz All the files packed in one tar.gz Bastien Loubet, 04/15/2013 03:14 PM
area_per_lipid.jpg (94.2 KB) area_per_lipid.jpg Various result with and without virtual sites Bastien Loubet, 05/08/2013 06:25 PM

Associated revisions

Revision 2c59fdcc (diff)
Added by Berk Hess over 7 years ago

added correction for PME virial with non-linear virtual sites

With non-linear virtual site constructs the virial nees to be corrected
when it is not calculated from the particle coordinates and masses.
Currently this is only a problem for the PME mesh contribution.
This would lead to an incorrect virial and pressure. The error was small
for unordered systems, but could be significant for e.g. an all-atom
bilayer with hydrogens replaced by virtual sites. Fixes #908

Change-Id: I27a06019aa09734a14f3630ef296aaaacb9c5f0b

History

#1 Updated by Berk Hess over 7 years ago

  • Status changed from New to Closed

Committed a proper fix to the release-4-5-patches branch for inclusion in the 4.5.6 release.

#2 Updated by Bastien Loubet over 6 years ago

We have tried the virtual site on POPC lipids with GROMACS 4.6.1. There is still an increase of the area per lipid, I am not sure if it is actually a bug. It is smaller than the one obtained with GROMACS 4.5.5 but it is still significant (~8%).

I join a figure demonstrating the problem (as a jpg).
I also join the mdp file used in the simulation, the modified CHARMM36 force field we used and a sample membrane with and without virtual sites.

Note that we have defined the virtual sites of the N-CH_3 groups in the POPC head groups using the OPLS parameters already present (both aminoacids.vsd and ffbonded.itp have been modified).
Also we did not get any significant rms change between our lipids with and without virtual sites.

If somebody could confirm if this is a bug or not that would be great.

#3 Updated by Berk Hess over 6 years ago

Do I see correctly that you used virtual sites without any bonds constraints?
Virtual sites can only be used with constraints=all-bonds.
Using no constraints at all is even problematic without v-sites.
Didn't grompp give you a warning about the too large time step?

#4 Updated by Bastien Loubet over 6 years ago

Do I see correctly that you used virtual sites without any bonds constraints?

Yes.

Virtual sites can only be used with constraints=all-bonds.

I did not know.

Using no constraints at all is even problematic without v-sites.
Didn't grompp give you a warning about the too large time step?

No, grompp gives me 4 notes, the same it gives me for the system without virtual sites. Note that in order to rule out an effect of the time step, both simulations with and without virtual sites have been done at the same 2 fs time step.
Here is the grompp output:

Ignoring obsolete mdp entry 'title'

Back Off! I just backed up mdout.mdp to ./#mdout.mdp.4#

NOTE 1 [file sim_eq_vs.mdp]:
nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting
nstcomm to nstcalcenergy

NOTE 2 [file sim_eq_vs.mdp]:
The Berendsen thermostat does not generate the correct kinetic energy
distribution. You might want to consider using the V-rescale thermostat.

NOTE 3 [file sim_eq_vs.mdp]:
The switch/shift interaction settings are just for compatibility; you
will get better performance from applying potential modifiers to your
interactions!

Generated 21528 of the 21528 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1
Generated 18355 of the 21528 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'POPC'
Excluding 2 bonded neighbours molecule type 'SOL'
Velocities were taken from a Maxwell distribution at 298 K
Cleaning up constraints and constant bonded interactions with virtual sites
Converted 23 Bonds with virtual sites to connections, 110 left
Removed 30 U-Bs with virtual sites, 226 left
Removed 4 Proper Dih.s with virtual sites, 461 left
Analysing residue names:
There are: 128 Other residues
There are: 5688 Water residues
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Number of degrees of freedom in T-Coupling group SOL is 34126.11
Number of degrees of freedom in T-Coupling group POPC is 19966.89
Largest charge group radii for Van der Waals: 0.079, 0.079 nm
Largest charge group radii for Coulomb: 0.079, 0.079 nm

NOTE 4 [file sim_eq_vs.mdp]:
The sum of the two largest charge group radii (0.157915) is larger than
rlist (1.300000) - rvdw (1.200000)

Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 56x56x72, spacing 0.116 0.116 0.107
Estimate for the relative computational load of the PME mesh part: 0.06
This run will generate roughly 885 Mb of data

There were 4 notes

gcq#15: "Don't Grumble, Give a Whistle !" (Monty Python)

Note that in the non virtual site run there are constraints on the hydrogen bonds, as in the CHARMM36 paper.
I am relaunching the virtual site calculation with all bonds constrained. I will update as soon as I have the new run.

#5 Updated by Berk Hess over 6 years ago

The vsite mdp you uploaded has no constraints at all.
If you actually used h-bonds constraints, than that explains why you didn't get a warning.
Warning for constraints=h-bonds combined with v-sites is hard. I can see if this can be done without to much effort.
Constraining all bonds will probably fix your issue.

#6 Updated by Bastien Loubet over 6 years ago

I am not sure this belong to a redmine issue but I said in my previous post that I was going to update the issue once I have the all atoms constrained simulation.
I join the plot of the area per lipid we obtained, there is still a difference between the virtual site and the non virtual site simulation. I am starting to suspect that maybe the bond length in the virtual site, which have been determined using OPLS, give some artifact for lipid simulation ? Any other idea ?
I guess at some point you must have a trade-off between performance and precision.

#7 Updated by Berk Hess over 6 years ago

  • Target version deleted (4.5.6)
  • Affected version set to 4.5.6

The difference of using virtual sites versus all bonds constrained should be negligible, I would think. Using incompatible parameters could have a significant effect, so I would check by how much the bond lengths and angles differ between OPLS and CHARMM.

Also available in: Atom PDF