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
