* Equilibration to allow solvent to equilibriate around the RNA
* Minimization and 20 ps NVT simulation with heavy atoms of the RNA
* harmonically constrained.
*

faster on

! identifiers
set id1 rr_sol
set id2 rr_equil
set top $CHARMMDATA
set out output

! cell parameters
set a 55.8562
set b 41.9526
set c 37.8237

! 1: topology
! 2: parameter
! 3: cutim
! 4: cutnb
! 5: ctonnb
! 6: ctofnb
! 7: electrostatic cutoff regime (shift or switch)
! 8: electrostatic pairing method (atom or group)
! 9: VDW pairing method (vatom or vgroup)
! VDW cutoff regime is vwsitch
set 1 @top/top_all27_na.rtf
set 2 @top/par_all27_na.prm
set 3 12.0 ! cutim
set 4 12.0 ! cutnb
set 5  8.0 ! ctonnb
set 6 10.0 ! ctofnb
set 7 switch
set 8 atom
set 9 vatom

! 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 20 read form name @out/@id1.psf 
read psf cards unit 20

open unit 20 read form name @out/@id1.crd
read coor cards unit 20

crystal define orthorhombic @a @b @c 90. 90. 90.

crystal build noper 0 cutoff @3

!kappa should be optimized based on the value of
!kappa at which the 1st derivative of the energy
!with respect to kappa goes to 0.0
set kap 0.36
set ord 6

!fftx dimensions should correspond to 1 A or less; must
!be power of 2, 3 or 5
update inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa @kap fftx 64 ffty 64 fftz 64 order @ord -
@7 @8 @9 vfswitch cutimg @3 cutnb @4 ctofnb @6 ctonnb @5

!turn on image centering of solvent
image byres xcen 0.0 ycen 0.0 zcen 0.0 -
     sele resn tip3 .or. resn sod .or. resn cl end

! harmonic constraints on all dna heavy atoms
cons harm force 5.0 mass sele segid rna1:rna2 end

shake bonh para

! minimize for initial 200 steps with constraints
mini SD nstep 200 nprint 5

set g 0

! open files for restart, trajectory and energies
open unit 31 write form name @out/@id2.rst
open unit 32 write unfo name @out/@id2.trj

! molecular dynamics, hot start
dyna leap strt nstep 100 timestep 0.002 ntrfrq 50 -
iprfrq 50 ieqfrq 50 inbfrq -1 imgfrq -1 ihbfrq 0 -
iunrea -1 iunwri 31 iuncrd 32 kunit -1 -
nprint 50 nsavc 5000 isvfrq 500 -
firstt 300.0 finalt 300.0 twindh 2.0 twindl 2.0 -
iasors 0 iasvel 1 iscvel 0 ichecw 1 iseed 37987473

open unit 20 write form name @out/@id2.crd
write coor cards unit 20
* RNA harmonically constrained and equilibrated
* @1
* @2
*
open unit 20 write form name @out/@id2.pdb
write coor pdb unit 20
* RNA harmonically constrained and equilibrated
* @1
* @2
*

stop