* Add dummy hydrogens and minimize structure
* Usage: charmm pdbname= < add_h_min.inp

!set pdbname = cpept
set segid = @pdbname
set topfile = toppar/top_all22_prot_cmap_phmd.inp
set parfile = toppar/par_all22_prot_cmap_phmd.inp
set gbsetup = toppar/setup-gbsw.str

! Read topology and parameter files
! ---------------------------------
open unit 1 read form name @topfile
read rtf card unit 1
close unit 1

open unit 1 read form name @parfile
read param card unit 1
close unit 1

! Read sequence
! -------------
open unit 1 read form name struct/@{pdbname}.pdb
read sequ pdb unit 1 

generate @segid first ACE last CT2 setup

! Read coordinates
!--------------------
open unit 1 read form name struct/@{pdbname}.pdb
read coor pdb unit 1

ic para
ic build

! Build hydrogens
! --------------
hbuild sele hydrogen end

! Patch dummy hydrogen atoms for ASP, GLU, and C-terminus
! -------------------------------------------------------
set segid = ?selsegi
set res1 = ?selresi         ! get the starting value of resid
calc reslast = ?selresi + ?nres -1

set ires @res1

label patch
 define residue sele segid @segid .and. resid @ires show end
 if ?selresn eq asp patch aspp2 @segid @ires
 if ?selresn eq glu patch glup2 @segid @ires
 incr ires
if ires le @reslast goto patch

! patch C-terminal carboxyl group
!patch ctrp2 @segid @reslast

autogen angles dihed

hbuild sele segi @segid .and. resn asp .and. type hd* end
hbuild sele segi @segid .and. resn glu .and. type he* end 

! add dummy hydrogens to C-ter carboxyl if needed
!hbuild sele segi @segid .and. resid @reslast .and. type hc* end 

! Edit dihedrals for Glupp2 and Aspp2 to syn- positions
! -----------------------------------------------------
ic gene
ic fill

set ires @res1

label edit
 define residue sele segid @segid .and. resid @ires show end
 if ?selresn eq asp then
    ic edit
       dihe @segid @ires cb @segid @ires cg @segid @ires od1 @segid @ires hd1 180.0
       dihe @segid @ires cb @segid @ires cg @segid @ires od2 @segid @ires hd2 180.0
    end
    coor init sele segi @segid .and. resid @ires .and. type hd* end
 endif
 if ?selresn eq glu then
    ic edit
       dihe @segid @ires cg @segid @ires cd @segid @ires oe1 @segid @ires he1 180.0
       dihe @segid @ires cg @segid @ires cd @segid @ires oe2 @segid @ires he2 180.0
    end
    coor init sele segi @segid .and. resid @ires .and. type he* end
 endif
 incr ires
if ires le @reslast goto edit

! These are for the C-terminal carboxyl group
!ic edit
!   dihe @segid @reslast ca @segid @reslast c @segid @reslast ot1 @segid @reslast hc1 180.0
!   dihe @segid @reslast ca @segid @reslast c @segid @reslast ot2 @segid @reslast hc2 180.0
!end
!coor init sele segid @segid .and. resi @reslast .and. type hc* end

ic build

! GB options
!----------
set ctonnb = 16
set ctofnb = 16
set cutnb = 20

stream @gbsetup

scalar wmain statistics select all end
define check select (.not type H* ) .and. ( property wmain .eq. 0.0 ) show end
if ?nsel ne 0  stop       !some heavy atom have a zero radius

set para atom switch cdie vdw vswitch -
         ctonnb @ctonnb ctofnb @ctofnb cutnb @cutnb

NBOND @para

define heavy sele .not. type H* end

cons harm force 5.0 sele heavy end

mini   sd nstep  50 nprint 10 step 0.005
mini abnr nstep  50 nprint 10 step 0.005

! Write coordinate file
! ---------------------

open unit 1 write card name struct/@{pdbname}_h.pdb
write coor pdb unit 1
close unit 1

open unit 1 card write name struct/@{pdbname}_h.psf
write psf card unit 1
close unit 1

open unit 1 card write name struct/@{pdbname}_h.crd
write coor card unit 1
close unit 1

stop