#6857 - 05/25/05 09:36 PM
lipid C-D order parameters using CORREL
|
Forum Member
Registered: 09/17/03
Posts: 8406
Loc: 39 03 48 N, 77 06 54 W
|
The following 3 scripts compute the carbon-deuterium order parameter, S(C-D), for all of the aliphatic C-H bond vectors in the lipid DPPC. The main assumption is that the lipids are oriented, i.e. in a bilayer with the normal aligned with the Z axis (a pre-requisite for bilayer simulations). An assumed consequence is that lipids in a bilayer rotate around their respective long axis easily, but rarely invert, which would usually involve switching leaflets. The X and Y axis projections of bond vectors are expected to be averaged due to the rotation around the long axis, and the formula for the order parameter simplifies to
S(CD) = (3*z**2 - 1)/2
where z is the Z axis projection of the normalized bond vector (a unit vector). Much of the complexity of the script is bookkeeping related to the fact that DPPC has 80 C-H bond vectors to consider, and that one must consider all of the lipid molecules in the simulation system. I broke it down into 3 scripts to keep each simpler, and for efficiency reasons as well; to compute S(CD) for chain position 5 for an 80 lipid system requires applying the above equation to 320 vectors.
Nvect = (2 C-H bonds) * (2 chains) * (80 lipids) = 320
One script computes S(CD) for the CH2 groups of the alpha chain (the ethanolamine fragment), a second does the 5 vectors for the glycerol, and and the third computes it for the 4 C-H vectors at a specific chain position. The chain script requires an extra command line argument, which is the the chain position number (2-15 for DPPC); all the scripts expect an argument for the number of files. Assuming 40 .trj files (also called .dcd files) and input filenames of scd-eth.inp, scd-gly.inp, and scd-chn.inp
charmm N:40 < scd-eth.inp >& scd-eth.out charmm N:40 C:5 < scd-chn.inp >& scd-chn5.out
_________________________
Rick Venable computational chemist
|
Top
|
|
|
|
#6858 - 05/25/05 10:01 PM
Re: lipid C-D order parameters using CORREL
[Re: rmv]
|
Forum Member
Registered: 09/17/03
Posts: 8406
Loc: 39 03 48 N, 77 06 54 W
|
scd-eth.inp; extra comments * axial avg order parameter for ethanolamine C:H vectors; @N files *
bomblev -1
! READ STD FILES FOR PROJECT; RTF, PARAM, PSF, COOR stream rtfprm.str stream psfcrd.str
! COMPUTE CORREL LIMITS calc mxa 80 * 6 * 4 calc mxs 10 + 80*2 * 4 calc mxt 200 * @N
! INVOKE CORREL correl maxtim @MXT maxa @MXA maxs @MXS ! ACCUMULATE AVERAGES IN THESE TIME SERIES; DECLARE IN THIS ORDER enter sr zero enter ss dupl sr enter sx dupl sr enter sy dupl sr ! LOOP OVER THE LIPIDS; USE ONLY Z AND R OF VECTOR set k 1 label elp enter w@K vect z L @K H11A L @K C11 enter a@K vect r L @K H11A L @K C11 enter x@K vect z L @K H11B L @K C11 enter b@K vect r L @K H11B L @K C11 enter y@K vect z L @K H12A L @K C12 enter c@K vect r L @K H12A L @K C12 enter z@K vect z L @K H12B L @K C12 enter d@K vect r L @K H12B L @K C12 incr k by 1 if k le 80 goto elp
! START FROM UNIT 8; ALLOWS UP TO CA. 90 FILES set k 1 label klp calc u @K + 7 open unit @U file read name dyn@K.trj incr k by 1 if k .le. @N goto klp traj firstu 8 nunit @N
! LOOP OVER 80 LIPIDS; NORM, CALC SCD, ACCUM calc r80 1. / 80. set k 1 label clp mantim w@K ratio a@K ! DIVIDE VECTOR Z BY R (z OF UNIT VECTOR) mantim w@K squa ! z**2 mantim w@K mult 1.5 ! * 3/2 mantim w@K shift -0.5 ! - 1/2 mantim w@K mult @R80 ! * 1/Nlipid mantim sr add w@K ! ACCUMLATE, H11A mantim x@K ratio b@K mantim x@K squa mantim x@K mult 1.5 mantim x@K shift -0.5 mantim x@K mult @R80 mantim ss add x@K ! ACCUMLATE, H11B mantim y@K ratio c@K mantim y@K squa mantim y@K mult 1.5 mantim y@K shift -0.5 mantim y@K mult @R80 mantim sx add y@K ! ACCUMLATE, H12A mantim z@K ratio d@K mantim z@K squa mantim z@K mult 1.5 mantim z@K shift -0.5 mantim z@K mult @R80 mantim sy add z@K ! ACCUMLATE, H12B incr k by 1 if k le 80 goto clp ! LOOP OVER LIPIDS ENDS; COMBINE THE FOUR TIME SERIES edit sr veccod 4 delta 0.001 skip 1000 offset 1. ! STASH DATA IN A SUBDIR NAMED scd; 5 COLS, TIME IN 1st COL open write unit 1 card name scd/eth.dat write sr unit 1 dumb time * dumb * end stop
|
Top
|
|
|
|
#6861 - 08/12/05 01:40 AM
Re: lipid C-D order parameters using CORREL
[Re: rmv]
|
Forum Member
Registered: 07/19/04
Posts: 165
|
Dear Dr Venable,
Please correct my understanding.
In case of plotting between scd and position of the carbon atoms, I need to make average; say sr and ss for C11 and so on.
Is it correct?
Thank you.
|
Top
|
|
|
|
#6862 - 08/12/05 11:07 AM
Re: lipid C-D order parameters using CORREL
[Re: beginner]
|
Forum Member
Registered: 09/17/03
Posts: 8406
Loc: 39 03 48 N, 77 06 54 W
|
In general, yes, the Scd for each of the 2 C-H vectors on a given C atom are averaged, as they cannot be distinguished experimentally in most cases. The exception is C2 of each chain (C22 and C32), especially the beta chain (C22); the C-H vectors are inequivalent and can be distinguished by the NMR quadrapole splitting experiment. Computing the Scd for the 2 vectors separately is a test of convergence; if <Scd> (average over time) for 2 C-H vectors at a given chain position don't agree that well, the simulation may need to be run longer.
_________________________
Rick Venable computational chemist
|
Top
|
|
|
|
#6863 - 07/31/06 11:42 AM
Re: lipid C-D order parameters using CORREL
[Re: rmv]
|
Forum Member
Registered: 07/28/04
Posts: 113
|
Let me ask you a question on S(CD) parameter. In order to compare to experimental results, the simulation setup should be exactly the same as in experiments, e.g, the same temperature, pressure, surface tension ?
Thanks.
|
Top
|
|
|
|
#6864 - 07/31/06 02:10 PM
Re: lipid C-D order parameters using CORREL
[Re: cola]
|
Forum Member
Registered: 09/17/03
Posts: 8406
Loc: 39 03 48 N, 77 06 54 W
|
Certainly, the temperature, hydration level, and the surface area per lipid need to be fairly close, as the order parameters are known to depend on these properties.
_________________________
Rick Venable computational chemist
|
Top
|
|
|
|
#6865 - 07/31/06 03:41 PM
Re: lipid C-D order parameters using CORREL
[Re: rmv]
|
Forum Member
Registered: 07/28/04
Posts: 113
|
Which one is the most important factor among the above three ? Is it necessary to check the convergence of the order parameter ?
Thanks.
|
Top
|
|
|
|
#6866 - 08/01/06 02:38 AM
Re: lipid C-D order parameters using CORREL
[Re: cola]
|
Forum Member
Registered: 01/03/05
Posts: 86
|
Dear Prof. Rick,
I have to analyse the order parameter for more than 100 files, the script mentioned here can take maximum of 90 files. Is there any way of reading more than 100 trj files. I have tried the close unit @u also but that too doesnt seem to be working.
waiting for your reply, regards, MzZt.
|
Top
|
|
|
|
|
|