* DNA, write out nucleic acid backbone dihedral time series
*                                    

bomlev -1

! identifiers
set id1 ecor_sod
set id2 test

set range 500_1000

set trajf 950  ! final traj identifier
set traji 50   ! traj increment value
set file1 500    ! initial traj file

set skip  500
set begin 250500   ! initial time frame
set begin2 250500
set stop 500000 ! final time frame

!total number of nucleotides in BOTH strands
set restot 24

! The following cutoffs are artificially low to speed up
! the processing of the CORREL command:
set 1 top_all27_na.rtf
set 2 par_all27_na.prm
set 3  6.0 ! cutim
set 4  6.0 ! cutnb
set 5  5.0 ! ctonnb
set 6  5.5 ! ctofnb
set 7 shift
set 8 atom
set 9 vatom

open unit 9 read form name @1
read rtf card unit 9

open unit 10 read form name @2 
read para card unit 10

open unit 11 read form name @id1_@id2.psf
read psf cards unit 11

open read unit 12 card name @id1_@id2_b.crd
read coor card unit 12

ic delete sele all end

set q 1

label ic_loop

! select between purine and pyrimidine glycosidac linkage
set m n9  ! purine
set n c4  ! purine

! reset atoms defining chi for pyrmidines
!
coor stat sele ires @q end
set resname ?SELRESN
if resname eq THY set m n1
if resname eq THY set n c2
if resname eq URA set m n1
if resname eq URA set n c2
if resname eq CYT set m n1
if resname eq CYT set n c2

! for first residue
correl maxser 48 maxt 1000

enter t0 dihe @q c4' @q o4' @q c1' @q c2'
enter t1 dihe @q o4' @q c1' @q c2' @q c3'
enter t2 dihe @q c1' @q c2' @q c3' @q c4'
enter t3 dihe @q c2' @q c3' @q c4' @q o4'
enter t4 dihe @q c3' @q c4' @q o4' @q c1'
! other sugar defined by Gelbin et al.
enter n1 dihe @q c5' @q c4' @q c3' @q c2'
enter n2 dihe @q c5' @q c4' @q c3' @q o3'
enter n3 dihe @q o4' @q c4' @q c3' @q o3'
enter n4 dihe @q c4' @q o4' @q c1' @q @m
enter n6 dihe @q o5' @q c5' @q c4' @q c3'
enter n7 dihe @q o4' @q c1' @q @m  @q @n

! loop to open trj files
prnlev 1
set a 50     ! unit number
set b @file1 ! trajectory start
set c 0      ! counter for total number of traj files

! loop over all the trajectory files
label traj_loop1
prnlev 3

! read the ith trajectory:
open read unit @a unform name @id1_@id2_@b.trj

        increase a by 1
        increase b by @traji
        increase c by 1

if b le @trajf goto traj_loop1 ! 695: last part of trajectory

trajectory firstu 50 nunit @c begin @begin2 stop @stop skip @skip

open write unit 20 form name dihe/@id1_@id2_@q_t0.dihe
write t0 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_t1.dihe
write t1 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_t2.dihe
write t2 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_t3.dihe
write t3 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_t4.dihe
write t4 unit 20 dumb time

open write unit 20 form name dihe/@id1_@id2_@q_n1.dihe
write n1 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_n2.dihe
write n2 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_n3.dihe
write n3 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_n4.dihe
write n4 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_n6.dihe
write n6 unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_n7.dihe
write n7 unit 20 dumb time

end

incr q by 1

if q le @restot goto ic_loop

stop