GROMACS

Thanks for the report. I've started GPU and CPU runs to see what I reproduce.

What is the source of your membrane parameters and protocol? Is that known to be stable ie. not undergo phase transitions?

Can you please let us also have the full topology file and supporting files, so we can build our own .tpr if we need to?

How long did you run a CPU run? Also more than 100 ns?
Membranes are/were not stable at constant surface tension in some version of the CHARMM force field, but I forgot with one(s) that was/were.

Issue 2007 affected this configuration, but that was fixed in version 5.1.4, in a patch release of 2016 and is not present in the 2018 release. But maybe this was not fixed correctly.

Can you attach a md.log file of a 2018 GPU run that produced incorrect results?

Berk Hess wrote:

How long did you run a CPU run? Also more than 100 ns?
Membranes are/were not stable at constant surface tension in some version of the CHARMM force field, but I forgot with one(s) that was/were.

Dear Berk,

- The CPU run was run for 500ns
- I have been using CHARMM with the same parameters for more than 6 years and never came across any problem like this. Since I use the CHARMM-GUI interface to generate my systems, the version I am using is always the last stable release (i.e. in this case, C36m).

Berk Hess wrote:

Issue 2007 affected this configuration, but that was fixed in version 5.1.4, in a patch release of 2016 and is not present in the 2018 release. But maybe this was not fixed correctly.

Unless the issue that caused #2007 was not correctly identified, to me it seems that the fix on 5.1 was appropriate.

Gerrit received a related patchset '1' for Issue #2845.
Uploader: Berk Hess (hess@kth.se)
Change-Id: gromacs~release-2019~Ib5c250a1f370d4c0a4fd652bf6efa03e70e18748
Gerrit URL: https://gerrit.gromacs.org/9066

I did actually test myself that this fix significantly reduces the difference in pressure between CPU and GPU at step 0.

Gerrit received a related patchset '1' for Issue #2845.
Uploader: Paul Bauer (paul.bauer.q@gmail.com)
Change-Id: gromacs~release-2016~Ib5c250a1f370d4c0a4fd652bf6efa03e70e18748
Gerrit URL: https://gerrit.gromacs.org/9083

Gerrit received a related patchset '1' for Issue #2845.
Uploader: Paul Bauer (paul.bauer.q@gmail.com)
Change-Id: gromacs~release-2018~Ib5c250a1f370d4c0a4fd652bf6efa03e70e18748
Gerrit URL: https://gerrit.gromacs.org/9084

Applied in changeset de27733a387922033b116bf7024c595e5b1e99a8.

Applied in changeset c131d790b48da10f5fcb9e4aabeef0d730b76724.

Could you attach a box plot for a long cpu run, so we can see what you think is normal and abnormal?

Have you continued the 2016.6 GPU run? It could be that the system with protein is just on the limit of the stability with the force field. So maybe the 2016.6 GPU results are just "unlucky", but the system remains at the same box size for most of the time?

Are you 100% sure that you are running 2016.6 with the fix?
I checked that the fix corrected the energy and the pressure, so the fix actually fixes a bug that occurred with your run setup. I can not exclude that there is a second bug, but I have not clue where that could be. The only special code with force-switching on GPUs is the few lines where the fix is.

Could you attach the md.log file for the GPU 2016.6 run?

Ramon Guixa-Gonzalez wrote:

Berk Hess wrote:

Are you 100% sure that you are running 2016.6 with the fix?
I checked that the fix corrected the energy and the pressure, so the fix actually fixes a bug that occurred with your run setup. I can not exclude that there is a second bug, but I have not clue where that could be. The only special code with force-switching on GPUs is the few lines where the fix is.

Could you attach the md.log file for the GPU 2016.6 run?

Well, I am running the 2016.6 version recently released (http://manual.gromacs.org/documentation/2016.6/download.html). There, you state this is fixed right? (http://manual.gromacs.org/documentation/2016.6/ReleaseNotes/index.html).

Please find attached the md.log file for the GPU 2016.6 run.

Hi there, have you figured what the problem is here? I only have access to a GPU cluster, so since this issue came up I am basically stuck...

This is an issue with high priority, as the setting says. But I have little clue where to look. Normally I would test things myself. But since this happens after so long time, and not without the protein, that is impractical.

The only suggestion I have is to run with a plain cut-off instead. Using a plain cut-off of 1.11 nm gives the same potential and force at distances between 0 and 1 nm and only differs in the switching region. This should not significantly change the properties of the system. Whether this is stable or not should also give us hints about this issue.

PS I meant a plain cut-off with shifted potential, which is the default potential modifier with the Verlet cutoff scheme.

Berk Hess wrote:

Do you already have results with a 1.11 nm cut-off without switch?

Sorry, I have not been able to launch this test yet.
Would you mind confirming the lines I need to have in my mdp in order to perform this run?

This is what I currently use:
cutoff-scheme = Verlet
coulombtype = pme
rcoulomb = 1.2

Should I just use this:
cutoff-scheme = Verlet
coulombtype = Cut-off
rcoulomb = 1.1

Would this be all?

Thanks

R

No, still use PME.
Change:
rcoulomb=1.11
rvdw=1.11
vdw-modifier=potential-shift (which is the default)

This should also run faster.

You seem to be running with the default PME fourier spacing. You could have used a coarser grid with rcoulomb=1.2. But I would leave it as is with the 1.11 nm cut-off.

Could you attach the log file as well?

Still no clue. Now we only know that it has nothing to do with the force-switch functionality.

Ramon Guixa-Gonzalez wrote:

Berk Hess wrote:

Could you attach the log file as well?

Sure, here you go Berk

Paul Bauer wrote:

@Berk and @Ramon, any progress with identifying what is happening there?
Bumped for now to next point release.

Berk Hess wrote:

Still no clue. Now we only know that it has nothing to do with the force-switch functionality.

Same thing. I don't know what else to check here...

I am (still) wondering as to why this happens only after 300 ns. Are there slow rearrangements in the system?
Can you restart the simulation from a checkpoint file or trajectory frame close to 300 ns to reproduce the issue much quicker and also check what happens with a CPU run? If so, we might be able to find what the issue is.

Hi there, ok, I will do that. Will keep you posted

I have seen similar simulation box deformations for POPC + protein simulations with Charmm36 FF. These deformations are observed using Gromacs 2018.6 in CPU version (I observe similar deformations for a POPC + CHOL + Protein MD runs in both CPU and GPU versions, while one repeat on CPU does not show deformations after 1000 ns with 15 restarts, two other CPU repeats show deformation as early as 220 ns. One GPU run shows deformations around 400 ns).

For the POPC + Protein system, in two CPU repeats I have noticed that these deformations (deformations begin near 600 ns or 1000 ns) with runs restarted every 24 hours (the maximum walltime available in one of the clusters I have access, approx. 18 restarts for 1000 ns). For a single GPU repeat of the same system with 96 hours max. walltime, I haven't seen these deformations after 1000 ns (5 restarts).

I don't know if this helps in any way, since sometimes I get lucky with some runs.

So it doesn't seem to be clear yet that there is an actual code issue here. It could be that the force field has some unwanted state that is close in free-energy to the intended state.
One question now is what the statistical certainty is of this happening in GPU and not is CPU, or even the other way around now?

As I said before, the most useful would be to have a checkpoint close before a fluctuation occurs.

Isn't this just the result of having a too short tau_P (5 ps) with semiisotropic scaling? Since relaxation is system size dependent the "common knowledge" that a 5 ps time constant is long enough may not hold. Try increasing it to 20 ps.

I missed that you had uploaded checkpoint files two weeks ago. I have reproduced the issue and have been running many simulations to find leads to what is wrong, but didn't find any hint up till now.

Do you know if the issue seems related to running non-bonded interactions or PME on the GPU?

David van der Spoel wrote:

Isn't this just the result of having a too short tau_P (5 ps) with semiisotropic scaling? Since relaxation is system size dependent the "common knowledge" that a 5 ps time constant is long enough may not hold. Try increasing it to 20 ps.

sorry, I was away, I'll give long tau_Ps a try, thanks!

Berk Hess wrote:

I missed that you had uploaded checkpoint files two weeks ago. I have reproduced the issue and have been running many simulations to find leads to what is wrong, but didn't find any hint up till now.

Do you know if the issue seems related to running non-bonded interactions or PME on the GPU?

Hey Berk, well, not a clue, really...

Berk Hess wrote:

I ran a GPU simulation from the checkpoint you provided and started runs from a later checkpoint at time 179558.6 ps with box_z=9.3686.
The plot shows that with GPU runs box_z fluctuates significantly and has spikes downward. But also, the CPU run, which runs much slower, shows a drop to 8.5. So this is likely not a GPU issue.
My hypothesis is still that the force field has a second free-energy minimum at state compressed in z and expanded along x/y.

Ok...

Berk Hess wrote:

My CPU continuation runs has dips down to box_z=8.2 nm.
So our conclusion is that this is a force field issue (which only appears after long simulation times, so likely slow rearrangements are involved finally leading to rapid changes).

Well, ok, this was also what I initially suspected. I will try to report it to the CHARMM force field community to see if they are aware. On the meanwhile, I will run more tests, namely:
1) one of the failing simulations on NAMD
2) one at longer tau_P (20)
3) one using an older release of the CHARMM force field

Thanks Berk

I can confirm this bug with GROMACS 2019.3. I have several simulations of a GPCR+Membrane system that get ripped apart after 100-200ns. But already before that unnaturally large fluctuations are seen that manage to recover however. After some time the system violently extends in x-y dimension (which is the membrane plane) and collapses in the z-dimension. I have two simulations of identical coordinates that were started with different random seeds but otherwise identical parameters and both get ripped apart. Interestingly, both systems have very similar box dimensions after the "event" although the disrupted membrane structure is very different between the two. Since I didn't notice the error immediately and the simulations did not crash, they continued for several hundert ns (first simulation for 450ns, second simulation for 900 ns). There still seem to be some pretty large fluctuations during that time but no violent changes as before.

I'll start a third version of the system today. A colleague of mine reported similar issues but I didn't have time to look at it yet.

I setup the simulation with charmm-gui and run them on an AMD 1700X + Nvidia RTX2080Ti. If more info is needed please tell me what you need.

Alexander Vogel