An added script that illustrates detection of C=C double bonds via the CEL1 chemical atom type; 72 DOPC lipids, segid D.
* the axial avg order parameter for an input C no. @C
* modified for C=C case; first file @FF last file @LF
* subdir @S
*
stream rtfprm.str
stream psfcrd.str
! number of lipids
set nl = 72
! number of files
calc n = 1 + ( @LF - @FF )
calc mxa @NL * 6 * 4
calc mxs 10 + @NL*2 * 4
calc mxt 500 * @N
! unsaturation check; beta chain
set ub = 0
define unsat sele atom D 1 C2@C .and. chem CEL1 end
if ub lt ?NSEL set ub ?NSEL
! unsaturation check; gamma chain
set ug = 0
define unsat sele atom D 1 C3@C .and. chem CEL1 end
if ug lt ?NSEL set ug ?NSEL
correl maxtim @MXT maxa @MXA maxs @MXS
enter sr zero
if ub eq 0 enter ss dupl sr
enter sx dupl sr
if ug eq 0 enter sy dupl sr
set k 1
label elp
enter w@K vect z D @K H@{C}R D @K C2@C
enter a@K vect r D @K H@{C}R D @K C2@C
if ub eq 0 enter x@K vect z D @K H@{C}S D @K C2@C
if ub eq 0 enter b@K vect r D @K H@{C}S D @K C2@C
enter y@K vect z D @K H@{C}X D @K C3@C
enter c@K vect r D @K H@{C}X D @K C3@C
if ug eq 0 enter z@K vect z D @K H@{C}Y D @K C3@C
if ug eq 0 enter d@K vect r D @K H@{C}Y D @K C3@C
incr k by 1
if k le @NL goto elp
! OPEN FILES, FILL SERIES VIA TRAJ
set k = @FF
set u = 8
label opnlp
open read unit @U file name @S/dyn@K.trj
incr u by 1
incr k by 1
if k le @LF goto opnlp
traj firstu 8 nunit @N
! LOOP OVER LIPIDS; NORM, CALC SCD, ACCUM
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 sr add w@K
if ub eq 0 then
mantim x@K ratio b@K
mantim x@K squa
mantim x@K mult 1.5
mantim x@K shift -0.5
mantim ss add x@K
endif
mantim y@K ratio c@K
mantim y@K squa
mantim y@K mult 1.5
mantim y@K shift -0.5
mantim sx add y@K
if ug eq 0 then
mantim z@K ratio d@K
mantim z@K squa
mantim z@K mult 1.5
mantim z@K shift -0.5
mantim sy add z@K
endif
incr k by 1
if k le @NL goto clp
! use final mantime to set AVER, FLUC for output
open write unit 1 card name scd/chn@{C}scd.txt
echu 1
mantim sr divi @NL
echo sr@C ?AVER ?FLUC
if ub eq 0 then
mantim ss divi @NL
echo ss@C ?AVER ?FLUC
endif
mantim sx divi @NL
echo sx@C ?AVER ?FLUC
if ug eq 0 then
mantim sy divi @NL
echo sy@C ?AVER ?FLUC
endif
close unit 1
! *** uncomment the next lines to write out the time series
! edit sr veccod 4 delta 0.001 skip 1000 offset 1.
! open write unit 1 card name scd/vchn@C.dat
! write sa unit 1 dumb time
! * dumb
! *
end
stop
