This example illustrates computing an angle theta wrt. the bilayer normal (Z axis), and getting the probability distribution, p(theta), for each leaflet.

* the angle from the Z proj for P--N; first file @FF last file @LF

*

stream rtfprm.str ! READ RTF,PARAM

stream psfcrd.str ! READ PSF,COOR

! CHANGE FOR OTHER SIZES; LIPID COUNT, AND HALFWAY POINT

set nlpd 72

set nlph 36

! CALC SIZES FOR CORREL

calc n = 1 + @LF - @FF

calc mxa = @NLPD * 6

calc mxs = 10 + @NLPD * 2

calc mxt = 100 * @N ! 100 FRAMES PER FILE

correl maxtim @MXT maxa @MXA maxs @MXS noupdate

enter sa zero

edit sa total 360 offset 0. delta 0.5 skip 1

enter sb dupl sa

! DECLARE Z, R VECTOR COMPONENTS

set k 1

label elp

enter z@K vect z L @K N L @K P

enter r@K vect r L @K N L @K P

incr k by 1

if k le @NLPD goto elp

! OPEN FILES

set m @FF

set u 8

label tlp

open read unit @U file name dyn@M.trj

incr u by 1

incr m by 1

if m le @LF goto tlp

traj firstu 8 nunit @N

! LOOP OVER LIPIDS; LEAFLET A = 1 THRU NLPH, B > NLPH

set k 1

label clp

set inv 1.0

if k gt @NLPH set inv -1.0 ! SIGN INVERTED FOR LEAFLET B

mantim z@K ratio r@K ! NORMALIZE (DIVIDE BY LENGTH)

mantim z@K mult @INV ! LEAFLET CORRECTION

mantim z@K acos ! COMPUTE ANGLE THETA

mantim z@K hist 0. 180. 360 ! COMPUTE p(THETA); LEAFLET AVG

if k le @NLPH mantim sa add z@K

if k gt @NLPH mantim sb add z@K

incr k by 1

if k le @NLPD goto clp

! NORMALIZE SO THAT p(theta) SUMS TO 1.0

mantim sa divi @NLPH

mantim sb divi @NLPH

edit sa veccod 2 delta 0.5 skip 1 offset 0. total 360

open write unit 1 card name pn-hst@LF.dat

write sa unit 1 dumb time

* dumb

*

end

stop