Calculating surface tension of lipid bilayer using gmx energy

User discussions
1

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 xy-plane, 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                  = Parrinello-Rahman
pcoupltype              = semiisotropic
tau_p                   = 5.0
compressibility         = 4.5e-5  4.5e-5
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_z-P_xx)=0.5L_z(P_z-P_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:

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

  2. P_xz is different from 0. Is this worrisome?

  3. 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 of gamma from P_xx or P_yy doesn’t match this quantity?

  4. 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?

  5. 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  Tot-Drift
-------------------------------------------------------------------------------
Box-X                       9.52134      0.053     0.1252  -0.371649  (nm)
Box-Y                       9.52134      0.053     0.1252  -0.371649  (nm)
Box-Z                       9.85503      0.096   0.232017   0.666881  (nm)
Pres-XX                    0.430572       0.59    179.734   -1.25861  (bar)
Pres-XY                   -0.589803       0.62    113.245   -3.38764  (bar)
Pres-XZ                    0.551734       0.14    121.067  -0.114979  (bar)
Pres-YX                   -0.588489       0.62    113.243   -3.39254  (bar)
Pres-YY                     1.06575       0.65    179.605    0.89764  (bar)
Pres-YZ                    0.114023       0.18    120.959   0.154326  (bar)
Pres-ZX                     0.54939       0.14    121.069  -0.125734  (bar)
Pres-ZY                    0.108121       0.18     120.96   0.153392  (bar)
Pres-ZZ                    0.997789       0.19    203.787   0.414094  (bar)
#Surf*SurfTen               1.48211        1.9    2104.26    6.19725  (bar nm)
T-MEMB                          300     0.0015    2.33974 -0.00786001  (K)
T-SOLV                      300.002     0.0012    1.46932 0.000656422  (K)
2

Pressures fluctuate a lot, as you see. In you formula you should use the average of the xx and yy components, not only zz. This is what the #Surf*SurfTen term does, but it doesn’t have the factor 0.5 as it doesn’t know how many interfaces there are.

gmx energy -h states that the error estimate is computed from the variance of block averages using 5 block (by default).

PS: note that the reported surface tension is zero withing the error estimate.

3

Thank you so much for the clarification.

I am now calculating the following for each frame:

2*gamma = L_z*(P_z-0.5(P_xx+P_yy))

and comparing it with #Surf*SurfTen value for each frame. I find that the standard deviation of #Surf*SurfTen is an order of magnitude smaller than that of the 2*gamma value I calculate. Further, the means are also off by almost a factor of 0.4. What could I be missing?

4

Usually energies are computed much more often, every nstcalcenergy steps, then frames are saved, every nstenergy steps. You should use the average values from the energy file. You can use the total average reported by gmx energy or get the average between frames using the -aver option.