Project

General

Profile

Bug #1339

Center of mass drift with Nose-Hoover, MTTK and md-vv

Added by ABolfazl Noorjahan about 6 years ago. Updated over 4 years ago.

Status:
New
Priority:
Normal
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Simulation of single long polymer chain shows huge center of mass drift during microsecond simulation.
The very same system has been simulated with md integrator and nose-hoover and parrinello-rahman with identical settings and no significant drift has been detected.
Here I have attached the simulation files and the architecture used for running simulation. The plot of the center of mass velocity is larger than 10 mb, so I could not attach it.
I have seen your discussion in bug#165, so thought yo might be interested in this one too.

md.mdp (764 Bytes) md.mdp ABolfazl Noorjahan, 09/19/2013 10:13 PM
run.pbs (738 Bytes) run.pbs ABolfazl Noorjahan, 09/19/2013 10:13 PM
md.gro (162 KB) md.gro ABolfazl Noorjahan, 09/19/2013 10:13 PM
topol.top (719 KB) topol.top ABolfazl Noorjahan, 09/19/2013 10:22 PM
1.ps (335 KB) 1.ps TotalTime ABolfazl Noorjahan, 09/21/2013 02:08 AM
2.ps (569 KB) 2.ps ZoomedIn ABolfazl Noorjahan, 09/21/2013 02:08 AM
3.ps (348 KB) 3.ps ZoomedIn ABolfazl Noorjahan, 09/21/2013 02:08 AM

Related issues

Related to GROMACS - Feature #1137: Proposal for integrator framework (do_md) in future GROMACSNew

History

#1 Updated by ABolfazl Noorjahan about 6 years ago

I am attaching the topol.top file in case if you need it.

#2 Updated by David van der Spoel about 6 years ago

Is it center of mass of the whole system? That is, is the box changing shape?

#3 Updated by Szilárd Páll about 6 years ago

Have you tried 4.6.3?

#4 Updated by ABolfazl Noorjahan about 6 years ago

Yes, the center of mass of the whole system (which is a polymer chain with 400 monomer) is drifting. The box shape stays cubic and with a bit oscillation due to pressure control (The density is constant with correct value for polymer chain at all specified temperatures.).

I did simulations with 4.6 on a cluster and the only version available was 4.6 and below. If you think that may cause problem I can request an upgrade and repeat the calculations.

#5 Updated by Michael Shirts about 6 years ago

What's the shortest time scale over which it is statistically significant?

ABolfazl Noorjahan wrote:

Yes, the center of mass of the whole system (which is a polymer chain with 400 monomer) is drifting. The box shape stays cubic and with a bit oscillation due to pressure control (The density is constant with correct value for polymer chain at all specified temperatures.).

I did simulations with 4.6 on a cluster and the only version available was 4.6 and below. If you think that may cause problem I can request an upgrade and repeat the calculations.

#6 Updated by ABolfazl Noorjahan about 6 years ago

Here I am attaching the Velocity of center of mass during time. I tried to zoom in to see the time scale on which the diversion starts.

#7 Updated by David van der Spoel about 6 years ago

Is your molecule cyclic?

#8 Updated by Michael Shirts about 6 years ago

David van der Spoel wrote:

Is your molecule cyclic?

It's a single polyethylene chain. I'm not sure what to make of that, other than perhaps something odd might happen with temperature groups. Note that it's using h-bonds, not all bonds.

Currently, it looks like an abrupt 1 nm vdw cutoff, which could give artifacts.

I'd like to see it run with (and I am trying a short version of this now).

nstcalctcoupl = 1
nstcalpcoupl = 1
nstcalcenergy = 1
nstcomm = 1

To remove any weird issues with properties being calculated only every N steps.

#9 Updated by Michael Shirts about 6 years ago

I started running this but realized that it may take a rather long time to check this -- it's running on my laptop at 2 ns/day, and it's not clear the issues would really be distinguishable for 1-2 days. Error cases that can be distinguished are usually better. . .

So far (20-30 ps in), it appears that the COM velocities are about 10x times larger for md-vv, using nstcomm = 100. For md/pr, it's down at 5^10-6, while for md-vv, it's more like 5^10-5, occasionally 1^10-4. This seems somewhat suspicious. I'll keep investigating.

#10 Updated by Michael Shirts about 6 years ago

Interestingly, it seems to be happening when nstcomm = 1 as well, and also when MTTK is off. If it's just that it's outputting the velocities at a different place in the loop, then that's not a big deal. I'll let it go overnight and see if the drift continues, and investigate the code further.

#11 Updated by Erik Lindahl over 5 years ago

  • Assignee changed from David van der Spoel to Michael Shirts
  • Target version changed from 4.6.x to 5.x

#12 Updated by Erik Lindahl over 4 years ago

  • Target version changed from 5.x to future

Any update on the investigation? ;-)

#13 Updated by Mark Abraham over 1 year ago

  • Related to Feature #1137: Proposal for integrator framework (do_md) in future GROMACS added

Also available in: Atom PDF