* imidazole minimization and imidazole-water
* interaction surfaces
* 

fast on

! output identifier                                                  
set id1 imid_1

open unit 90 write card name @id1_wat.ene

! 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
* imidazole
*  
1
imia

generate a first none last none setup warn

read coor card
* imidazole
*
    9
    1    1 IMIA CG     1.41005  -0.09516   0.01975 A    1      0.00000
    2    1 IMIA HG     1.93888  -1.03723   0.04795 A    1      0.00000
    3    1 IMIA CD2    1.89866   1.17143  -0.03758 A    1      0.00000
    4    1 IMIA HD2    2.92468   1.51776  -0.06770 A    1      0.00000
    5    1 IMIA ND1    0.03669   0.03098   0.03706 A    1      0.00000
    6    1 IMIA HD1   -0.62870  -0.71351   0.07688 A    1      0.00000
    7    1 IMIA CE1   -0.23842   1.35741  -0.01024 A    1      0.00000
    8    1 IMIA HE1   -1.25075   1.76593  -0.00980 A    1      0.00000
    9    1 IMIA NE2    0.85865   2.08270  -0.05630 A    1      0.00000


!assign atom numbers to variable for quick commands
coor stat sele type cg end
set cg ?selatom
coor stat sele type hg end
set hg ?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

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

!fully minimize structure to insure that minimum
!energy structure is used in the interactions
!with water

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

!calculate dipole moment
!experiment: 3.67
!note that c22 dipole moments are typically larger than experiment
!due to the implicit inclusion of electronic polarization in the
!force field
coor dipole sele all end

write title unit 90
* imidazole
* experimental dipole moment
* 3.67
* empirical dipole moment
* X=  ?xdip Y=  ?ydip Z=  ?zdip Tot=   ?rdip
*______________________________________________
*imidazole-water interactions
*

!imidazole-water interaction surfaces
!

! variables for averages etc.
set count 0
set asum 0.0
set sum 0.0
set sum2 0.0

!calculaiton "outside" ce1-ne2-cd2 angle for later use
quick @ce1 @ne2 @cd2
calc angne2 = ( 360. - ?thet ) / 2.0

! calculate "outside" angle for tip3 water
calc angwat = (360.0-104.52)/2.0

! delete all IC tables before starting
ic delete sele all end

!1) Ce1-He1...OHH

!adjust d and maximum distance to sample
!only the minimum energy region of the
!surface, thereby saving cpu time.

set d 2.20  ! initial interaction distance

read sequence card
* water
*
1
tip3

generate wat first none last none setup warn noangle nodihedral

join a wat renumber

read sequence card
* dummy residue
*
2
dum dum

generate dum first none last none setup warn

join a dum renumber

! build initial IC's that must get info from
! cartesian coordinates
!residue numbers
!1: imidazole
!2: water
!3: dummy
!4: dummy
ic edit
dihe 1 ne2 1 ce1 1 he1 3 dum   0.0
dihe 1 ce1 1 he1 3 dum 4 dum   0.0
end
ic fill preserve
ic print
ic edit
dihe 4 dum 3 dum 1 he1 2 oh2 180.0
dihe 3 dum 1 he1 2 oh2 2 h1    0.0
dihe 3 dum 1 he1 2 oh2 2 h2  180.0
bond 1 he1 3 dum 1.0
bond 3 dum 4 dum 1.0
bond 1 he1 2 oh2 @d   !distance being varied
bond 2 oh2 2 h1 0.9572
bond 2 oh2 2 h2 0.9572
angl 1 ce1 1 he1 3 dum 90.0
angl 1 he1 3 dum 4 dum 90.0
angl 3 dum 1 he1 2 oh2 90.0
angl 1 he1 2 oh2 2 h1  @angwat
angl 1 he1 2 oh2 2 h2  @angwat
end

ic print
coor print
ic build
coor print

!write initial coordinates to allow for 
!visual inspection of structure
open unit 30 write form name imid_wat_he1.pdb
write coor pdb unit 30

!create nonbond list, use of infinite cutoff is essential
update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

set m 1000. ! required for determination minimum energy geometry

!uncomment to write out energy surfaces
!open unit 50 write form name @id1_wat_he1.surf

label loop_1

!set distance prior to build
ic edit
bond 1 he1 2 oh2 @d
end

!build the cartesian coordinates
ic build

!calculate interaction energy (this command avoids a nonbond list
!update, thereby saving cpu time)
inter -
   sele resn imia end -
   sele resn tip3 end

!uncomment to write out energy surfaces
!write title unit 50
!* @d ?ener ?elec ?vdw
!*

! check for lowest energy and save data
if ?ener le @m then set emin ?ener
if ?ener le @m then set dmin @d
if ?ener le @m then set elem ?elec
if ?ener le @m then set vdwm ?vdw
if ?ener le @m then set m ?ener

!initialize water and dummy atom coordinates to
!allow them to be rebuilt
coor init sele resn tip3 .or. resn dum show end

!increase distance 
incr d by 0.01

if d le 2.50 goto loop_1

calc count = @count + 1
set ai -2.44
set aidist 2.45
calc diff = @emin - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

write title unit 90
* 1) Ce1-He1...OHH
* a.i: @aidist, @ai, ai-emp diff:@diff
* emp: @dmin @emin @elem @vdwm
*

! remove old water and dummy atoms and associated IC's
delete atom sele resn tip3 .or. resn dum end
ic delete sele all end

!2) Nd1-Hd1...OHH

!adjust d and maximum distance to sample
!only the minimum energy region of the
!surface, thereby saving cpu time.

set d 1.80  ! initial interaction distance

read sequence card
* water
*
1
tip3

generate wat first none last none setup warn noangle nodihedral

join a wat renumber

read sequence card
* dummy residue
*
2
dum dum

generate dum first none last none setup warn

join a dum renumber

! build initial IC's that must get info from
! cartesian coordinates
!residue numbers
!1: imidazole
!2: water
!3: dummy
!4: dummy
ic edit
dihe 1 ce1 1 nd1 1 hd1 3 dum   0.0
dihe 1 nd1 1 hd1 3 dum 4 dum   0.0
end
ic fill preserve
ic print
ic edit
dihe 4 dum 3 dum 1 hd1 2 oh2 180.0
dihe 3 dum 1 hd1 2 oh2 2 h1    0.0
dihe 3 dum 1 hd1 2 oh2 2 h2  180.0
bond 1 hd1 3 dum 1.0
bond 3 dum 4 dum 1.0
bond 1 hd1 2 oh2 @d   !distance being varied
bond 2 oh2 2 h1 0.9572
bond 2 oh2 2 h2 0.9572
angl 1 nd1 1 hd1 3 dum 90.0
angl 1 hd1 3 dum 4 dum 90.0
angl 3 dum 1 hd1 2 oh2 90.0
angl 1 hd1 2 oh2 2 h1  @angwat
angl 1 hd1 2 oh2 2 h2  @angwat
end

ic print
coor print
ic build
coor print

!write initial coordinates to allow for 
!visual inspection of structure
open unit 30 write form name imid_wat_hd1.pdb
write coor pdb unit 30

!create nonbond list, use of infinite cutoff is essential
update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

set m 1000. ! required for determination minimum energy geometry

!uncomment to write out energy surfaces
!open unit 50 write form name @id1_wat_hd1.surf

label loop_2

!set distance prior to build
ic edit
bond 1 hd1 2 oh2 @d
end

!build the cartesian coordinates
ic build

!calculate interaction energy (this command avoids a nonbond list
!update, thereby saving cpu time)
inter -
   sele resn imia end -
   sele resn tip3 end

!uncomment to write out energy surfaces
!write title unit 50
!* @d ?ener ?elec ?vdw
!*

! check for lowest energy and save data
if ?ener le @m then set emin ?ener
if ?ener le @m then set dmin @d
if ?ener le @m then set elem ?elec
if ?ener le @m then set vdwm ?vdw
if ?ener le @m then set m ?ener

!initialize water and dummy atom coordinates to
!allow them to be rebuilt
coor init sele resn tip3 .or. resn dum show end

!increase distance 
incr d by 0.01

if d le 2.00 goto loop_2

calc count = @count + 1
set ai -6.65
set aidist 2.07
calc diff = @emin - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

write title unit 90
* 2) Nd1-Hd1...OHH
* a.i: @aidist, @ai, ai-emp diff:@diff
* emp: @dmin @emin @elem @vdwm
*

! remove old water and dummy atoms and associated IC's
delete atom sele resn tip3 .or. resn dum end
ic delete sele all end

!3) Cg-Hg...OHH

!adjust d and maximum distance to sample
!only the minimum energy region of the
!surface, thereby saving cpu time.

set d 2.40  ! initial interaction distance

read sequence card
* water
*
1
tip3

generate wat first none last none setup warn noangle nodihedral

join a wat renumber

read sequence card
* dummy residue
*
2
dum dum

generate dum first none last none setup warn

join a dum renumber

! build initial IC's that must get info from
! cartesian coordinates
!residue numbers
!1: imidazole
!2: water
!3: dummy
!4: dummy
ic edit
dihe 1 nd1 1 cg  1 hg  3 dum   0.0
dihe 1 cg  1 hg  3 dum 4 dum   0.0
end
ic fill preserve
ic print
ic edit
dihe 4 dum 3 dum 1 hg  2 oh2 180.0
dihe 3 dum 1 hg  2 oh2 2 h1    0.0
dihe 3 dum 1 hg  2 oh2 2 h2  180.0
bond 1 hg  3 dum 1.0
bond 3 dum 4 dum 1.0
bond 1 hg  2 oh2 @d   !distance being varied
bond 2 oh2 2 h1 0.9572
bond 2 oh2 2 h2 0.9572
angl 1 cg  1 hg  3 dum 90.0
angl 1 hg  3 dum 4 dum 90.0
angl 3 dum 1 hg  2 oh2 90.0
angl 1 hg  2 oh2 2 h1  @angwat
angl 1 hg  2 oh2 2 h2  @angwat
end

ic print
coor print
ic build
coor print

!write initial coordinates to allow for 
!visual inspection of structure
open unit 30 write form name imid_wat_hg.pdb
write coor pdb unit 30

!create nonbond list, use of infinite cutoff is essential
update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

set m 1000. ! required for determination minimum energy geometry

!uncomment to write out energy surfaces
!open unit 50 write form name @id1_wat_hg.surf

label loop_3

!set distance prior to build
ic edit
bond 1 hg  2 oh2 @d
end

!build the cartesian coordinates
ic build

!calculate interaction energy (this command avoids a nonbond list
!update, thereby saving cpu time)
inter -
   sele resn imia end -
   sele resn tip3 end

!uncomment to write out energy surfaces
!write title unit 50
!* @d ?ener ?elec ?vdw
!*

! check for lowest energy and save data
if ?ener le @m then set emin ?ener
if ?ener le @m then set dmin @d
if ?ener le @m then set elem ?elec
if ?ener le @m then set vdwm ?vdw
if ?ener le @m then set m ?ener

!initialize water and dummy atom coordinates to
!allow them to be rebuilt
coor init sele resn tip3 .or. resn dum show end

!increase distance 
incr d by 0.01

if d le 2.70 goto loop_3

calc count = @count + 1
set ai -2.40
set aidist 2.51
calc diff = @emin - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

write title unit 90
* 3) Cg-Hg...OHH
* a.i: @aidist, @ai, ai-emp diff:@diff
* emp: @dmin @emin @elem @vdwm
*

! remove old water and dummy atoms and associated IC's
delete atom sele resn tip3 .or. resn dum end
ic delete sele all end

!4) Cd2-Hd2...OHH

!adjust d and maximum distance to sample
!only the minimum energy region of the
!surface, thereby saving cpu time.

set d 2.40  ! initial interaction distance

read sequence card
* water
*
1
tip3

generate wat first none last none setup warn noangle nodihedral

join a wat renumber

read sequence card
* dummy residue
*
2
dum dum

generate dum first none last none setup warn

join a dum renumber

! build initial IC's that must get info from
! cartesian coordinates
!residue numbers
!1: imidazole
!2: water
!3: dummy
!4: dummy
ic edit
dihe 1 ne2 1 cd2 1 hd2 3 dum   0.0
dihe 1 cd2 1 hd2 3 dum 4 dum   0.0
end
ic fill preserve
ic print
ic edit
dihe 4 dum 3 dum 1 hd2 2 oh2 180.0
dihe 3 dum 1 hd2 2 oh2 2 h1    0.0
dihe 3 dum 1 hd2 2 oh2 2 h2  180.0
bond 1 hd2 3 dum 1.0
bond 3 dum 4 dum 1.0
bond 1 hd2 2 oh2 @d   !distance being varied
bond 2 oh2 2 h1 0.9572
bond 2 oh2 2 h2 0.9572
angl 1 cd2 1 hd2 3 dum 90.0
angl 1 hd2 3 dum 4 dum 90.0
angl 3 dum 1 hd2 2 oh2 90.0
angl 1 hd2 2 oh2 2 h1  @angwat
angl 1 hd2 2 oh2 2 h2  @angwat
end

ic print
coor print
ic build
coor print

!write initial coordinates to allow for 
!visual inspection of structure
open unit 30 write form name imid_wat_hd2.pdb
write coor pdb unit 30

!create nonbond list, use of infinite cutoff is essential
update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

set m 1000. ! required for determination minimum energy geometry

!uncomment to write out energy surfaces
!open unit 50 write form name @id1_wat_hd2.surf

label loop_4

!set distance prior to build
ic edit
bond 1 hd2 2 oh2 @d
end

!build the cartesian coordinates
ic build

!calculate interaction energy (this command avoids a nonbond list
!update, thereby saving cpu time)
inter -
   sele resn imia end -
   sele resn tip3 end

!uncomment to write out energy surfaces
!write title unit 50
!* @d ?ener ?elec ?vdw
!*

! check for lowest energy and save data
if ?ener le @m then set emin ?ener
if ?ener le @m then set dmin @d
if ?ener le @m then set elem ?elec
if ?ener le @m then set vdwm ?vdw
if ?ener le @m then set m ?ener

!initialize water and dummy atom coordinates to
!allow them to be rebuilt
coor init sele resn tip3 .or. resn dum show end

!increase distance 
incr d by 0.01

if d le 2.70 goto loop_4

calc count = @count + 1
set ai -0.99
set aidist 2.66
calc diff = @emin - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

write title unit 90
* 4) Cd2-Hd2...OHH
* a.i: @aidist, @ai, ai-emp diff:@diff
* emp: @dmin @emin @elem @vdwm
*

! remove old water and dummy atoms and associated IC's
delete atom sele resn tip3 .or. resn dum end
ic delete sele all end

!5) Ne2...HoH

!adjust d and maximum distance to sample
!only the minimum energy region of the
!surface, thereby saving cpu time.

set d 1.80  ! initial interaction distance

read sequence card
* water
*
1
tip3

generate wat first none last none setup warn noangle nodihedral

join a wat renumber

read sequence card
* dummy residue
*
2
dum dum

generate dum first none last none setup warn

join a dum renumber

! build initial IC's that must get info from
! cartesian coordinates
!residue numbers
!1: imidazole
!2: water
!3: dummy
!4: dummy
ic edit
dihe 1 nd1 1 ce1 1 ne2 2 h1  180.0
dihe 1 ce1 1 ne2 2 h1  3 dum 180.0
dihe 1 ne2 2 h1  3 dum 4 dum 180.0
end
ic fill preserve
ic print
ic edit
dihe 4 dum 3 dum 2 h1  2 oh2   0.0
dihe 3 dum 2 h1  2 oh2 2 h2    0.0
bond 1 ne2 2 h1  @d   !distance being varied
bond 2 h1  3 dum 1.
bond 3 dum 4 dum 1.
bond 2 h1  2 oh2 0.9572
bond 2 oh2 2 h2  0.9572
angl 1 ce1 1 ne2 2 h1  @angne2
angl 1 ne2 2 h1  3 dum  90.
angl 2 h1  3 dum 4 dum  90.
angl 3 dum 2 h1  2 oh2  90.
angl 2 h1  2 oh2 2 h2  104.52
end

ic print
coor print
ic build
coor print

!write initial coordinates to allow for 
!visual inspection of structure
open unit 30 write form name imid_wat_ne2.pdb
write coor pdb unit 30

!create nonbond list, use of infinite cutoff is essential
update inbfrq -1 ihbfrq 0 -
@7 @8 @9 vfswitch cutnb @4 ctofnb @6 ctonnb @5

set m 1000. ! required for determination minimum energy geometry

!uncomment to write out energy surfaces
!open unit 50 write form name @id1_wat_ne2.surf

label loop_5

!set distance prior to build
ic edit
bond 1 ne2 2 h1 @d
end

!build the cartesian coordinates
ic build

!calculate interaction energy (this command avoids a nonbond list
!update, thereby saving cpu time)
inter -
   sele resn imia end -
   sele resn tip3 end

!uncomment to write out energy surfaces
!write title unit 50
!* @d ?ener ?elec ?vdw
!*

! check for lowest energy and save data
if ?ener le @m then set emin ?ener
if ?ener le @m then set dmin @d
if ?ener le @m then set elem ?elec
if ?ener le @m then set vdwm ?vdw
if ?ener le @m then set m ?ener

!initialize water and dummy atom coordinates to
!allow them to be rebuilt
coor init sele resn tip3 .or. resn dum show end

!increase distance 
incr d by 0.01

if d le 2.00 goto loop_5

calc count = @count + 1
set ai -7.26
set aidist 2.12
calc diff = @emin - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

write title unit 90
* 1) Ne2...HOH
* a.i: @aidist, @ai, ai-emp diff:@diff
* emp: @dmin @emin @elem @vdwm
*

! remove old water and dummy atoms and associated IC's
delete atom sele resn tip3 .or. resn dum end
ic delete sele all end

!calculate statistics
! average difference
calc avediff = @sum/@count
calc avediff2 = @sum2/@count
calc rmsdif = sqrt(@avediff2 -(@avediff*@avediff))
!AAE; average absolute error
calc aae = @asum/@count

write title unit 90
* avediff, rmsdiff, average absolute error
* @avediff @rmsdif @aae
*

stop