
Joined: Feb 2009
Posts: 83
Forum Member

OP
Forum Member
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). ! 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,861 Likes: 10
Forum Member

Forum Member
Joined: Sep 2003
Posts: 4,861 Likes: 10 
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 centerofmass 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,623 Likes: 24
Forum Member

Forum Member
Joined: Sep 2003
Posts: 8,623 Likes: 24 
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 centerofmass (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 threadNote 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
Forum Member

OP
Forum Member
Joined: Feb 2009
Posts: 83 
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. 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 centerofmass (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 threadNote 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,623 Likes: 24
Forum Member

Forum Member
Joined: Sep 2003
Posts: 8,623 Likes: 24 
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
Forum Member

OP
Forum Member
Joined: Feb 2009
Posts: 83 
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. 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. 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. 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 reinvent the wheel. Thanks a lot! Mail has been dispatched. Thanks a lot! Kind regards Bjoern




Joined: Sep 2003
Posts: 8,623 Likes: 24
Forum Member

Forum Member
Joined: Sep 2003
Posts: 8,623 Likes: 24 
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
Forum Member

OP
Forum Member
Joined: Feb 2009
Posts: 83 
COOR ANAL can handle PBC with orthogonal lattice types, but only computes MSD for single atoms, not the molecular COM. I see. 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 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,623 Likes: 24
Forum Member

Forum Member
Joined: Sep 2003
Posts: 8,623 Likes: 24 
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
Forum Member

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




