* generate the full drug (model compound 2) and
* perform energy surfaces
*

set step 15.

!set identifier
set id1 model_2
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.00001

! 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 2
*
1
mod2

generate a first none last none setup warn

!generate preliminary IC table
!ic generate
!stop

ic param
ic seed 1 cg  1 cd1 1 ne1
ic build
ic print

coor print

open unit 26 write form name @id1_gen.pdb
write coor pdb unit 26

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

open unit 26 write form name @id1_min.pdb
write coor pdb unit 26

ic fill
ic print

!set minimum values to offset surface to 0.0
set totmin ?ener
set elemin ?elec
set vdwmin ?vdw
set dihmin ?dihe
calc intmin = ?bond + ?angl + ?urey + ?impr

!get dihe a 1 ce2 a 1 cz2 a 1 c2  a 1 n3 = 0 starting structure
!for surface 1

coor init

ic edit
dihe 1 CE2 1 CZ2 1 C2 1 N3 10.0
end

ic print 
ic seed 1 cg  1 cd1 1 ne1
ic build
coor print

!initial minimization to get correct starting structure. Note
!use of 10 deg. in IC edit to "bump" the structure out of a local
!minimum

cons dihe a 1 ce2 a 1 cz2 a 1 c2  a 1 n3   force 10000.0 min 0.0
cons dihe a 1 cz2 a 1 c2  a 1 n3  a 1 n4   force 10000.0 min 180.0
cons dihe a 1 c2  a 1 n3  a 1 n4  a 1 c5   force 10000.0 min 180.0
cons dihe a 1 n3  a 1 n4  a 1 c5  a 1 cpg  force 10000.0 min 180.0
cons dihe a 1 n4  a 1 c5  a 1 cpg a 1 cpd1 force 10000.0 min 180.0

mini SD nstep 100
mini conj nstep 100
cons cldh

open unit 30 write form name @id1_surf1.map

set dihecons 0.

label loop_1
cons dihe a 1 ce2 a 1 cz2 a 1 c2  a 1 n3   force 10000.0 min @dihecons
cons dihe a 1 cz2 a 1 c2  a 1 n3  a 1 n4   force 10000.0 min 180.0
cons dihe a 1 c2  a 1 n3  a 1 n4  a 1 c5   force 10000.0 min 180.0
cons dihe a 1 n3  a 1 n4  a 1 c5  a 1 cpg  force 10000.0 min 180.0
cons dihe a 1 n4  a 1 c5  a 1 cpg a 1 cpd1 force 10000.0 min 180.0

mini conj nstep 100 inbfrq 0 inbfrq 0
mini abnr nstep 100 inbfrq 0 inbfrq 0
mini nraph tolstp 0.000001 tolgrd 0.000001 -
nstep 50 nprint 100 inbfrq 0 inbfrq 0

cons cldh
energy

calc tote = ?ener - @totmin
calc elee = ?elec - @elemin
calc vdwe = ?vdw - @vdwmin
calc dihe = ?dihe - @dihmin
calc inte = ( ?bond + ?angl + ?urey + ?impr ) - @intmin

write title unit 30
* @dihecons @tote @elee @vdwe @dihe @inte
*

ic fill
ic print

!open unit 26 write form name mod1_surf1_@dihecons.pdb
!write coor pdb unit 26

calc dihecons = @dihecons + @step

!if @dihecons le 180 goto loop_1
if @dihecons le 360 goto loop_1

!read in starting coordinates
open unit 26 read form name @id1_min.pdb
read coor pdb unit 26

open unit 30 write form name @id1_surf2.map

set dihecons 180.

label loop_2
cons dihe a 1 ce2 a 1 cz2 a 1 c2  a 1 n3   force 10000.0 min 180.0
cons dihe a 1 cz2 a 1 c2  a 1 n3  a 1 n4   force 10000.0 min 180.0
cons dihe a 1 c2  a 1 n3  a 1 n4  a 1 c5   force 10000.0 min 180.0
cons dihe a 1 n3  a 1 n4  a 1 c5  a 1 cpg  force 10000.0 min 180.0
cons dihe a 1 n4  a 1 c5  a 1 cpg a 1 cpd1 force 10000.0 min @dihecons

mini conj nstep 100 inbfrq 0 inbfrq 0
mini abnr nstep 100 inbfrq 0 inbfrq 0
mini nraph tolstp 0.000001 tolgrd 0.000001 -
nstep 50 nprint 100 inbfrq 0 inbfrq 0

cons cldh
energy

calc tote = ?ener - @totmin
calc elee = ?elec - @elemin
calc vdwe = ?vdw - @vdwmin
calc dihe = ?dihe - @dihmin
calc inte = ( ?bond + ?angl + ?urey + ?impr ) - @intmin

write title unit 30
* @dihecons @tote @elee @vdwe @dihe @inte
*

ic fill
ic print

!open unit 26 write form name mod2_surf1_@dihecons.pdb
!write coor pdb unit 26

calc dihecons = @dihecons - @step

if @dihecons ge 0 goto loop_2

stop

!generate and minimize global energy minimum
coor init

ic edit
dihe 1 c1 1 c2 1 n3 1 n4   0.0
dihe 1 c2 1 n3 1 n4 1 c5 180.0
dihe 1 n3 1 n4 1 c5 1 c6 180.0
end

ic seed 1 h11 1 c1 1 c2
ic build
ic print

coor print

open unit 26 write form name @id1_gen_2.pdb
write coor pdb unit 26

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

open unit 26 write form name @id1_min_2.pdb
write coor pdb unit 26

ic fill
ic print

!generate and minimize global energy minimum
coor init

ic edit
dihe 1 c1 1 c2 1 n3 1 n4 180.0
dihe 1 c2 1 n3 1 n4 1 c5   0.0
dihe 1 n3 1 n4 1 c5 1 c6 180.0
end

ic seed 1 h11 1 c1 1 c2
ic build

coor print

open unit 26 write form name @id1_gen_3.pdb
write coor pdb unit 26

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

open unit 26 write form name @id1_min_3.pdb
write coor pdb unit 26

ic fill
ic print

!generate and minimize global energy minimum
coor init

ic edit
dihe 1 c1 1 c2 1 n3 1 n4   0.0
dihe 1 c2 1 n3 1 n4 1 c5   0.0
dihe 1 n3 1 n4 1 c5 1 c6 180.0
end

ic seed 1 h11 1 c1 1 c2
ic build

coor print

open unit 26 write form name @id1_gen_4.pdb
write coor pdb unit 26

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

open unit 26 write form name @id1_min_4.pdb
write coor pdb unit 26

ic fill
ic print

stop

!perform surface
open unit 26 write form name pacp_@id1.surf

!set minimum values to offset surface to 0.0

set a ?ener
set b ?elec
set c ?vdw
set d ?dihe
set e ?angl
set f ?bond
set g ?urey

set 1 180

if stop eq 1 then STOP

label loop_pacp
!CONS DIHE A 1 c23 A 1 C22 A 1 N21 A 1 H211 FORCE 10000.0 MIN 120.0
CONS DIHE A 1 c23 A 1 c22 A 1 N21 A 1 C15 FORCE 10000.0 MIN @1.0

mini abnr nstep 50 inbfrq 0 inbfrq 0
mini abnr nstep 50 inbfrq 0 inbfrq 0
mini nraph tolstp 0.000001 tolgrd 0.000001 -
nstep 500 nprint 100 inbfrq 0 inbfrq 0

set 5 ?ener
decr 5 by ?cdih
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
set 9 ?angl
decr 9 by @e
set 0 ?bond
decr 0 by @f
set 4 ?urey
decr 4 by @g

calc inter = @9 + @0 + @4
write title unit 26
* @1.0 @5 @6 @7 @8 @inter
*

!* @1.0 @5 @6 @7 @8 @9 @0 @4
!*

cons cldh

ic fill
ic print

!open unit 26 write form name pacp.@1.pdb
!write coor pdb unit 26

decr 1 by 10
if 1 ge 0 goto loop_pacp

close unit 26

stop