* calculate rms differences with respect to the NMR structure
*

faster on

! identifiers
set id1 rr_sol
set id2 rr
set top $CHARMMDATA 
set psfcrd ../output !directory where psf and crd files exist
set trjdir ../output ! directory for the trajectories

! 1: topology
! 2: parameter
set 1 @top/top_all27_na.rtf
set 2 @top/par_all27_na.prm

! read topology and parameters
open unit 9 read form name @1
read rtf card unit 9

open unit 9 read form name @2
read para card unit 9

open unit 10 read form name @psfcrd/@id1.psf
read psf cards unit 10

open unit 10 read form name @psfcrd/@id1.crd
read coor cards unit 10
coor copy comp

open unit 11 write form name @id2.rms

set t 1.0  ! time

set begin 100
set stop 5000
set skip 100
calc time = @begin * 0.002
calc stime = @skip * 0.002
calc time_max = @nstep * 0.002

set trajs 1
set traji 1
set trajf 3
set unit 13
set trajn 0

label trj_loop

open unit @unit unform read name @trjdir/@id2_prod_@trajs.trj

calc trajs = @trajs + @traji
calc trajn = @trajn + @traji
calc unit = @unit + 1

if @trajs le @trajf goto trj_loop

traj iread 13 nread @trajn nskip @skip begin @begin stop @stop

label rms_loop

traj read

coor orient rms sele (segid RNA1 .and. resid 2:9 .and. .not. hydrogen) -
                .or. (segid RNA2 .and. resid 2:9 .and. .not. hydrogen) show end

write title unit 11
* @time ?rms
*

calc time = @time + @stime
if time le @time_max goto rms_loop

stop