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

bomlev -1

faster on

! 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

calc nfile = (@stop - @begin + @skip)/@skip

!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

! We read coordinates of ecor1 to allow correl to work
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 a1 angl @q c1' @q c2' @q c3'
enter a2 angl @q c2' @q c3' @q c4'
enter a3 angl @q c3' @q c4' @q o4'
enter a4 angl @q c4' @q o4' @q c1'
enter a5 angl @q o4' @q c1' @q c2'
!enter a6 angl @q c1' @q c2' @q o3'
!enter a7 angl @q c3' @q c2' @q o3'
enter a8 angl @q c2' @q c3' @q o3'
enter a9 angl @q c4' @q c3' @q o3'
enter a10 angl @q c5' @q c4' @q c3'
enter a11 angl @q c5' @q c4' @q o4'
enter a12 angl @q o4' @q c1' @q @m
enter a13 angl @q c2' @q c1' @q @m 

! 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 angl/@id1_@id2_@q_a1.angl
write a1 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a2.angl
write a2 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a3.angl
write a3 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a4.angl
write a4 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a5.angl
write a5 unit 20 dumb time
!open write unit 20 form name angl/@id1_@id2_@q_a6.angl
!write a6 unit 20 dumb time
!open write unit 20 form name angl/@id1_@id2_@q_a7.angl
!write a7 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a8.angl
write a8 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a9.angl
write a9 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a10.angl
write a10 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a11.angl
write a11 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_a12.angl
write a12 unit 20 dumb time
open write unit 20 form name angl/@id1_@id2_@q_13.angl
write a13 unit 20 dumb time

end

incr q by 1

if q le @restot goto ic_loop

stop