[BUG] MembThickness returns wrong value

[BlueSky2311]

When computing membrane thickness on a DPPC bilayer (140 lipids per leaflet), lipyphilic.analysis.memb_thickness.MembThickness produces unexpected results depending on the selected atom. For DPPC, the choline nitrogen (N) surface should lie further from the bilayer midplane than phosphate (P), and certainly further out than a tail carbon like C32. So, thickness(N) > thickness(P) > thickness(C32) is expected

However, I observed the opposite: C32 gives the largest thickness, and N is smaller than P.

I used latest version of lipyphilic (just install today), python version=3.13 (conda) on Linux.

This is the script to reproduce the issue:

import MDAnalysis as mda
import lipyphilic as lpp
from lipyphilic.analysis.memb_thickness import MembThickness
import numpy as np
from pandas import DataFrame
import argparse, os, sys

parser = argparse.ArgumentParser()
parser.add_argument("structure")
parser.add_argument("trajectory")
parser.add_argument("--sel", default="P")  # try: P, N, C32
parser.add_argument("--output", default="mem_thickness.csv")
args = parser.parse_args()

u = mda.Universe(args.structure, args.trajectory)
all_res = np.unique(u.residues.resnames)

lipids = []
for r in all_res:
    sel = u.select_atoms(f"resname {r} and name {args.sel}")
    if sel:
        lipids.append(r)

if not lipids:
    print(f"No residues with atom name '{args.sel}'")
    sys.exit(1)

lipid_sel = " or ".join([f"resname {r} and name {args.sel}" for r in lipids])

# Leaflets are assigned on the SAME selection used for thickness
leaflet = lpp.AssignLeaflets(universe=u, lipid_sel=lipid_sel)
leaflet.run(verbose=False)
filtered = leaflet.filter_leaflets(" or ".join([f'resname {r}' for r in lipids]))

mt = MembThickness(universe=u, leaflets=filtered, lipid_sel=lipid_sel)
mt.run(verbose=False)

thick = mt.memb_thickness
valid = ~np.isnan(thick)
DataFrame({"frame": np.arange(len(thick))[valid], "thickness": thick[valid]}).to_csv(args.output, index=False)

print(f"Mean: {thick[valid].mean():.2f} Å; Std: {thick[valid].std():.2f} Å; "
      f"Min: {thick[valid].min():.2f} Å; Max: {thick[valid].max():.2f} Å; "
      f"N: {valid.sum()}")