* generate model B and perform energy surfaces
*

!note that this input does NOT give energy surfaces with the
!lowest energy being zero.  To obtain these, run this script
!first and follow with model_b_surf_all_two.inp, from which
!properly zeroed energy surfaces will be obtained.

set step 15

!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

coor print

!obtain all trans minimum
update cutnb 99.0 ctofnb 98.0 ctonnb 97.0 switch vswitch
mini abnr nstep 200 nprint 10
mini nraph nstep 50

!write coordinates to be read later to initiate
!surfaces 2 and 3
open unit 26 write form name @id1_min.crd
write coor card unit 26
close unit 26

!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

!1)  c1 c2 n3 n4 surface

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

!output coordinates at each step if required
!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

open unit 27 write form name mod1_surf1_final.crd
write coor card unit 27
close unit 27

!2) second surface c2 n3 n4 c5

!read all-trans starting structure

open unit 26 read form name @id1_min.crd
read coor card 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

open unit 27 write form name mod1_surf2_final.crd
write coor card unit 27
close unit 27

!3) third surface n3 n4 c5 c6

!read all-trans structure
open unit 26 read form name @id1_min.crd
read coor card unit 26

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

open unit 27 write form name mod1_surf3_final.crd
write coor card unit 27
close unit 27

stop