Interpretation of Pressure Tensor in log file

User discussions
1

GROMACS version: 2021.3
GROMACS modification: Yes
(version modified, but same question applies to 2021.3 without modification. Our modification should not alter the pressure calculation.)

I have a hopefully quick question about interpretation of the pressure tensor as printed on the .log file output by mdrun. Our particular system is a liquid crystal with a custom force field --our gromacs is modified to accept this alteration with significant direction from Hess (thank you!). The molecules of the LC are oriented along the z-axis of the simulation box and the simulation is run as pcoupltype = semiisotropic. Most of the simulation has been with the z-axis compressibility set to zero, but I’ve lifted this and run continuations where pcoupltype is set to isotropic and where it’s still semiisotropic but with the compressibility for z shifted from zero to 100-fold and 1000-fold smaller than x/y.

The oddity that I’m trying to understand is that when the pcoupltype is set to semiisotropic, the xx and yy elements of the pressure tensor in the log appear to have enforced opposite signs… after much equilibration, to where the volume of the simulation box is nearly unchanging, the magnitudes are close, but the signs are always opposite. And, I’ve seen this running a semiisotropic simulation of the same general sort in an unaltered version of gromacs as well. With the current simulation if I lift the forced semiisotropic type and make it fully isotropic, x and y become able to have the same sign. I thought that this might be a side effect of zeroing the compressibility along z, so I set z to be merely a small non-zero value and I still see the same effect.

This is troubling to me because the referencing value has no hope of being achieved in the x or y directions simultaneously and not really ever in z because the lack of compressibility there doesn’t allow it to change since x and y are canceling each other out. This is on the assumption that the system is trying to come to Pref supposing that P = (1/3)Trace(P-tensor).

Is there something here that I’m missing?

Thanks!
Greg

Quick Addendum: I’ve tried this with Berendsen and P-R barostats and seen the same in both cases. I am trying an experiment with anisotropic pcoupltype for completeness.

2

I’ve done some big work here and since I couldn’t find information about this all on a quick search prior to my original posting, I thought I would document my observations in better detail for the next time someone starts looking. I believe that I’m looking at a bug in pressure coupling for pcoupltype=semiisotropic. I don’t know if this is a bug that has been discussed elsewhere that I simply didn’t find information about, but on the off chance that the devs take interest, this is what I’m seeing. First, the system doesn’t appear to ever quite reach an equilibirum:

image
image736×456 68.9 KB

This plot is 1 ns window averages of a series of 10 ns runs I made with Berendsen barostat, varying tau-p trying to find a value that suppresses the fluctuations. Of note, this simulation is a continuation from a previous simulation that was at 450 ns length, where the density of the system is in an apparent equilibrium. Here are the pressure averages for 20 ns intervals on that 450 ns simulation, the error bars are the calculated error on the 20 ns blocks. The reference pressure was 1 bar:

image

There are two separate simulations in this plot, one as orange and one as blue; the details of these are immaterial. As you can see, the pressure seems to wander up and down and the error only rarely encompasses the reference pressure. Further, in my simulations, the pressure tensor is always showing this bizarre feature where the x and y dimensions are restricted to have opposite signs:

Pressure (bar)
-1.50627e+01 4.89732e+01 -7.04428e+00
4.89732e+01 1.81663e+01 6.42916e+00
-7.04427e+00 6.42916e+00 -9.72516e+01

And, I as I had run for a long time with the z-axis compressibility set to zero, I tested and found that it is independent of that setting, still showing the same feature regardless of whether the z-axis is turned off or not.

I concluded that the bug is in the semiisotropic pcoupltype because it is insensitive to barostat choice, reference choice and time constant choice. I was able to run my system as pcoupltype=anisotropic with the compressibilities set to mimic a semiisotropic run and the pressure became much more stable:

image

This may not look like much, but it’s a continuation across the same time region as the very first plot, again with the same 1 ns window average, but now anisotropic set to mimic semiisotropic. The average across the whole 10 ns region is 6.7 bar with an error of 7.5 bar where the reference value was 1.0 bar. I can probably tweak the time constant response, but this is basically right on track now in only 10 ns. Moreover, the pressure tensor, while clearly not totally at equilibrium, releases the x and y values so that they can have the same sign:

Pressure (bar)
5.29416e-01 4.75237e+01 -1.69066e+01
4.75237e+01 5.49666e-01 -4.25935e+00
-1.69065e+01 -4.25937e+00 1.90526e+01

For anybody else who has my problem, this was my solution, pending hope that the devs get a chance to do something with the semiisotropic method to correct it in some upcoming version of Gromacs.

Greg

3

Hi – I’m curious if there has been any follow-up about semi-isotropic pressure coupling. Though I am using an older version of Gromacs sitll (2018.3), I have had issues getting semi-isotropic pressure coupling to work. I see papers published that have used it with this same version of gromacs, and so I’m curious if it has been in fact validated and/or if there are known bugs and workarounds.

Also, if there is a known bug, how do you use anisotropic coupling to mimic semisotropic coupling? I would still like the reference pressures to be the same in all three directions even though I want semi-isotropic coupling (I’m doing a membrane simulation), so I’m not sure how to “fool” gromacs to mimic semi-isotropic via an anisotropic setting.

Thanks for any information!!

Mala

4

Hi Mala,

I’ve been continuing to work on my system. I’ve had some experience more recently that suggests that this can be involved with how the version of gromacs was compiled. I typically use some local HPC cluster resources and the resources allocated to a simulation job can be quite flexible; I’ve found that I’m getting more accurate pressure coupling results if the resources allocated to run a particular job are the same as those where the gromacs version was compiled. This has seemed quite subtle to me because the cluster system where I’ve been running is effectively housing virtual machines on the nodes; meaning that the node with 64 cpus and 3 gpus can actually house multiple virtual machines, say 10 cpus and 1 gpu. Gromacs can seemingly run if the allocated resources are different from those used when it was compiled, but I’m getting way better results if they match.

My understanding is that gromacs has been altered quite extensively since the 2018 version. My running has typically been with 2020 through 2023 versions (though I’ve been mainly using 2022 recently due to some compilation issues with 2023). The pressure coupling and barostat systems between 2022 and 2018 may not have the same issues.

In my experience, the semi-isotropic and anisotropic settings are somewhat similiar-ish in how they function, but anisotropic allows all the dimensions to vary with respect to one another while semi-isotropic is keeping x-y coupled together and varying x-y as a unit independently from z. Often times, I am fixing one dimension of the simulation box to a specific size by setting the compressibility for that direction to zero and then relying on the other compressibility values to vary their dimensions in the barostat while they chase the reference pressure. I’ve noticed that when I want it to work the best, it’s important to be using the Parinello-Rahman barostat. Some of the other barostats can be really useful for achieving a desired density quickly, but they seem to have weird effects if you continue to use them after you’re at an equilibrium density. I’ve seen C-rescale super-compress a material several times and claim it’s near the reference pressure.

For a system like a membrane, I would put the long axis of the lipid along the z-axis and have the membrane stretch x-y (which is probably what you’re doing). If you’re compressing only along the z-axis for the simulation, setting x and y compressibilities to zero would make the anisotropic version the same as the semiisotropic version if the x-y setting is set to zero; this would allow it to compress from z but not permit any expansion or contraction in x and or y. If you’re trying to bunch the lipids together by compressing in x-y and keeping z fixed, anisotropic and semiisotropic would be similar, but the anisotropic has the capacity to deform the box aspect ratio where x can be quite different from y. It doesn’t seem to run around much if you’re near equilibrium, but it can and does wander around some.

Those are my experiences, anyway.

Greg