GROMACS version: 2021.1
GROMACS modification: No
I’m trying to calculate the surface tension of a lipid bilayer I have simulated. I have a bilayer parallel to the xyplane, with a semiisotropic NP_zP_xyT ensemble. I have another simulation of the same setup with the NPgammaT ensemble as well. They produce different results, even though they are both set to be at zero tension, so I’m trying to figure out what the problem is by actually calculating the tension. (By different results, I mean that at 300K, the DMPC membrane I have in the semiisotropic case forms what looks like the ripple phase, while the bilayer in the NPgammaT case remains fluid. The latter is physically what I would expect.)
In the semiisotropic case, I have
pcoupl = ParrinelloRahman
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e5 4.5e5
ref_p = 1.0 1.0
which means that on average I should have P_xx=P_yy=1.0 bar, and P_z=1.0 bar, which should result in
gamma = 0.5L_z(P_zP_xx)=0.5L_z(P_zP_xx)=0
.
I figured I would use gmx energy
to spit out the pressure tensor and the box size so that I can calculate if the surface tension is in fact zero on average, for either simulation. At this point I encounter a number of confusing things. See the bottom of this post for the energy averages output. My questions are:

Averages of P_xx and P_yy are different from each other, although their error intervals seem to overlap. Is this expected?

P_xz is different from 0. Is this worrisome?

When I run
gmx energy
, there seems to be an option called#Surf*SurfTen
. What does this calculate, exactly? The equation I write above? If so, is it worrisome that my calculation ofgamma
from P_xx or P_yy doesn’t match this quantity? 
How is the error estimate in the output of
gmx energy
calculated? Is the number of independent snapshots somehow calculated (via some autocorrelation analysis, for example)? Or is there bootstrapping involved? 
Is it normal for the fluctuations of any element of the pressure tensor to be so large? For instance, P_xx fluctates between 600 and 600 bar.
The averages output by gmx energy
looks like this:
Statistics over 50000001 steps [ 0.0000 through 100000.0000 ps ], 15 data sets
All statistics are over 500001 points
Energy Average Err.Est. RMSD TotDrift

BoxX 9.52134 0.053 0.1252 0.371649 (nm)
BoxY 9.52134 0.053 0.1252 0.371649 (nm)
BoxZ 9.85503 0.096 0.232017 0.666881 (nm)
PresXX 0.430572 0.59 179.734 1.25861 (bar)
PresXY 0.589803 0.62 113.245 3.38764 (bar)
PresXZ 0.551734 0.14 121.067 0.114979 (bar)
PresYX 0.588489 0.62 113.243 3.39254 (bar)
PresYY 1.06575 0.65 179.605 0.89764 (bar)
PresYZ 0.114023 0.18 120.959 0.154326 (bar)
PresZX 0.54939 0.14 121.069 0.125734 (bar)
PresZY 0.108121 0.18 120.96 0.153392 (bar)
PresZZ 0.997789 0.19 203.787 0.414094 (bar)
#Surf*SurfTen 1.48211 1.9 2104.26 6.19725 (bar nm)
TMEMB 300 0.0015 2.33974 0.00786001 (K)
TSOLV 300.002 0.0012 1.46932 0.000656422 (K)