* To calculate RMS deviations between simulated DNA and ideal 
* A and B forms of ecor1
*                                    

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

set 1 top_all27_na.rtf
set 2 par_all27_na.prm

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

!merging the trajectory files into a single trj file, to allow
!selection of subsets of atoms

! define the subset of atoms to keep in the merged trajectory
define comp_set sele ((segid DNA1 .and. resid 2:11) .or. -
       (segid DNA2 .and. resid 2:11)) .and. .not. (type H*) -
       end
 
! We open the trajectory files through a loop:

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

label traj_loop

	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_loop 

! open file to which trajectory will be written
open unit 30 write unform name @id2_entire.trj

! merge the pieces of a trajectory into a single file:
merge coor sele comp_set end -
      first 50 nunit @c skip @skip nfile @nfile begin @begin -
      stop @stop output 30

!goto skip_b

!=================== Comparaison vs B DNA ======================

!read ideal B conformation coordinates
open read unit 12 card name ecor1_ideal_b.crd
read coor card unit 12

! Keep only selected atoms
delete atom sele .not. comp_set end

coor copy comp

correl maxseries 6 maxtime 20000 maxatoms 800

enter rmsb rms orient 

! read the merged trajectory:
open read unit 30 unform name @id2_entire.trj
trajectory firstu 30 nunit 1 begin @begin stop @stop skip @skip 

!edit rmsb offset 501. skip 500 delta 0.002

open write unit 31 form name @id1_@id2_vs_b.rms
write rmsb unit 31 dumb time

end ! end of correl

!=================== Comparaison vs A DNA ======================

delete atom sele all end

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

!label skip_b

! We read coordinates of ecor1 in ideal A conformation
open read unit 13 card name ecor1_ideal_a.crd
read coor card unit 13

! Keep only selected atoms
delete atom sele .not. comp_set end

coor copy comp

correl maxseries 6 maxtime 20000 maxatoms 800

enter rmsa rms orient

! read the merged trajectory:
open read unit 30 unform name @id2_entire.trj
trajectory firstu 30 nunit 1 begin @begin stop @stop skip @skip

!edit rmsa offset 501. skip 500 delta 0.002

open write unit 31 form name @id1_@id2_vs_a.rms
write rmsa unit 31 dumb time

end ! end of correl

stop