* To analyse the puckering in ecor1 residu by residue,
* in both strands
* For convenience, we read a single (merged) trajectory,
* obtained during analysis of RMS deviations. 
* Puckering is analyzed every ps
*                                    

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 one 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  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

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

correl maxseries 48 maxtime 50000 maxatoms 7000

! We loop through all the residues of the first strand,
! and then again through all the residues of the second strand

set segid dna1 
set resid 1
set u 30  ! unit number where to write the puckering data

!============== Analysing the puckering in strand 1 ===============

set a 50  ! unit number
set b 0 ! 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

if b le @trajf goto traj_loop1 

!loop to enter timeseries
label strand1

enter p@resid puck resi @resid segi @segid

incr resid by 1

if resid le @restot goto strand1

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

!loop to write timeseries
set resid 1
label strand1b 

open write unit @u form name @id1_@id2_@segid_@resid_@range.puc
write p@resid unit @u dumb time 

increase resid by 1
increase u by 1

if resid le @restot goto strand1b

end

!============== Analysing the puckering in strand 2 ==============

set segid dna2
set resid 1
set u 30  ! unit number where to write the puckering data

correl maxseries 48 maxtime 50000 maxatoms 7000

set a 50  ! unit number
set b 0 ! trajectory start
set c 0   ! counter for total number of traj files

! loop over all the trajectory files
label traj_loop2

! 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 ! 695: last part of trajectory

label strand2
enter p@resid puck resi @resid segi @segid
incr resid by 1
if resid le @restot goto strand2

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

set resid 1
label strand2b

open write unit @u form name @id1_@id2_@segid_@resid_@range.puc
write p@resid unit @u dumb time

increase resid by 1
increase u by 1

if resid le @restot goto strand2b

end ! end of correl

stop