[QUESTION]Units of diffusion coefficient

[ErikPoppleton]

What is the topic of your question: Usage / Documentation, maybe Bug

Add your question below: When I run the LyPyphilic lateral_diffusion.MSD module on simulations of Amber Lipid21 lipids in GROMACS and then calculate the diffusion coefficient, I get an answer that is the wrong order of magnitude for standard lipids. The following code is based off the example notebook.

import matplotlib.pyplot as plt
from lipyphilic.lib import lateral_diffusion
import MDAnalysis as mda
from MDAnalysis import transformations as mda_transformations
from lipyphilic import transformations as lpy_transformations

def get_diffusion(u, sel_str):
  msd = lateral_diffusion.MSD(
    universe=u,
    lipid_sel=sel_str,
    com_removal_sel=sel_str,
    dt=0.2
  )

  msd.run(verbose=True)
  return msd

# Load DOPC
u_DO = mda.Universe('../big_DOPC/md.tpr', '../big_DOPC/md.xtc')
lipids_DO = u_DO.select_atoms("moltype OOPC")
t_DO = [lpy_transformations.triclinic_to_orthorhombic(lipids_DO), mda_transformations.NoJump(lipids_DO)]
u_DO.trajectory.add_transformations(*t_DO)

# Load DLPC
u_DL = mda.Universe('../DLPC/md.tpr', '../DLPC/md.xtc')
lipids_DL = u_DL.select_atoms("moltype LLPC")
t_DL = [lpy_transformations.triclinic_to_orthorhombic(lipids_DL), mda_transformations.NoJump(lipids_DL)]
u_DL.trajectory.add_transformations(*t_DL)

# Compute msd
msd_DO = get_diffusion(u_DO, 'moltype OOPC')
msd_DL = get_diffusion(u_DL, 'moltype LLPC')
mean_msd_DL = np.mean(msd_DL.msd, axis=0)
mean_msd_DO = np.mean(msd_DO.msd, axis=0)

# Make plot of msd/time
fig, ax = plt.subplots()
ax.loglog(msd_DL.lagtimes, mean_msd_DL, label="DLPC")
ax.loglog(msd_DO.lagtimes, mean_msd_DO, label='DOPC')
ax.set_xlabel("Lagtime (ns)")
ax.set_ylabel(r"MSD$_{xy}\ \rm{(nm^2)}$")
plt.legend()
plt.show()

# Calculate diffusion coefficient
d_DL, sem_DL = msd_DL.diffusion_coefficient()
d_DO, sem_DO = msd_DO.diffusion_coefficient()
print("lipid\tDiff\tstdev")
print(f"DLPC {d_DL:.2e} {sem_DL:.2e}")
print(f"DOPC {d_DO:.2e} {sem_DO:.2e}")

The resulting graph shows really high MSD:

And the diffusion coefficients are ~10000x too large

lipid Diff stdev
DLPC 6.11e-04 6.09e-05
DOPC 2.40e-03 3.60e-05

Lipid diffusivity should be on the order of 1e-8 cm^2/sec, and Amber Lipid14 did capture this (table 9) (there isn't diffusion data in the Lipid21 paper)

I have tried messing with the dt argument of MSD. My simulation uses a tilmestep of 0.002 ps (2 fs) and saves a trajectory frame every 100000 steps (200 ps = 0.2 ns). Therefore I wouldn't expect that setting dt to 0.2 would change anything, but setting it does change the calculated diffusion coefficient (2 orders of magnitude higher when set). The documentation says:

dt (float, optional) – The time, in nanoseconds, between consecutive frames in universe.trajectory. The defualt is None, in which case dt is taken to be universe.trajectory.dt divided by 1000.

In my code, u_DO.trajectory.dt is 200, so if I were to not set dt, the result should also be 0.2.

Is there something I am not understanding, or is there maybe a unit conversion error somewhere in the code?