* Generate NAD+ and propionaldehyde for QM/MM calculations
* minimize and perform dynamics using 2 umbrella constraints
*

faster on

set id1 nad_pral_1

! 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 a top_all27_prot_na_qm.rtf
set b par_all27_prot_na_qm.prm
set 3 999.0 ! cutim
set 4 999.0 ! cutnb
set 5 800.0 ! ctonnb
set 6 950.0 ! ctofnb
!set 7 fshift
set 7 switch
set 8 atom
set 9 vatom

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

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

!generate propionaldehyde
read sequence card
* propionaldehyde
*
1
pral

generate alde first none last none setup warn

!generate NAD+
read sequence card
* nad
*
1
nad

generate nada first none last none setup warn

open unit 19 read form name nad_pral_qm_min.crd
read coor card unit 19 

!define the QM atoms, including propionaldehyde, NAD
!nicotinamide atoms and NA nc1' carbon in the sugar
!
define qms sele resn pral .or. (resn nad .and. (type NN1 .or. -
type NC6 .or. type NH6 .or. type NC5 .or. type NH5 .or. -
type NC4 .or. type NH4 .or. type NC3 .or. type NC2 .or. -
type NH2 .or. type NC7 .or. type NO7 .or. type NN7 .or. -
type NH71 .or. type NH72 .or. type nc1')) show end

quantum sele qms end glnk sele resn nad .and. type nc1' end -
am1 charge 1 scfc 0.00001

!constrain MM region of NAD+, excluding hydrogens
cons harm force 5.0 sele .not. (qms .or. hydrogen) end

update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

!constraint distances
set cdist1 3.0
set cdist2 1.0

!loop over constraint region using langevin dynamics

!identifiers that avoid decimal points
calc cdistid1 = @cdist1 * 10
calc cdistid2 = @cdist2 * 10

!set up distance constraints
MMFP
!NAD nc4 to propionaldehyde H distance
Geo sphere distance harm symmetric force 100.0 droff @cdist1 -
sele atom nada 1 nc4 end sele atom alde 1 h11 end

!propionaldehyde C to H distance
Geo sphere distance harm symmetric force 100.0 droff @cdist2 -
sele atom nada 1 h11 end sele atom alde 1 c1 end

END

!initial minimization to satisfy constraints
mini abnr nstep 100 nprint 50

!langevin dynamics, fbeta to mimic aqueous solution
scalar fbeta set 60.0 sele all end

open unit 30 write form name dyn_@id1_@cdistid1_@cdistid2_eq.rst

!2 ps of equilibration
DYNA LEAP LANGEVIN STRT  NSTEP 100 TIMESTEP 0.001 -
    IPRFRQ 100 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0  -
    IUNREA -1 IUNWRI 30 IUNCRD -1 IUNVEL -1 KUNIT -1 -
    NPRINT 100 NSAVC 0 NSAVV 0 -
    ILBFRQ 1000 RBUFFER 0.0 TBATH 300.0 -
    FIRSTT 300.0 FINALT 300.0  -
    IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0

open unit 30 write form name dyn_@id1_@cdistid1_@cdistid2_1.rst
open unit 31 write unfo name dyn_@id1_@cdistid1_@cdistid2_1.trj
open unit 29 read form name dyn_@id1_@cdistid1_@cdistid2_eq.rst

!10 ps of sampling
DYNA LEAP LANGEVIN RESTRT NSTEP 100 TIMESTEP 0.001 -
    IPRFRQ 100 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0  -
    IUNREA 29 IUNWRI 30 IUNCRD 31 IUNVEL -1 KUNIT 70 -
    NPRINT 100 NSAVC 10 NSAVV 0 -
    ILBFRQ 1000 RBUFFER 0.0 TBATH 300.0 -
    FIRSTT 300.0 FINALT 300.0  -
    IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0

!write final coordinate frame
open unit 20 write form name dyn_@id1_@cdistid1_@cdistid2.crd
write coor card unit 20
* final time frame from sampling window
* QM section: propionaldehyde, nicotinamide and nc1' atom
* 500 abnr w/ cons harm 5 not qms .or. hydrogen
* nc4-h11: @cdist1
* c1-h11: @cdist2
*

stop