* generate model B and perform energy surfaces
*

!run model_one_surf_all_one.inp to create the minimum energy
!structures for each surface; these are read into this
!input to properly zero the surfaces

set step 5

!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.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 1
*
1
mod1

generate a first none last none setup warn

!generate preliminary IC table
!ic generate
!stop

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

!1) c1 c2 n3 n4 surface

!read previous minimum energy coordinates
open unit 27 read form name mod1_surf1_final.crd
read coor card unit 27
close unit 27

cons dihe a 1 c1  a 1 c2  a 1 n3  a 1 n4  force 10000.0 min 0.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 c6  force 10000.0 min 180.0

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

!check dihedrals
!c1-c2-n3-n4
quick 1 5 7 9
!c2-n3-n4-c5
quick 5 7 9 10
!n3-n4-c5-c6
quick 7 9 10 12

!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

!read the all-trans starting geometry
open unit 26 write form name @id1_min.crd
write coor card unit 26
close unit 26

open unit 30 write form name @id1_surf1.map

set dihecons 180.

label loop_1
cons dihe a 1 c1  a 1 c2  a 1 n3  a 1 n4  force 10000.0 min @dihecons
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 c6  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 ge 0 goto loop_1

!2) c2 n3 n4 c5 surface

!read previous minimum energy coordinates
open unit 27 read form name mod1_surf2_final.crd
read coor card unit 27
close unit 27

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

!check dihedrals
!c1-c2-n3-n4
quick 1 5 7 9
!c2-n3-n4-c5
quick 5 7 9 10
!n3-n4-c5-c6
quick 7 9 10 12

!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

!read the all-trans starting geometry
open unit 26 write form name @id1_min.crd
write coor card unit 26
close unit 26

open unit 30 write form name @id1_surf2.map

set dihecons 180.

label loop_2
cons dihe a 1 c1  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 @dihecons
cons dihe a 1 n3  a 1 n4  a 1 c5  a 1 c6  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_surf2_@dihecons.pdb
!write coor pdb unit 26

calc dihecons = @dihecons - @step

if @dihecons ge 0 goto loop_2

!3) n3 n4 c5 c6 surface

!read the all-trans starting geometry
!this is close to the minimum structure for surface 3
open unit 26 write form name @id1_min.crd
write coor card unit 26
close unit 26

cons dihe a 1 c1  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 c6  force 10000.0 min 180.0

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

!check dihedrals
!c1-c2-n3-n4
quick 1 5 7 9
!c2-n3-n4-c5
quick 5 7 9 10
!n3-n4-c5-c6
quick 7 9 10 12

!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

open unit 30 write form name @id1_surf3.map

set dihecons 180.

label loop_3
cons dihe a 1 c1  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 c6  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 mod1_surf2_@dihecons.pdb
!write coor pdb unit 26

calc dihecons = @dihecons - @step

if @dihecons ge 0 goto loop_3

stop