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

! identifiers
set id1 ecor_sod
set id2 test
set id3 500_1000

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

set nfile @stop
divi nfile by @skip

bomlev -1

faster on

! 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

ic delete sele all end

! generate nucleic acid backbone and glycosidac linkage ic's

! set variables based on particular sequence
set restot 10
set first1 1
calc first2 = @restot + 1
calc last1 = @restot
calc last2 = @restot * 2
calc total = @restot * 2

set q 1
set r 2
set s 3

! second counter for individual strands
set t 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

if q eq @first1 goto ic_first
if q eq @first2 goto ic_first
goto skip_first
label ic_first
! for first residue

correl maxser 48 maxt 2000

enter chi dihe @q o4' @q c1' @q @m  @q @n 
enter gam dihe @q o5' @q c5' @q c4' @q c3' 
enter del dihe @q c5' @q c4' @q c3' @q o3' 
enter eps dihe @q c4' @q c3' @q o3' @r p 
enter zet dihe @q c3' @q o3' @r p   @r o5' 

! 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_chi.dihe
write chi unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_gam.dihe
write gam unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_del.dihe
write del unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_eps.dihe
write eps unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_zet.dihe
write zet unit 20 dumb time

end

goto skip_all

label skip_first

if q eq @last1 goto ic_last
if q eq @last2 goto ic_last
goto skip_last
label ic_last
set s @q
decr s by 1
ic edit
! for last residue in strand
dihe @q o4' @q c1' @q @m  @q @n  0.0 ! chi, glycosidac linkage
dihe @s o3' @q p   @q o5' @q c5' 0.0 ! alpha
dihe @q p   @q o5' @q c5' @q c4' 0.0 ! beta
dihe @q o5' @q c5' @q c4' @q c3' 0.0 ! gamma
dihe @q c5' @q c4' @q c3' @q o3' 0.0 ! delta
end

correl maxser 48 maxt 2000

enter chi dihe @q o4' @q c1' @q @m  @q @n 
enter alp dihe @s o3' @q p   @q o5' @q c5'
enter bet dihe @q p   @q o5' @q c5' @q c4'
enter gam dihe @q o5' @q c5' @q c4' @q c3' 
enter del dihe @q c5' @q c4' @q c3' @q o3' 

! 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_loop3
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_loop3

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

open write unit 20 form name dihe/@id1_@id2_@q_chi.dihe
write chi unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_alp.dihe
write alp unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_bet.dihe
write bet unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_gam.dihe
write gam unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_del.dihe
write del unit 20 dumb time

end

goto skip_all
label skip_last

set s @q
decr s by 1
! for central residues

correl maxser 48 maxt 2000

enter chi dihe @q o4' @q c1' @q @m  @q @n 
enter alp dihe @s o3' @q p   @q o5' @q c5'
enter bet dihe @q p   @q o5' @q c5' @q c4'
enter gam dihe @q o5' @q c5' @q c4' @q c3' 
enter del dihe @q c5' @q c4' @q c3' @q o3' 
enter eps dihe @q c4' @q c3' @q o3' @r p 
enter zet dihe @q c3' @q o3' @r p   @r o5' 

! 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_loop2
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_loop2

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

open write unit 20 form name dihe/@id1_@id2_@q_chi.dihe
write chi unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_alp.dihe
write alp unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_bet.dihe
write bet unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_gam.dihe
write gam unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_del.dihe
write del unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_eps.dihe
write eps unit 20 dumb time
open write unit 20 form name dihe/@id1_@id2_@q_zet.dihe
write zet unit 20 dumb time

end

label skip_all

incr q by 1
incr r by 1

if q le @total goto ic_loop

stop