* generate model compound B, minimize and calculate 
* minimum interaction energies and distances with
* water for select orientations
*

!set identifier
set id1 model_b
set toppar $CHARMMDATA

set top @toppar/top_all27_prot_lipid.rtf
set par @toppar/par_all27_prot_lipid.prm
set toppar top_mmtsb_example.str

! tolgrd setting: RMS force tolerance to exit minimizer
set a 0.0001

! read the topology
open unit 9 read form name @top
read rtf card unit 9
close unit 9

! read the parameters
open unit 9 read form name @par
read para card unit 9
close unit 9

!stream toppar file with additional rtf/param information
stream @toppar

read sequence cards
* cmpd model 1
*
1
mod1

generate a first none last none setup warn

!create coordinates based on the IC table and minimize the
!structure in the gas phase

ic param
ic seed 1 h11 1 c1 1 c2
ic build

update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch
mini abnr nstep 200 nprint 10
mini nraph nstep 50


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

!QM dipole moment
!hf/6-31g* = 5.7
coor dipole

!output file for dipole moments and interactions with water
!data

open unit 50 write form name @id1_wat.ene
write title unit 50
* model compound 1, charge determination based on
* dipole and interactions with water
* dipole(hf/6-31g*) = 5.69
* empirical  tot=?rdip  x=?xdip  y=?ydip  z=?zdip
*

!The following sections build water molecules in different
!ideal h-bonding orientations around the model compound.
!The minimum interaction energy and distance is then obtained
!and output for comparison with HF/6-31G* values.  Note that
!orientation of the water is built using the IC commands, which
!is similar to a Z-matrix.  Thus, the IC table may be built in 
!a number of ways.  I tend to use extra dummy atoms to allow
!easier visualization of the IC table.

! start interactions with water
!1) O6...HOH

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
*
4
dum dum dum dum

generate dum first none last none setup warn
join a dum renumber

set d 1.70   !starting hbond distance
set u 10000. !variable to identify lowest energy conformation

ic delete sele all end

!create initial ICs that require information from cartesian
!coordinates
ic edit
dihe 1 n3  1 c2  1 o2  3 dum    0.0
dihe 1 c2  1 o2  3 dum 4 dum    0.0
end 
ic fill preserve
!build remaining ICs, include all necessary internal coordinate info
ic edit
dihe 4 dum 3 dum 1 o2  2 h1   180.0
dihe 3 dum 1 o2  2 h1  5 dum    0.0
dihe 1 o2  2 h1  5 dum 6 dum    0.0
dihe 6 dum 5 dum 2 h1  2 oh2  180.0
dihe 5 dum 2 h1  2 oh2 2 h2   180.0
bond 1 o2  3 dum 1.
bond 3 dum 4 dum 1.
bond 1 o2  2 h1  @d  ! distance being optimized
bond 2 h1  5 dum 1.
bond 5 dum 6 dum 1.
bond 2 h1  2 oh2 0.9572
bond 2 oh2 2 h2  0.9572
angl 1 c2  1 o2  3 dum  90.0
angl 1 o2  3 dum 4 dum  90.0
angl 3 dum 1 o2  2 h1   90.0
angl 1 o2  2 h1  5 dum  90.0
angl 2 h1  5 dum 6 dum  90.0
angl 5 dum 2 h1  2 oh2  90.0
angl 2 h1  2 oh2 2 h2  104.52
end

ic print
ic build
coor print

open unit 30 write form name mod1_wat1.pdb
write coor pdb unit 30
* model 1 to water, O2...HOH
*

set m 1000. ! minimum energy

!infinite cutoff to avoid additional updates
update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch

label loop_1

inter -
   sele resn mod1 end -
   sele resn tip3 end

! check for lowest energy and store data
if ?ener lt @m set n @d
if ?ener lt @m set o ?elec
if ?ener lt @m set p ?vdw
if ?ener lt @m set m ?ener

!exit of minimum interaction energy has been passed
if @u lt @m goto done_1
set u @m

!reset hbond length
ic edit
bond 1 o2 2 h1 @d                   ! distance of interest
end

! initialize water coordinates and rebuild with new
! hbond distance
coor init sele resn tip3 .or. resn dum end
ic build

! lower print level to decrease output file size
prnlev -5

incr d by 0.01

if d le 4.0 goto loop_1
label done_1

prnlev 5

!accumulate statistics

calc count = @count + 1
set ai -6.124     ! scaled by 1.16 more neutral compound
set aidist 2.055
calc diff = @m - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

calc distdiff = @n - @aidist

! write ab initio and empirical information to unit 50
write title unit 50
* 1) O2...HOH interaction
* a.i. @ai @aidist
* emp. @m @o @p @n
* ene diff: @diff dist diff: @distdiff
*

!delete dummy atoms and water
delete atom sele resn tip3 .or. resn dum end

!2) N3-H3...OHH (out of plane)

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

set d 1.80   !starting hbond distance
set u 10000. !variable to identify lowest energy conformation

ic delete sele all end

!create initial ICs that require information from cartesian
!coordinates
ic edit
dihe 1 c2  1 n3  1 h3  3 dum    0.0
dihe 1 n3  1 h3  3 dum 4 dum    0.0
end 
ic fill preserve
!build remaining ICs, include all necessary internal coordinate info
ic edit
dihe 4 dum 3 dum 1 h3  2 oh2  180.0
dihe 3 dum 1 h3  2 oh2 2 h1    90.0
dihe 3 dum 1 h3  2 oh2 2 h2   -90.0
bond 1 h3  3 dum 1.
bond 3 dum 4 dum 1.
bond 1 h3  2 oh2 @d  ! distance being optimized
bond 2 oh2 2 h1  0.9572
bond 2 oh2 2 h2  0.9572
angl 1 n3  1 h3  3 dum  90.0
angl 1 h3  3 dum 4 dum  90.0
angl 3 dum 1 h3  2 oh2  90.0
angl 1 h3  2 oh2 2 h1  127.74
angl 1 h3  2 oh2 2 h2  127.74
end

ic print
ic build
coor print

open unit 30 write form name mod1_wat2.pdb
write coor pdb unit 30
* model 1 to water, N3-H3...OHH out of plane
*

set m 1000. ! minimum energy

!infinite cutoff to avoid additional updates
update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch

label loop_2

inter -
   sele resn mod1 end -
   sele resn tip3 end

! check for lowest energy and store data
if ?ener lt @m set n @d
if ?ener lt @m set o ?elec
if ?ener lt @m set p ?vdw
if ?ener lt @m set m ?ener

!exit of minimum interaction energy has been passed
if @u lt @m goto done_2
set u @m

!reset hbond length
ic edit
bond 1 h3 2 oh2 @d                   ! distance of interest
end

! initialize water coordinates and rebuild with new
! hbond distance
coor init sele resn tip3 .or. resn dum end
ic build

! lower print level to decrease output file size
prnlev -5

incr d by 0.01

if d le 2.0 goto loop_2
label done_2

prnlev 5

!accumulate statistics

calc count = @count + 1
set ai -7.267     ! scaled by 1.16 more neutral compound
set aidist 2.118
calc diff = @m - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

calc distdiff = @n - @aidist

! write ab initio and empirical information to unit 50
write title unit 50
* 2) N3-H2...OHH interaction (out of plane)
* a.i. @ai @aidist
* emp. @m @o @p @n
* ene diff: @diff dist diff: @distdiff
*

!delete dummy atoms and water
delete atom sele resn tip3 .or. resn dum end

!3) N5...HOH (out of plane)

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

set d 1.9    !starting hbond distance
set u 10000. !variable to identify lowest energy conformation

ic delete sele all end

!get n3-n4-c5 angle and calculate needed term
quick 7 9 10
calc angn4 = ( 360. - ?thet ) / 2.0

!create initial ICs that require information from cartesian
!coordinates
ic edit
dihe 1 c2  1 n3  1 n4  2 h1     0.0
dihe 1 n3  1 n4  2 h1  3 dum    0.0
dihe 1 n4  2 h1  3 dum 4 dum  180.0
end 
ic fill preserve
!build remaining ICs, include all necessary internal coordinate info
ic edit
dihe 4 dum 3 dum 2 h1  2 oh2    0.0
dihe 3 dum 2 h1  2 oh2 2 h2    90.0  !note out of plane
bond 1 n4  2 h1  @d  ! distance being optimized
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 n3  1 n4  2 h1   @angn4
angl 1 n4  2 h1  3 dum  90.0
angl 2 h1  3 dum 4 dum  90.0
angl 3 dum 2 h1  2 oh2  90.0
angl 2 h1  2 oh2 2 h2  104.52
end

ic print
ic build
coor print

open unit 30 write form name mod1_wat3.pdb
write coor pdb unit 30
* model 1 to water, N4...HOH out of plane
*

set m 1000. ! minimum energy

!infinite cutoff to avoid additional updates
update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch

label loop_3

inter -
   sele resn mod1 end -
   sele resn tip3 end

! check for lowest energy and store data
if ?ener lt @m set n @d
if ?ener lt @m set o ?elec
if ?ener lt @m set p ?vdw
if ?ener lt @m set m ?ener

!exit of minimum interaction energy has been passed
if @u lt @m goto done_3
set u @m

!reset hbond length
ic edit
bond 1 n4 2 h1 @d                   ! distance of interest
end

! initialize water coordinates and rebuild with new
! hbond distance
coor init sele resn tip3 .or. resn dum end
ic build

! lower print level to decrease output file size
prnlev -5

incr d by 0.01

if d le 2.3 goto loop_3
label done_3

prnlev 5

!accumulate statistics

calc count = @count + 1
set ai -5.217     ! scaled by 1.16 more neutral compound
set aidist 2.329
calc diff = @m - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

calc distdiff = @n - @aidist

! write ab initio and empirical information to unit 50
write title unit 50
* 3) N4...HOH interaction
* a.i. @ai @aidist
* emp. @m @o @p @n
* ene diff: @diff dist diff: @distdiff
*

!delete dummy atoms and water
delete atom sele resn tip3 .or. resn dum end

!4) C5-H5...OHH (out of plane)

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

set d 2.2    !starting hbond distance
set u 10000. !variable to identify lowest energy conformation

ic delete sele all end

!create initial ICs that require information from cartesian
!coordinates
ic edit
dihe 1 n4  1 c5  1 h5  3 dum    0.0
dihe 1 c5  1 h5  3 dum 4 dum    0.0
end 
ic fill preserve
!build remaining ICs, include all necessary internal coordinate info
ic edit
dihe 4 dum 3 dum 1 h5  2 oh2  180.0
dihe 3 dum 1 h5  2 oh2 2 h1    90.0
dihe 3 dum 1 h5  2 oh2 2 h2   -90.0
bond 1 h5  3 dum 1.
bond 3 dum 4 dum 1.
bond 1 h5  2 oh2 @d  ! distance being optimized
bond 2 oh2 2 h1  0.9572
bond 2 oh2 2 h2  0.9572
angl 1 c5  1 h5  3 dum  90.0
angl 1 h5  3 dum 4 dum  90.0
angl 3 dum 1 h5  2 oh2  90.0
angl 1 h5  2 oh2 2 h1  127.74
angl 1 h5  2 oh2 2 h2  127.74
end

ic print
ic build
coor print

open unit 30 write form name mod1_wat4.pdb
write coor pdb unit 30
* model 1 to water, C5-H5...OHH out of plane
*

set m 1000. ! minimum energy

!infinite cutoff to avoid additional updates
update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch

label loop_4

inter -
   sele resn mod1 end -
   sele resn tip3 end

! check for lowest energy and store data
if ?ener lt @m set n @d
if ?ener lt @m set o ?elec
if ?ener lt @m set p ?vdw
if ?ener lt @m set m ?ener

!exit of minimum interaction energy has been passed
if @u lt @m goto done_4
set u @m

!reset hbond length
ic edit
bond 1 h5 2 oh2 @d                   ! distance of interest
end

! initialize water coordinates and rebuild with new
! hbond distance
coor init sele resn tip3 .or. resn dum end
ic build

! lower print level to decrease output file size
prnlev -5

incr d by 0.01

if d le 2.7 goto loop_4
label done_4

prnlev 5

!accumulate statistics

calc count = @count + 1
set ai -3.858     ! scaled by 1.16 more neutral compound
set aidist 2.460
calc diff = @m - @ai
calc sum = @sum + @diff
calc diff2 = abs(@diff*@diff)
calc sum2 = @sum2 + @diff2
calc asum = @asum + abs(@diff)

calc distdiff = @n - @aidist

! write ab initio and empirical information to unit 50
write title unit 50
* 4) C5-H2...OHH interaction (out of plane)
* a.i. @ai @aidist
* emp. @m @o @p @n
* ene diff: @diff dist diff: @distdiff
*

!delete dummy atoms and water
delete atom sele resn tip3 .or. resn dum 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 50
* avediff, rmsdiff, average absolute error
* @avediff @rmsdif @aae
*

stop