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

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

scd-gly.inp; see comments in scd-eth.inp

* the axial avg order parameter for glycerol C:H vectors; @N files
*

bomblev -2
stream rtfprm.str
stream psfcrd.str

calc mxa 80 * 6 * 5
calc mxs 10 + 80*2 * 5
calc mxt 200 * @N

correl maxtim @MXT maxa @MXA maxs @MXS
enter sa zero
enter sb dupl sa
enter ss dupl sa
enter sx dupl sa
enter sy dupl sa
set k 1
label elp
enter v@K vect z L @K HS  L @K C2
enter r@K vect r L @K HS  L @K C2
enter w@K vect z L @K HA  L @K C1
enter a@K vect r L @K HA  L @K C1
enter x@K vect z L @K HB  L @K C1
enter b@K vect r L @K HB  L @K C1
enter y@K vect z L @K HX  L @K C3
enter c@K vect r L @K HX  L @K C3
enter z@K vect z L @K HY  L @K C3
enter d@K vect r L @K HY  L @K C3
incr k by 1
if k le 80 goto elp

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
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim w@K mult @R80
mantim sa add w@K
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 sb add x@K
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
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
mantim v@K ratio r@K
mantim v@K squa
mantim v@K mult 1.5
mantim v@K shift -0.5
mantim v@K mult @R80
mantim ss add v@K
incr k by 1
if k le 80 goto clp
edit sa veccod 5 delta 0.001 skip 1000 offset 1.
open write unit 1 card name scd/gly.dat
write sa unit 1 dumb time
* dumb
*
end

scd-chn.inp; bonus comments included; see also scd-eth.inp

* the axial avg order parameter for an input C no.; @N files
* additional required arg C: chain carbon number @C
*

bomblev -2
stream rtfprm.str
stream psfcrd.str

calc mxa 80 * 6 * 4
calc mxs 10 + 80*2 * 4
calc mxt 200 * @N

correl maxtim @MXT maxa @MXA maxs @MXS
enter sr zero 
enter ss dupl sr
enter sx dupl sr
enter sy dupl sr
set k 1
! TRAILING LETTERS IN THE 'H' ATOM NAMES REQUIRE @{C}, NOT @C
label elp
enter w@K vect z L @K H@{C}R  L @K C2@C
enter a@K vect r L @K H@{C}R  L @K C2@C
enter x@K vect z L @K H@{C}S  L @K C2@C
enter b@K vect r L @K H@{C}S  L @K C2@C
enter y@K vect z L @K H@{C}X  L @K C3@C
enter c@K vect r L @K H@{C}X  L @K C3@C
enter z@K vect z L @K H@{C}Y  L @K C3@C
enter d@K vect r L @K H@{C}Y  L @K C3@C
incr k by 1
if k le 80 goto elp

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 LIPIDS; NORM, CALC SCD, ACCUM
calc r80 1. / 80.
set k 1
label clp
mantim w@K ratio a@K
mantim w@K squa
mantim w@K mult 1.5
mantim w@K shift -0.5
mantim w@K mult @R80
mantim sr add w@K
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
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
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
incr k by 1
if k le 80 goto clp
edit sr veccod 4 delta 0.001 skip 1000 offset 1.
! MIXED CASE IS IGNORED; @C IS THE CHAIN CARBON NO.
open write unit 1 card name scd/v@C.dat
write sr unit 1 dumb time
* dumb
*
end
stop

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.

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.

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.

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.

Which one is the most important factor among the above three ?
Is it necessary to check the convergence of the order parameter ?

Thanks.

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.