CHARMM Development Project Forums User Discussion & Questions Setup, I/O, and Basic questions Lipid diffusion with IMSD?

Hi,

i was wondering if I can use IMSD to calculate the lipid diffusion or if the algorithm is based on some assumptions for water (All I found was IMSD in the context with water).



Cheers and thanks
Bjoern

I think you mean to ask if you can use COOR ANAL to calculate the lipid diffusion.
There is no assumption about water, but you have to be aware that the MSD calculation is made for each individual atom in the selection (no center-of-mass or other positional averaging is done). It is not clear to me why you don'tn need PBC (support for which is limited in COOR ANAL).

If the lipids were subject to image centering during the simulation, then you must account for that; but that's not the biggest problem.

I don't believe COOR ANAL is suitable for lipid diffusion. The rigorous approach to compute diffusion is to use the center-of-mass (COM), and for lipids, that's often done as a 2D diffusion, not 3D. COOR ANAL is not set up to use the COM to compute r^2 vs t, or to compute r in 2D instead of 3D.

Also, coupling between neighbors leads to size effects; see

Dynamical motions of lipids and a finite size effect in simulations of bilayers
Jeffery B. Klauda, Bernard R. Brooks, and Richard W. Pastor
J. CHEMICAL PHYSICS 125(14) Article No. 144710 (2006)

See also this thread

Note that the finite size correction discussed in that thread only applies to cubic systems at present, and would not help with the coupling problem.

Somewhere I read that in an constant pressure simulation the box size is taken from the trajectory and there is no need to add "xbox @x ybox @y zbox @z"

The simulation itself was run with PBC.




I see. So I have to come up with an own code to get the proper diffusion rate comparable to lab experiments.

But lets make the assumption I only want to observe, if the lipid diffusion constant changes during the simulation. Selection a single atom from the headgroup and calculating the diffusion constant for this atom - could this reflect those changes, couldn't it? The quick and dirty way to see if all the afford of coming up with a code to calculate the proper diffusion constant pays off.

Thanks for the link and the paper.

Cheers
Bjoern

If the lipids were subject to image centering during the simulation, then you must account for that. Otherwise, the r^2 vs. t data will be contaminated by large jumps across the unit cell.

In order to have the box size taken from the stored trajectory frames, you must both set up CRYSTAL and use XBOX etc.; the lattice must be orthogonal (all angles 90 deg) which means only CUBIC, TETRAGONAL, or ORTHORHOMBIC lattice types. (The same restriction applies to MERGE COOR UNFOLD; only orthogonal lattice types are supported.)

Calculating diffusion accurately for just one molecule is highly problematic, regardless of the approach; standard procedure is to average over molecules and/or replicate simulations. Lipids pose bigger challenges due to the presence of short time diffusion based on "rattling in a cage", and long time diffusion based a jump model with lipids changing places. Several replicate simulations might be needed to get statistically meaningful results for a small number of lipids.


While I haven't had time to clean up them up for the Script Archive, I do have CHARMM inputs, csh batch scripts, and a Fortran90 program which do these steps:

(1) image unfold the coordinate trajectory to a scratch volume (MERGE)
(2) for each molecule, compute the r^2 vs. t difference correlation function (CORREL), in both 3D and 2D
(3) compute D from input r^2 vs t data, in 3D or 2D, with optional finite size correction (cubic only); provides standard errors based on grouping (F90 prog)

Available via email request.

I guess this is best done vie unwrapping the trajectory. Or is there any build in code which is able to detect those jumps. The only way, except unwrapping, I can imagine right now is to check frame for frame if a large translation occurs. But this would be rather time consuming.


That's fine for me, my lattice is orthogonal.


Mhh, well my system consists of 86 POPC, 43 POPS, and 43 Cholesterols per leaflet so I guess this counts as small number.




I guess using and modifying your code yields faster a result then writing my own and since you are willing to share it I don't have to re-invent the wheel. Thanks a lot! Mail has been dispatched.

Thanks a lot!
Kind regards
Bjoern

COOR ANAL can handle PBC with orthogonal lattice types, but only computes MSD for single atoms, not the molecular COM.

By "small", I meant less than ten or so lipids; your system seems okay.

Note that since there is both a short time and long time diffusion, one must care in the fitting to discard either the short time or long time data. Plots of r^2 vs t, esp. averaged over molecules, should make it clear where the slope change occurs.

I see.



I am interested in the long time diffusion so I'll drop the short time data. Mhhh, well thinking twice about it, I am probably interested in both parts so I have to calculate it twice Also i will cut my simulation into some larger parts and calculate the diffusion only for these parts to compare them against each other.

Thanks for the notes.

Kind regards
Bjoern

Diffusion on both timescales has been measured experimentally.

You should compare the two leaflets of the bilayer as well.

Already on my ToDo list.
It is one of the most interesting part since I have a membrane binding protein on the lower leaflet.

Thanks again for your help!
Cheers
Bjoern