* 4-ethylimidazole minimization and torsional energy surface
* 

fast on

! output identifier                                                  
set id1 imid_ethyl_1

! data output file for energies and forces
open unit 90 write card name @id1_surf.map
! data output file for selected angles
open unit 91 write card name @id1_surf.ang

! 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_all22_model.inp
set b par_all22_prot.inp
!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 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 
read sequence card
* 4-ethylimidazole
*  
1
eimi

generate a first none last none setup warn

! autogenerate the ic table and put into topology (reformat as required)
ic generate

!fill ic table with bond lengths and angles from parameter file
ic para

!build molecule
ic seed 1 ha1 1 ca 1 cb
ic build

coor print

!assign atom numbers to variable for quick commands
coor stat sele type cg end
set cg ?selatom
coor stat sele type cd2 end
set cd2 ?selatom
coor stat sele type hd2 end
set hd2 ?selatom
coor stat sele type nd1 end
set nd1 ?selatom
coor stat sele type hd1 end
set hd1 ?selatom
coor stat sele type ce1 end
set ce1 ?selatom
coor stat sele type he1 end
set he1 ?selatom
coor stat sele type ne2 end
set ne2 ?selatom
coor stat sele type cb end
set cb ?selatom
coor stat sele type ca end
set ca ?selatom

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

!minimization

!apply dihedral constraint to get close to global minima
!and then perform full minimization

cons dihe a 1 ca a 1 cb a 1 cg a 1 nd1 min 60. force 10000.
mini abnr nstep 200 nprint 50

!clear dihedral constraint
cons cldh
mini abnr nstep 50 nprint 10
mini nraph nstep 20 nprint 5 tolgrd 0.0001

!set global minimum energies to allow for
!offset of energy surface 0.0
set etotmin ?ener
set elecmin ?elec
set vdwmin ?vdw
set dihemin ?dihe
set bondmin ?bond
set anglmin ?angl
set ureymin ?urey
set imprmin ?impr
calc intmin = ?bond + ?angl + ?urey + ?impr

!obtain value of dihedral to be investigated
quick 1 8 7 5

coor print

!force dihedral to 0.0, which is the value at
!which the surface will be initiated.
!cons dihe a 1 ca a 1 cb a 1 cg a 1 nd1 min 90. force 10000.
!mini abnr nstep 200 nprint 50
!cons cldh
cons dihe a 1 ca a 1 cb a 1 cg a 1 nd1 min 0. force 10000.
mini abnr nstep 200 nprint 50
cons cldh

set dihe 0.0

label loop

!constraint to force dihedral to desired value
cons dihe a 1 ca a 1 cb a 1 cg a 1 nd1 min @dihe force 10000.

!minimization
mini abnr nstep 50 nprint 10
mini nraph nstep 20 nprint 5

!save value of gradient for output below
set force ?grms

!clear dihedral constraint and get energy
cons cldh
energy

!offset energy terms
calc etot = ?ener - @etotmin
calc eele = ?elec - @elecmin
calc evdw = ?vdw - @vdwmin
calc edih = ?dihe - @dihemin
calc eint = ( ?bond + ?angl + ?urey + ?impr ) - @intmin

!write offset energies and force to .map file
write title unit 90
* @dihe @etot @eele @evdw @edih @eint @force
*

!obtain selected angles and write to unit 91
quick @ca @cb @cg
set ang1 ?thet
quick @cb @cg @nd1
set ang2 ?thet
quick @cb @cg @cd2
set ang3 ?thet

write title unit 91
* @dihe @ang1 @ang2 @ang3
*

calc dihe = @dihe + 15.

if dihe le 360. goto loop

stop