* Constant pH Langevin simulation for C-peptide at pH 4 with ionic strength of 150 mM
* Titrating residues: Glu, His, Lys, Arg

set topfile = toppar/top_all22_prot_cmap_phmd.inp
set parfile = toppar/par_all22_prot_cmap_phmd.inp
set gbradius = toppar/radius_gbsw.str
set nstage = 2
set stage = 0
set pdbname = cpept
set ph = 4.0
set salt = 0.15
set segid = prot
set temp = 298.0
set nstep = 1000
set phmdpar = toppar/phmd.in

faster on 
random uniform

! 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 in psf and coordinate files
! --------------------------------
open read unit 1 form name struct/@{pdbname}_h.psf
read psf card unit 1 
close unit 1

open read unit 1 form name struct/@{pdbname}_h.pdb
read coor pdb unit 1 
close unit 1

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

stream @gbradius
gbsw reset
gbsw sgamma 0.005 nang 50 conc @salt

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

! Turn on Shake
! -------------
shake bonh para tol 1.0e-8

! Run Langevin dynamics
! ---------------------
scalar fbeta set 5.0 sele (.not. type h*) end ! set friction coefficients

label run_dyna

open unit 23 read form name @phmdpar
open unit 25 write form name data/@{pdbname}-@{stage}.lambda
phmd par 23 wri 25 ph @ph npri 250 barr 1.75 bartau 2.5 temp @temp lam -
! if sele not present all titrable residues will be titrated
sele (resn Glu .or. resn Hsp .or. resn Lys .or. resn Arg) end 


if @stage eq 0 then
  set cmd start
  cons harm force 1.0 sele (.not. type h*) end
else
  cons harm clear
  set cmd restart
  calc stagem1 = @stage - 1
  open unit 29 read form name data/@{pdbname}-@{stagem1}.res
endif

open unit 30 write unform name data/@{pdbname}-@{stage}.dcd
open unit 31 write form name data/@{pdbname}-@{stage}.res
open unit 32 write form name data/@{pdbname}-@{stage}.eng

if @cmd .eq. start then
 dynamics leap lang @cmd timestep 0.002 iseed 314159 -
   nstep   @nstep -
   firstt  @temp finalt  @temp  tstruc @temp tbath @temp  -
   iasors  0 iasvel  1 inbfrq -1 nsavv 0 ntrfrq 1000 -
   nprint  250 iprfrq 10000 isvfrq 25000 nsavc 250 -
   iunrea -1 iuncrd 30 iunwrit 31 kunit 32
else
 dynamics leap lang @cmd timestep 0.002 -
   nstep   @nstep -
   firstt  @temp finalt  @temp  tstruc @temp tbath @temp  -
   iasors  0 iasvel  1 inbfrq -1 nsavv 0 ntrfrq 1000  -
   nprint  250 iprfrq 10000 isvfrq 25000 nsavc 250 -
   iunrea 29 iuncrd 30 iunwrit 31 kunit 32
endif

open unit 1 write format name data/@{pdbname}-@{stage}.pdb
write coor pdb unit 1

if @stage lt @nstage then
  incr stage by 1 
  goto run_dyna
endif

stop