* MD on 5'-cgcgaattcgcg-3' ecor1 dna duplex in a box of water and counterions
* Particle Mesh Ewald to treat long range electrostatics
* NPT ensemble via the Langevin Piston
*

faster on

! identifiers
set id1 ecor_sod
set id2 test

! cell parameters
set a 60.8529
set b 36.9108
set c 36.9108

! 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_all27_na.rtf
set 2 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 @id1_@id2.psf 
read psf cards unit 20

open unit 20 read form name @id1_@id2_b.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 32 fftz 32 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

shake bonh para

!minimization prior to NPT dynamics
mini ABNR nstep 100 

! identifiers
set g 0     ! current dynamics identifier

! open files for restart, trajectory and energies
open unit 31 write form name @id1_@id2_@g.rst
open unit 32 write unfo name @id1_@id2_@g.trj
open unit 33 write form name @id1_@id2_@g.ene

shake bonh para

set nstep 100

! molecular dynamics, hot start
dynamics cpt leap start time 0.002 nstep @nstep iseed 314159 -
        firstt 300.0 finalt 300.0 tstruc 300.0 -
        inbfrq -1 imgfrq -1 iprfrq 50 nprint 50 -
        iasors 1 iasvel 1 ieqfrq 0 ntrfrq 50 -
        iunrea -1 iunwri 31 iuncrd 32 kunit 33 nsavc 500 isvfrq 5000 -
        pconstant pmass 600.0 pref 1.0 isot pgamma 0.0 tbath 300.0 -
        hoover reft 300.0 tmass 1000.0 echeck 100000.

set f 0
set g 50

label dyn_loop

! open files for restart, trajectory and energies
open unit 31 write form name @id1_@id2_@g.rst
open unit 32 write unfo name @id1_@id2_@g.trj
open unit 33 write form name @id1_@id2_@g.ene
open unit 30 read form name @id1_@id2_@f.rst

shake bonh para

! molecular dynamics, hot start
dynamics cpt leap restart time 0.002 nstep @nstep iseed 314159 -
        firstt 300.0 finalt 300.0 tstruc 300.0 -
        inbfrq 10 imgfrq 10 iprfrq 50 nprint 50 -
        iasors 1 iasvel 1 ieqfrq 0 ntrfrq 50 -
        iunrea 30 iunwri 31 iuncrd 32 kunit 33 nsavc 500 isvfrq 5000 -
        pconstant pmass 600.0 pref 1.0 isot pgamma 0.0 tbath 300.0 -
        hoover reft 300.0 tmass 1000.0 echeck 100000.

incr g by 50
incr f by 50

if g le 100 goto dyn_loop

stop