CHARMM Development Project
Posted By: blubbi Lipid diffusion with IMSD? - 10/08/10 04:16 PM
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).

Code:
! All POPC lipids (check with MSD vs. one Atome selection)
define LOI sele RESN POPC .and. type P* end

coor anal sele LOI end -
firstu @fU nunit @N begin @F skip @SK stop @L -
imsd 21 rspin 0.0 rspout 999.9 ncors @nsteps

!NPgT simulation no PBC setup required


Cheers and thanks
Bjoern
Posted By: lennart Re: Lipid diffusion with IMSD? - 10/08/10 05:11 PM
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).
Posted By: rmv Re: Lipid diffusion with IMSD? - 10/08/10 06:37 PM
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.
Posted By: blubbi Re: Lipid diffusion with IMSD? - 10/09/10 10:49 AM
Originally Posted By: lennart
It is not clear to me why you don't need PBC (support for which is limited in COOR ANAL).


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.


Originally Posted By: rmv
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.


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
Posted By: rmv Re: Lipid diffusion with IMSD? - 10/09/10 04:56 PM
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.
Posted By: blubbi Re: Lipid diffusion with IMSD? - 10/09/10 05:54 PM
Originally Posted By: rmv
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.

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.

Originally Posted By: rmv
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.)

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

Originally Posted By: rmv
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.

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


Originally Posted By: rmv
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 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
Posted By: rmv Re: Lipid diffusion with IMSD? - 10/12/10 05:40 PM
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.
Posted By: blubbi Re: Lipid diffusion with IMSD? - 10/14/10 08:08 AM
Originally Posted By: rmv
COOR ANAL can handle PBC with orthogonal lattice types, but only computes MSD for single atoms, not the molecular COM.

I see.

Originally Posted By: rmv

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 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 wink 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
Posted By: rmv Re: Lipid diffusion with IMSD? - 10/14/10 08:33 PM
Diffusion on both timescales has been measured experimentally.

You should compare the two leaflets of the bilayer as well.
Posted By: blubbi Re: Lipid diffusion with IMSD? - 10/15/10 01:15 PM
Originally Posted By: rmv
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
Posted By: blubbi Re: Lipid diffusion with IMSD? - 11/19/10 02:16 PM
RMV, thanks again for your support and Fortran code (which I have successfully ported to python with the help of the cited papers [Fortran code is not that straight forward when you are used to python :-) ]).

Still I am wondering about one thing.
The long time diffusion should be calculated from a region where the r2com is more or less linear.

From the attached picture (r2com for 86 POPC) I would tend to use frames between ~6000 and ~13000.

Is this reasonable? I am somewhat unsure if it is okay to discard so much points.

Currently I am calculating r2com for different time windows to get a clue how this effects the shape of the curve.

Kind regards,
Bjoern

Attached picture r2com.png
Posted By: rmv Re: Lipid diffusion with IMSD? - 11/19/10 05:41 PM
I tend to use the averaged data to make a judgment about the slope change between the two diffusion time scales for lipids in a bilayer. It's not always clear from the r^2 vs t plots for the individual lipids. Also, try plotting just 0-5000 ps to look for the slope change.
© CHARMM forums