* 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 per strand
set restot 12

! 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  1.0 ! cutim
set 4  1.0 ! cutnb
set 5  0.5 ! ctonnb
set 6  0.7 ! 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

! coordinates required for correl
open read unit 12 card name @id1_@id2_b.crd
read coor card unit 12

update -
inbfrq -1 imgfrq -1 ihbfrq 0 -
@7 @8 @9 vswitch cutimg @3 cutnb @4 ctofnb @6 ctonnb @5

! set input residue numbers
set 1 1          ! residue number 1
set 2 dna1       ! segid 1
set 3 @restot    ! residue number 2
set 4 dna2       ! segid 2

! loop for the individual residues
label loop

! initial atoms for G-C
set 5 o6
set 6 n4

!reset atoms if necessary

coor stat sele segid dna1 .and. resi @1 end
set resname ?SELRESN
if resname eq CYT set 5 n4
if resname eq CYT set 6 o6
if resname eq THY set 5 n3
if resname eq THY set 6 n1
if resname eq URA set 5 n3
if resname eq URA set 6 n1
if resname eq ADE set 5 n1
if resname eq ADE set 6 n3

correl maxser 6 maxt 20000

enter dis1 vect r @2 @1 @5 @4 @3 @6  ! cyt n4 to gua o6
                                     ! ade n1 to thy n3

!offset time series if reading from a trajectory
!file that does not start from 0 time
!edit dis1 offset 1001. skip 500 delta 0.002

open unit 21 write form name dist/@id1_@id2_@1_@5_@6.dist

! initial atoms for G-C
set 5 n1
set 6 n3

if resname eq CYT set 5 n3
if resname eq CYT set 6 n1
if resname eq THY set 5 o4
if resname eq THY set 6 n6
if resname eq URA set 5 o4
if resname eq URA set 6 n6
if resname eq ADE set 5 n6
if resname eq ADE set 6 o4

enter dis2 vect r @2 @1 @5 @4 @3 @6  ! cyt n3 to gua n1
                                     ! ade n6 to thy o4

!edit dis2 offset 1001. skip 500 delta 0.002

open unit 22 write form name dist/@id1_@id2_@1_@5_@6.dist

if resname eq THY goto skip_third ! for A-T
if resname eq ADE goto skip_third ! for A-T
if resname eq URA goto skip_third ! for A-T

! initial atoms for G-C
set 5 n2
set 6 o2

if resname eq CYT set 5 o2
if resname eq CYT set 6 n2

enter dis3 vect r @2 @1 @5 @4 @3 @6  ! cyt o2 to gua n2

!edit dis3 offset 1001. skip 500 delta 0.002

open unit 23 write form name dist/@id1_@id2_@1_@5_@6.dist

label skip_third

! 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

! 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

prnlev -5
if b le @trajf goto traj_loop1 
prnlev 5

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

write dis1 unit 21 dumb time
write dis2 unit 22 dumb time

if resname eq THY goto skip_third2 ! for A-T
if resname eq ADE goto skip_third2 ! for A-T
if resname eq URA goto skip_third2 ! for A-T

write dis3 unit 23 dumb time
label skip_third2

end

incr 1 by 1
decr 3 by 1

if 1 le @restot goto loop

stop