Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
#25574 10/08/10 04:16 PM
Joined: Feb 2009
Posts: 83
B
blubbi Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Feb 2009
Posts: 83
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

Joined: Sep 2003
Posts: 4,840
Likes: 3
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 4,840
Likes: 3
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).


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Joined: Sep 2003
Posts: 8,597
Likes: 12
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,597
Likes: 12
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.


Rick Venable
computational chemist

Joined: Feb 2009
Posts: 83
B
blubbi Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Feb 2009
Posts: 83
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

Joined: Sep 2003
Posts: 8,597
Likes: 12
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,597
Likes: 12
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.


Rick Venable
computational chemist

Joined: Feb 2009
Posts: 83
B
blubbi Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Feb 2009
Posts: 83
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

Joined: Sep 2003
Posts: 8,597
Likes: 12
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,597
Likes: 12
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.


Rick Venable
computational chemist

Joined: Feb 2009
Posts: 83
B
blubbi Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Feb 2009
Posts: 83
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

Last edited by blubbi; 10/14/10 08:11 AM.
Joined: Sep 2003
Posts: 8,597
Likes: 12
rmv Offline
Forum Member
Offline
Forum Member
Joined: Sep 2003
Posts: 8,597
Likes: 12
Diffusion on both timescales has been measured experimentally.

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


Rick Venable
computational chemist

Joined: Feb 2009
Posts: 83
B
blubbi Offline OP
Forum Member
OP Offline
Forum Member
B
Joined: Feb 2009
Posts: 83
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

Page 1 of 2 1 2

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.017s Queries: 35 (0.011s) Memory: 0.7892 MB (Peak: 0.8889 MB) Data Comp: Off Server Time: 2022-07-03 12:41:54 UTC
Valid HTML 5 and Valid CSS