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
