Regarding membrane curvature analysis of bilayer

[SoumyaCYZ]

@ojeda-e @hmacdope @orbeckst Ma'am/sir, I have calculated the 2D map of the membrane of curvature using the PO4 index from my CG system. It gives me a nice contour plot. What is the unit of the x- and y-axis in the contour plot? Suppose my membrane's surface area is 25 * 25 nm2. Is the range of x-axis is 0 to 25 and y is 0 to 25? It'll be a nice help if you clear these doubts.

[orbeckst]

Can you show the code that you used?

[SoumyaCYZ]

I'm attaching here my code

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from membrane_curvature.tests.datafiles import (MEMB_GRO, 
                                                MEMB_XTC)

u = mda.Universe("mem_cg.gro", "mem_cg.xtc")

mc_upper = MembraneCurvature(u, 
                             select="resid 1-960 and name PO4",
                             n_x_bins=12, 
                             n_y_bins=12).run()

mc_lower = MembraneCurvature(u, 
                             select='resid 961-1920 and name PO4', 
                             n_x_bins=12, 
                             n_y_bins=12).run()
# surface
surf_upper = mc_upper.results.average_z_surface
surf_lower = mc_lower.results.average_z_surface

# mean curvature
mean_upper = mc_upper.results.average_mean
mean_lower = mc_lower.results.average_mean

# Gaussian curvature
gauss_upper = mc_upper.results.average_gaussian
gauss_lower = mc_lower.results.average_gaussian
import matplotlib.pyplot as plt
from scipy import ndimage
import numpy as np

results = [mean_lower, mean_upper]
leaflets = ["Lower", "Upper"]

fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(4,3.5))
max_ = max([np.max(abs(h)) for h in results])

for ax, rs, lf in zip((ax1, ax2), results, leaflets):
    rs = ndimage.zoom(rs, 4, mode='wrap', order=1)
    levs = np.linspace(-max_, max_, 40)
    im = ax.contourf(rs, cmap='bwr', 
                     origin='lower', levels=levs,
                     vmin=-max_, vmax=max_)

    ax.set_aspect('equal')
    ax.set_title('{} Leaflet'.format(lf), fontsize=6)
    ax.axis('off')

cbar = plt.colorbar(im, ticks=[-max_, 0, max_], 
                    orientation='vertical', 
                    ax=[ax1, ax2], shrink=0.4, 
                    aspect=10, pad=0.05)
cbar.ax.tick_params(labelsize=4, width=0.5)
cbar.set_label('H (nm$^{-1}$)', fontsize=6, labelpad=2) 

[orbeckst]

It looks as if membrane-curvature does not explicitly keep track of coordinates and just uses the bins. Thus, if you use n_x_bins = n_y_bins = 12 then you have to look at your system x-y dimensions to calculate the size and extent of each bin. If you're simulating with constant area then you can just get the dimensions LX and LY from LX, LY, LZ, alpha, beta, gamma = u.dimensions at the first step; then one bin in X is LX/n_x_bins wide. If you are simulating with constant pressure you probably should average LX and LY over the whole trajectory and take the average. LX and LY are in Ångstrom (which is the standard length unit in MDA).

Once you know the actual bin sizes I would calculate the midpoints for each bin and then use np.meshgrid to generate the proper X, Y coordinate 2D arrays as input for contourf.

I don't know what ndimage.zoom() does and what rs contains. That's a part that you have to investigate yourself.

[SoumyaCYZ]

@orbeckst sir, can you please respond?

[orbeckst]

Btw, if you mark something as an answer then you indicate that this answered your question. If you think that you answered your own question then mark that as an answer; however, given that you're asking for a follow-up here I assume that you didn't actually think that you answered your question so I suggest you unmark your "answer" so that others aren't confused.