*  Run and analyze dynamics of the protein G domain
*  swapped dimer in TIP4P water using the fluctuating charge
*  force field.
*

!***REQUIRED INPUT VARIBLES***!
!
if @?nwat eq 0     set nwat = 2623         !number of water molecules
if @?watseg eq 0   set watseg = w00w       !segid for water
if @?protbase eq 0 set protbase = pro      !base segment name for protein
if @?length eq 0   set length = 58.413132  !dimension of solvation box

if @?seed ne 0 goto gotseed
  system "date +%H%M%S | awk '{seed=$0*2+1;print "* Title"; print "*"; print "set seed = "seed}' > seed.stream"
  stream seed.stream
  system "rm seed.stream"
label gotseed

! Read in the default CHEQ topology and parameter files
open unit 1 read form name top_allcheq06_prot_cmap.rtf
read rtf card unit 1
close unit 1

open unit 1 read form name par_allcheq06_prot_cmap.prm
read param card unit 1
close unit 1

! We use the MMTSB Tool Set to assist in a number of these
! steps. First we generate the separate chains using the commands:
!convpdb.pl -center -model 1 -out charmm22 -segnames 1Q10.pdb > 1q10_a+b.pdb
!convpdb.pl -center -chain A -out charmm22 1q10_a+b.pdb > 1q10_a.pdb
!convpdb.pl -center -chain B -out charmm22 1q10_a+b.pdb > 1q10_b.pdb
!
!
! Then we solvate the system and minimize it initially using 
! a nonpolarizable force field. See the tutorial example at
! http://mmtsb.scripps.edu/workshops/mmtsb-ctbp_2006/Tutorials 
! for details.
!

open unit 1 read form name 1q10_a.pdb
read sequ pdb unit 1
generate gb1a first ace last ct3 setup
close unit 1

open unit 1 read form name 1q10_b.pdb
read sequ pdb unit 1
generate gb1b first ace last ct3 setup
close unit 1

! The system was solvated with @nwat water molecules
read sequ TIP4 @nwat !Solvate with TIP4FQ water model
generate wat  first none last none setup noangl nodihe

! set up tip4p bisector
lonepair bisector dist 0.15 angle 0.0 dihe 0.0 -
       sele atom wat * OM end -
       sele atom wat * OH2 end -
       sele atom wat * H1 end -
       sele atom wat * H2 end

rename segid @{protbase}a select segid gb1a end
rename segid @{protbase}b select segid gb1b end
rename segid @watseg select segid wat end
rename resname tip3 select segid @watseg end
!minimized solvated coordiantes
open unit 1 read form name 1q10_a+b_solvmin.pdb  
read coor pdb unit 1 resid
close unit 1
rename segid gb1a select segid @{protbase}a end
rename segid gb1b select segid @{protbase}b end
rename segid wat select segid @watseg end
rename resname tip4 select segid wat end

! Add OM site using coor shake
coor shake

coor copy comp

! set the usual fluctuating charge parameters
!  Masses are from Rick, Stuart, and Berne paper.
cheq norm byres sele segid wat end
cheq norm byres sele .not. segid wat end
cheq tip4 sele segid wat end
cheq flex select .not. segid wat end
cheq qmas cgma 0.000069 tsta 0.01 sele segid wat end
cheq qmas cgma 0.000069 tsta 0.01 sele .not. segid wat end

set angle =  109.4712206344907

calc cutoff = @length / 2
let cutoff = mini  @cutoff   12 

set cutim = @cutoff
Calc ctonnb = @cutoff - 3
Calc ctofnb = @cutoff - 2

!  Set-up crystal parameters
crystal defi octa @length @length @length -
    109.4712206344907 109.4712206344907 109.4712206344907
crystal build cutoff @cutoff Noper 0
imag byres xcen 0 ycen 0 zcen 0 select segid wat end
imag byseg xcen 0 ycen 0 zcen 0 select .not. segid wat end

! Set-up required non-bonded options
faster 1
update -
bycb inbfrq -1 imgfrq -1 nbscale 1.5 -
e14fac 0 nbxmod -5 -  ! These are required for CHEQ dynaamics
cdiel eps 1.0 cutnb @cutoff cutim @cutim ctofnb @ctofnb ctonnb @ctonnb vswi -  
Ewald kappa 0.320 pmEwald order 4 fftx 64 ffty 64 fftz 64 !PM

! set this up so we could branch to analysis of we so choose
if @?analysis ne 0 goto analysis
label dodynam

shake bonh tol 1.0e-8 param mxit 2000

! minimize configuration w/o fluctuating charge
mini sd nstep 100 nprint 10

! With the scale commands below, we can calculate the shift in
! the charges from initial values to values that are "equilibrated"
! to the new environment.

! Turn cheq on and minimize again
cheq on

scalar charge store 1  ! charges before equilibraiton

! First just the charges
mini sd nstep 100 nprint 10 -
cheq cgmd 1 noco

! Now the whole configuration
mini sd nstep 500 nprint 10 -
cheq cgmd 1

scalar charge store 2

! Loop over protein atoms and write data in a file to plot
open unit 1 write form name 1q10_a+b_min_0.dat
set return = frommin
set unit = 1   
goto getdipole
label frommin

! Put a copy of the charges into the wmain array
! and we can use these to color the atoms.
scalar wmain = charge
open unit 1 write form name 1q10_a+b_fqmin.pdb
write coor pdb unit 1

open unit 1 write form name 1q10_a+b_fqmin_no-om.pdb
write coor pdb unit 1 select .not. type om end

! set up the Nose-Hoover baths for the charges
! the calling procedure is identical to the NOSE 
! call for applying the bath to nuclear degrees of freedom
!
!  the "NOSE" is replaced by "FQBA"
!!!
! Here we use one bath for the solvent and one for the peptide
!!!
FQBA 3                           ! set up '3' baths
CALL 1 sele segid wat end        ! select the atoms that will be coupled 
                                 ! to bath 1
COEF 1  QREF  0.005  TREF 1.0    ! bath parameters (usual meaning)
CALL 2 sele segid gb1a end       ! select the atoms that will be coupled 
                                 ! to bath 2
COEF 2  QREF  0.005  TREF 1.0    ! bath parameters (usual meaning)
CALL 3 sele segid gb1b end       ! select the atoms that will be coupled 
                                 ! to bath 3
COEF 3  QREF  0.005  TREF 1.0    ! bath parameters (usual meaning)
END

open unit 10 write form name 1q10_a+b_fq_0.res
open unit 11 write unform name 1q10_a+b_fq_0.dcd
open unit 12 write form name 1q10_a+b_fq_0.kun

! Set-up the dynamcis run
dynamics cpt leap start timestep 0.00050 nstep 1000 nprint 100 iprfrq 100 -
   cheq cgmd 1 cgeq 1 fqint 1 -   ! Turn on CHEQ stuff for dynamics
   firstt 298.0 finalt 298.0  twindl -5.0 twindh 5.0 iseed @seed -
   ichecw 1 teminc 0.0 ihtfrq 0 ieqfrq 0 ilbfrq 0 -
   iasors 1 iasvel 1 iscvel 0 isvfrq 10 -
   iunwri 10 nsavc 100 nsavv 0 iunvel -1 -
   iunread -1 iuncrd 11 kunit 12 - !{* Nonbond options *}
   ntrfrq 2000 echeck 100000 -
   pconstant pmass 100.0 pref 1.0 pgamma 20.0 - !  Constant pressure
   Hoover tmass 50.0 refT 298.0


scalar wmain = charge
open unit 1 write form name 1q10_a+b_fq_dyn.pdb
write coor pdb unit 1

open unit 1 write form name 1q10_a+b_fq_dyn_no-om.pdb
write coor pdb unit 1 select .not. type om end

! Loop over protein atoms and write data in a file to plot
open unit 1 write form name 1q10_a+b_dyn_0.dat
set return = fromdyn
set unit = 1
goto getdipole
label fromdyn

stop

label analysis

! open the relevant files for simulation output
open unit 10 read unform name 1q10_a+b_fq_0.dcd
open unit 11 write form name 1q10_a+b_fq_0.ene

traj iread 10

traj query unit 10
set cnt =  1
set last = ?nfile
set skip = ?skip

! Use the access to system to create the desired drectory
system "mkdir pdb"
label readtrj

    traj read

    open unit 1 write form name pdb/1q10_a+b_fq_@cnt.pdb
    write coor pdb unit 1 select .not. type om end

    energy cheq cgmd 1

    calc step = @cnt * @skip
    write title unit 11
    * @step ?ENER
    *

    incr cnt by 1
if cnt le @last goto readtrj

close unit 10
close unit 11
 
stop

!###################Subroutine getdipole####################!
!!!!
!!!!
!!!!   This routine assumes the charges have been stored in
!!!!    scalar arrays 1 (old charges) and 2 (new charges)
!!!!   and the the dipole data will be written to unit.
!!!!   The routine returns to label return.
!!!!   The dipole moment for the backbone and side chains are
!!!!   computed for each residue and written to file unit.
!!!!
label getdipole

define co select type c .or. type o end
define nh select type ca .or. type ha .or. - 
       type n .or. type hn end
define sidechain select .not. ( co .or. nh .or. segid wat ) end

set icnt = 1
label dodipole

      coor dipole select co .and. ires @icnt end
      set co = ?rdip
      scalar charge recall 1
      coor dipole select co .and. ires @icnt end
      set coo = ?rdip
      scalar charge recall 2

      coor dipole select nh .and. ires @icnt end
      set nh = ?rdip
      scalar charge recall 1
      coor dipole select nh .and. ires @icnt end
      set nho = ?rdip
      scalar charge recall 2

      coor dipole select sidechain .and. ires @icnt end
      set sc = ?rdip
      scalar charge recall 1
      coor dipole select sidechain .and. ires @icnt end
      set sco = ?rdip
      scalar charge recall 2
      write title unit 1
      * @icnt @co @coo @nh @nho @sc @sco
      *

      incr icnt by 1
if icnt le 112 goto dodipole
scalar charge recall 2

goto @return

stop