* Generate NAD and perform energy surface about the
* glycosidac linkage
* QM
*

faster on

! 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
!use infinite cutoffs
set 3 999.0 ! cutim
set 4 999.0 ! cutnb
set 5 800.0 ! ctonnb
set 6 900.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 NAD+
read sequence card
* nad
*
1
nadq

generate nada first none last none setup warn

!read coordinates with dihedral at starting values
!for the surface
open unit 19 read form name nad_qm_min.crd
read coor card unit 19

!define QM atoms.  these will also be used for
!deining the constrained region in the MM surface
define qms sele resn nadq .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

!get atom numbers for quick commands
coor stat sele type nc2 end
set a1 ?selatom
coor stat sele type nn1 end
set a2 ?selatom
coor stat sele type nc1' end
set a3 ?selatom
coor stat sele type no4' end
set a4 ?selatom

!obtain initial value of glycosidac linkage
quick @a1 @a2 @a3 @a4

!define QM region and set NC1' to be the link atom
!using the GHO frontier orbital method
quantum sele qms end glnk sele resn nadq .and. type nc1' end -
am1 charge 1 scfc 0.00001

!harmonic constraints except on nicotinamide (QM atoms)
!these are required to avoid complications in the 
!energy surface associated with additional degrees
!of freedom in the molecule
cons harm force 5.0 mass sele .not. qms end

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

!perform initial minimization to get energies to 
!offset energies in the surface
mini abnr nstep 500 nprint 100

set a ?ener
set b ?elec
set c ?vdw
set d ?dihe
set e ?angl
set f ?bond
set g ?urey
set h ?impr
set i ?qmel
set j ?qmvd

!put bond, angle, urey-bradley and improper terms together
calc itot = ?angl + ?bond + ?urey + ?impr

!perform mulliken population analysis
mulliken

!read coordinates with dihedral at starting values
!for the surface
open unit 19 read form name nad_qm_min_n180.crd
read coor card unit 19

open unit 22 write form name nad_qm_surf_1.map

!reset harmonic constraints 
cons harm force 5.0 mass sele .not. qms end

set i -180.0

label loop

!set dihedral constraint
cons dihe bynu @a1 @a2 @a3 @a4 min @i force 10000.

!minimize, remove constraint and get energy
mini abnr nstep 500 nprint 100 tolgrd 0.01
set force ?grms
cons cldh
energy

set 5 ?ener
decr 5 by @a
set 6 ?elec
decr 6 by @b
set 7 ?vdw
decr 7 by @c
set 8 ?dihe
decr 8 by @d
calc itot2 = ( ?angl + ?bond + ?urey + ?impr ) - @itot
calc qme = ?qmel - @i
calc qmv = ?qmvd - @j

!write out energies to data file
write title unit 22
* @i @5 @force @6 @7 @8 @itot2 @qme @qmv
*

calc i = @i + 15.

if i le 180. goto loop

stop