The following script uses a vector from C3 to C13 of each lipid chain to estimate the overall chain tilt.

* DPPC 80 lipid gel; chain tilt time series
*
bomblev -2

! READ RTF, PARAM, PSF, AND COOR FILES
stream rtfprm.str
stream psfcrd.str

! COMPUTE TIME SERIES LIMITS; 1 NS SAVED AT 1 PS INCR IN EACH FILE
calc mxt = 1000 * @N
! FOUR VECTORS PER LIPID + TEN MORE
calc mxs = 10 + ( 4 * 80 )

correl maxt @MXT maxs @MXS maxa ?NATOM noupdate
! DECLARE ACCUMULATION SERIES (FOR AVGS) FIRST
enter aav zero
enter bav zero
! LOOP OVER LIPIDS; DEFINE VECTOR LENGTH (R) AND Z COMPONENT
set k 1
label elp
enter r@K vect r L @K C23 L @K C213 ! BETA CHAIN
enter s@K vect r L @K C33 L @K C313 ! GAMMA CHAIN
enter z@K vect z L @K C23 L @K C213
enter w@K vect z L @K C33 L @K C313
incr k by 1
if k le 80 goto elp

! OPEN FILES (COPIED TO /scratch FOR EFFICIENCY)
set j 1
label olp
calc u = 7 + @J
open unit @U file name /scratch/dyn@J.trj
incr j by 1
if j le @N goto olp
traj first 8 nunit @N

set k 1
label tlp
! NORMALIZE Z COMPONENT (DIVIDE BY R)
mantim z@K ratio r@K
mantim w@K ratio s@K
! IGNORE SIGN; NEEDED TO AVG OVER BOTH LEAFLETS
mantim z@K abs
mantim w@K abs
! SUM OVER CHAINS, EACH SEPARATELY
mantim aav add z@K
mantim bav add w@K
incr k by 1
if k le 80 goto tlp

! DIVIDE BY N; RESULT IS AVG Z PROJECTION, WHICH IS COS(THETA)
mantim aav divi 80.
mantim bav divi 80.
! CONVERT COS TO AN ANGLE
mantim aav acos
mantim bav acos

! OUTPUT FILE WITH (TIME, BETA-TILT, GAMMA-TILT) FOR EACH TIME POINT
edit aav veccod 2
open unit 2 write card name tiltavg.dat
write aav unit 2 dumb time
* dummy
*

end
stop

Dear Dr Venable,

I would like to calculate the tilt angles in the upper/lower leaflet separately. Could you please guide me how to do that?

Thank you very much.

Most of the script would be the same, only how the individual chain vectors are averaged changes. The loop with the label 'tlp' becomes

set k 1
label tlp
mantim z@K ratio r@K
mantim w@K ratio s@K
mantim z@K abs
mantim w@K abs
if k gt 40 goto skp1
mantim aav add z@K
mantim aav add w@K
goto skp2
label skp1
mantim bav add z@K
mantim bav add w@K
label skp2
incr k by 1
if k le 80 goto tlp

In this case, it is assumed that lipids 1 thru 40 are in one leaflet, and lipids 41 thru 80 are in the other.

A final note on this script for chain tilt-- the value is likely only meaningful for systems with considerable collective tilt, usually with mostly trans alkane chains packed together, as in a gel phase. For the more disordered fluid phase, this value may not be very useful at all.

Dear Dr Venable,

Thank you very much for your advice.

Unfortunately, I got this following error message
fmt: end of file
apparent state: unit 5 (unnamed)
last format: (A)
lately reading sequential formatted external IO
Abort

I doubt it came from reading trajectories. However, it went through when I computed both leaflet at the same time. How could I solve this problem?

Thank you very much.

It's hard to say-- although UNIT 5 is standard input, which indicates a problem reading the input itself, and not some other file. Note that because of the loops, there may be problems using an MPI-based version for this; two solutions are: [1] use a non-parallel version of CHARMM, or [2] use the STREAM command to read the input.

Dear Charmmers,

I am interested in calculating the angle between the helical axis of peptide and the bilayer normal.
If I replace all the DMPCs with the peptide in which I am interested, I get the warning mesg. Am I going in the right direction of using this above script for calculating angle between helix axis and membrane normal?

regards,
MzZt.

Probably not-- the above is for averaging many vectors, and a helical axis is only a single vector. A time series of the principal axis of the helix backbone atoms from an INERTIA time series would be more what you need; the angle can be computed from the unit vector Z projection.

dear rick,

as I would like to focus on the distribution of that angle,
let say, I would like to get the histogram from 0-180 degree on the x-axis and probability on the y-axis.

and I did it via following script (based on your script), is it correct?

enter aav zero
edit aav total 180 delta 1. offset 0 skip 1
enter bav dupl aav

! LOOP OVER LIPIDS; DEFINE VECTOR LENGTH (R) AND Z COMPONENT
set k 1
label elp
enter r@K vect r L @K C23 L @K C213 ! BETA CHAIN
enter s@K vect r L @K C33 L @K C313 ! GAMMA CHAIN
enter z@K vect z L @K C23 L @K C213
enter w@K vect z L @K C33 L @K C313
incr k by 1
if k le 224 goto elp
! OPEN FILES (COPIED TO /scratch FOR EFFICIENCY)
open unit 8 file name lipid.dcd
traj first 8 nunit 1 skip 2500
set k 1
label tlp
! NORMALIZE Z COMPONENT (DIVIDE BY R)
mantim z@K ratio r@K
mantim w@K ratio s@K
! IGNORE SIGN; NEEDED TO AVG OVER BOTH LEAFLETS
mantim z@K acos
mantim w@K acos
mantim z@K hist 0 180 180
mantim w@K hist 0 180 180
mantime aav add z@K
mantime bav add w@K
incr k by 1
if k le 224 goto tlp
mantim aav divi 224.
mantim bav divi 224.
!mantim aav acos
!mantim bav acos
! OUTPUT FILE WITH (TIME, BETA-TILT, GAMMA-TILT) FOR EACH TIME POINT
edit aav total 180 delta 1 offset 0 skip 1 veccod 2
open unit 2 write card name cahin_distribu.dat
write aav unit 2 dumb time
* dummy
*

I don't see any obvious errors, but I usually allow the output to guide me.

Hi Dr. Venerable,

why not just use coor helix?

COOR HELIX does not process trajectory files; in general, an analysis using CORREL will be much more efficient than a loop over frames.

There often can be a couple different ways to do a particular analysis, and sometimes it's useful to try more than one and compare the results.

There is a helix timeseries in correl, but this is (just as COOR HELIx) intended for helical structures, and may not be appropriate for a fatty acid hydrocarbon tail.

In either case, it should be noted that the tilt calculation is only meaningful when the chains are fairly ordered, such as more gel-like states.

Dear rick,

Hello, I have a question:
Actually you know, DPPC lipid has 16 carbons in each of chain.
But, script in this archive can calculate tilt angle from C3 to C13.
Are they have specific reason?
Why you don't calculate from C1 to C16?

The free ends of the chains are highly disordered, and the upper ends near the glycerol are influenced by the ester binding geometry and headgroup packing.

Dear Rick,

Thanks, I understood that in case of C3.
However, I think that setting the end to end vector (C3-C16) is more reasonable to check the tilt angle.
Is it does not matter if all values are calculated from C3 to C13?

What matters most is that you are consistent in your methods, and clearly state what you've done in reports and publications.

Because of the flexibility of the free chain ends, I do not agree that using the terminal methyl is more reasonable; if you are concerned, do the calculation both ways and see how much difference it makes.