Interpretation of Pressure Tensor in log file

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.

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:

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