Project

General

Profile

Bug #2020

Possible issue with md-vv integrator

Added by Miguel Caro over 3 years ago. Updated almost 3 years ago.

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

Description

I am filing this issue as advised by Mark Abraham. The used version is 4.6.5. I attach my .tpr file and paste the email I sent to the mailing list below:

Hi all,

I am simulating a mixture of water and acetonitrile with the OPLS force
field. I am using the md-vv integrator because I want my velocities to
be in synchrony with my positions (md outputs velocities and positions
half a time step apart).

My box is relatively small, 2400 atoms, and I am using 1.3 nm for all
cutoffs. The properties of the mixture evolve smoothly with composition
until the point where I reach the pure acetonitrile extreme, which also
corresponds to the largest box size in the series. L increases with
acetonitrile mole fraction from about 2.86 nm for pure water to 3.39 nm
for pure acetonitrile (T = 323K). For pure acetonitrile the smooth trend
breaks down abruptly and I can observe a severe kinetic energy transfer
from the translational degrees of freedom to the vibrational degrees of
freedom. However, if I increase the cutoff to 1.6 nm the issue goes away.

Also, if I use md as integrator instead of md-vv, with the original set
of cutoffs at 1.3 nm, the issue also goes away. I would like to
understand how/why/whether the choice of integrator and the box size
affect the effect of the cutoffs on the distribution of kinetic energy
in the system. I've gone as far as pinpointing the possible origin of
the issue and a workable solution (i.e. just switch to md), but I would
like to understand the underlying cause of the problem.

Many thanks,
Miguel

md.tpr (82.4 KB) md.tpr Miguel Caro, 07/29/2016 10:19 AM
gromacs_test.zip (150 KB) gromacs_test.zip Miguel Caro, 08/01/2016 05:06 PM

Related issues

Related to GROMACS - Feature #2070: Physical validation testingIn Progress

History

#1 Updated by Mark Abraham over 3 years ago

Thanks. What is L? What is in this particular tpr? What should one be able to observe with it? Why does it have periodic molecules in it if it's a simple liquid mixture? Are you using suitable virtual sites to model your linear nitrile moiety?

The 4.6.x series has been out of maintenance for about two years now, and many subsequent bugs fixed could be relevant, so I would strongly encourage you to use a version that's up to date. We strive to have no simulation types that seem to work with one integrator and not with another, but we're human :-)

I would also encourage you to use the Verlet scheme, rather than unbuffered cutoffs in the group scheme, though it should not have any bearing on the issue at hand.

#2 Updated by Miguel Caro over 3 years ago

Hi.

L is the box side length.

In this tpr file is the end point of the simulation (acetonitrile molar fraction = 1), with 400 pre-equilibrated acetonitrile molecules, there is no water here, although in some other simulations issues also appear with water+acetonitrile mixtures that go away when switching from md-vv to md.

The periodic molecules are in it for no particular reason, I simply edited the mdp file from a previous simulation where I did need it (I had a surface), I guess I can remove this option (thanks for bringing it to my attention).

I am using virtual sites with the OPLS force field, but with GAFF there are no VS (itps from virtualchemistry.org) and I am also seeing strange behavior that goes away with md.

I will switch to a newer version, although with the large amount of scripts that would need modification, I think the best option for me at the moment is to simply switch to md for this project.

Many thanks for your help and the great work developing Gromacs.

#3 Updated by Mark Abraham over 3 years ago

Miguel Caro wrote:

Hi.

L is the box side length.

I presume you observed L to increase under constant-pressure situations? The above .tpr is constant-volume, however, so it's hard to do anything useful with it.

In this tpr file is the end point of the simulation (acetonitrile molar fraction = 1), with 400 pre-equilibrated acetonitrile molecules, there is no water here, although in some other simulations issues also appear with water+acetonitrile mixtures that go away when switching from md-vv to md.

The periodic molecules are in it for no particular reason, I simply edited the mdp file from a previous simulation where I did need it (I had a surface), I guess I can remove this option (thanks for bringing it to my attention).

Periodic molecules would be a prime suspect for an inadvertent bug - probably nobody has ever tested periodic + md-vv.

I am using virtual sites with the OPLS force field, but with GAFF there are no VS (itps from virtualchemistry.org) and I am also seeing strange behavior that goes away with md.

One has to use virtual sites to get a stable linear moiety, so whether you have one is a pertinent question - I don't know what various forcefields might have :-)

I will switch to a newer version, although with the large amount of scripts that would need modification, I think the best option for me at the moment is to simply switch to md for this project.

Understood.

Can we please have a set of your grompp inputs, e.g. mdp+gro+top for e.g. pure acetonitrile and some mix for which you found a problem. This will give us the necessary scope to tweak and test. Thanks!

#4 Updated by Miguel Caro over 3 years ago

Apologies for the late response.

A brief explanation of what I'm doing, in case it helps. I generate acetonitrile+water mixtures with different compositions, then I equilibrate at constant pressure to get the evolution of the density with compositions. Then, I fix the box size according to the equation I have fitted from the previous NPT simulations, and do constant volume equilibration. Then, I use the equilibrated structures and generate short 20ps runs (again NVT) which I post-process with my own implementation of the 2PT method. The file I gave you was for one of the latter runs.

When I calculate the entropy of the mixture, I can see strange behaviors and sharp discontinuities at some values. These disappear when I switch from md-vv to md. I also calculated the kinetic energy in each type of degree of freedom, and with md-vv there was too much kinetic energy in the vibrational degrees of freedom and too little in both translational and rotational degrees of freedom. This problem also goes away with md and is most likely directly related to the problem in the entropy.

Now I am attaching 3 sets of files (all with the GAFF as taken from virtualchemistry.org). One is pure acetonitrile (400 molecules), another is 370 acetonitrile + 60 water, and the other is 360 acetonitrile and 80 water. With md-vv I saw a sharp discontinuity in the entropy of the mixture when going from 370+60 to 360+80, so perhaps something fishy is going on there, as the box size increases. Use the files named "dynamics.gro" as starting points for the runs.

I have also done calculations, e.g. for dimethylformamide (DMF) with the same settings and the problem does not show. DMF has a density which is more similar to that of water, so I don't have to change the box size so much between simulations.

#5 Updated by Mark Abraham over 3 years ago

Thanks, that sounds great. I probably don't have time this week, and am on holiday next week, but when I'm back I'll prioritise seeing whether/what issue we might have still to fix!

#6 Updated by Michael Shirts over 3 years ago

Can I ask what happens when you increase the lincs order with md-vv? Since it works for some things and not for others, it would indicate that it may be something that changes from molecule to molecule, and the way the constraints behave is a clear potential example.

Alternatively, you can also try with constraints off (though I realize that's not quite the same problem, it would be interesting to see the results with md vs md-vv).

#7 Updated by Michael Shirts over 3 years ago

Michael Shirts wrote:

Can I ask what happens when you increase the lincs order with md-vv? Since it works for some things and not for others, it would indicate that it may be something that changes from molecule to molecule, and the way the constraints behave is a clear potential example.

Alternatively, you can also try with constraints off (though I realize that's not quite the same problem, it would be interesting to see the results with md vs md-vv).

And by increase, I would increase to lincs order 8.

#8 Updated by Miguel Caro over 3 years ago

Hi Michael, I also tried to run this simulations with fully flexible molecules, i.e. no constraints, and the issue still shows. It doesn't really show "for one molecule and not for another", since it also appears out-of-the-blue sometimes for the same set of molecules when the size of the box is changed. From my tests so far, box size and cutoffs are the quantities influencing this behavior. In a nut shell, increasing the box size or decreasing the cutoffs both have a negative effect.

#9 Updated by Michael Shirts over 3 years ago

Very strange. There should be very little difference between the two when running w/o constraints.

Could you post a new set of files that fails even with constraints off, including a README indicating which lines of .edr data show the problem? Remember, if we are not looking at the exact same files you are looking at, we're not going to be able to see the problem, which means we can't fix it.

Also, do you see any odd things when running in NVE?

#10 Updated by Michael Shirts almost 3 years ago

Did this get resolved?

Regardless, bumping so I remember to look at this when we are testing the physical consistency automated tests.

#11 Updated by Mark Abraham almost 2 years ago

Also available in: Atom PDF