*  This is a general template script to equilibrate a solvated 
*  simulation "box" using periodic boundaries with a few specific volumes.
*  The script will read in a psf file for a solute system already generated
*  in a:
*  1) cubic volume w/ crystal
*  2) truncated octahedral box w/ crystal (reoriented from images)
*  3) carry-out an equilibration npt dynamics run w/ PME
***
*  Required files: top_all22_prot.inp, par_all22_prot.inp
**                 rhomba.xtl, output from protein-generate_solvate.inp
*

!**_______________________________________________________________________**
!**                    WRITTEN BY C.L. BROOKS, III                        **
!**                           OCTOBER, 1999                               **
!**                                                                       **
!**     COPYRIGHT C.L. BROOKS, III, THE SCRIPPS RESEARCH INSTITUTE        **
!**                                                                       **
!**     Modified M.F.Crowley May 2000                                     **
!**                                                                       **
!**_______________________________________________________________________**
!********
!*  The usage requires that several variables either be set in the
!*  input script or on the CHARMM command line at execution.
!*  The needed variables are:
!**  Variable   Default               Function
!*__________________________________________________________________________
!*   toppar   ($toppar)  ==> path to parameter and topology files
!**
!*   pdbfile     (...)   ==> name of charmm (crd) coordinate file containing 
!**                          coordinates of the solvated system
!**                          (pdbfile_solv-cube(octa).crd)
!**                          and the pre-generated psf file (pdbfile.psf))
!*   shape      (cubic)  ==> variable specifying whether truncated octahedral or
!**                          cubic volume used (octahedral/cubic)
!*   crystal    (false)  ==> logical variable true/false specifying whether
!**                          to use crystal orientation of simulation volume
!*   cutoff      (9)    ==> value of cutoff to be used in establishing how
!**                          far images are to be from volume "edges", often set
!**                          at or near nonbonded list cutoff value
!*   T           (298)   ==> Temperature for run.  Default 298K.
!**                          at or near nonbonded list cutoff value
!*   begin       (1)     ==> Window to begin this production run from
!**                 
!*   end         (10)    ==> Window to end this production run at
!**                 
!*   current     (begin) ==> Window last run
!**                 
!

if @?pdbfile eq 0 goto iamdead

stream @pdbfile_dimens.stream

!*********************************************************
!     get parameters set for this run
!*********************************************************
!
   goto getparams
   label doneparams

!*********************************************************
!     get topology and force field for this run
!*********************************************************
!
!prnlev 0 node 0
wrnlev 0
open unit 1 read form name @toppar/top_all22_prot.inp
read rtf card unit 1
close unit 1

open unit 1 read form name @toppar/par_all22_prot.inp
read param card unit 1
close unit 1

prnlev 5 node 0
wrnlev 2

!*********************************************************
!     get psf for this system
!         This psf does not have the solvent
!*********************************************************
!
  open unit 1 read form name @pdbfile.psf
  read psf card unit 1
  close unit 1

!*********************************************************
!           Add water
!*********************************************************
  read sequ tip3 @nwat
  generate wat setup noangl nodihe

goto @task
!*********************************************************
!                  EQUILIBRATION
!*********************************************************
!
label equil
!*********************************************************
!         set up periodic boundary cond.
!*********************************************************
!
 if shape eq cubic open unit 1 read form name @pdbfile_solv-cube.crd
 if shape eq octahedral open unit 1 read form name @pdbfile_solv-octa.crd
 read coor card unit 1
 close unit 1

 if shape ne  cubic goto doocta
  Crystal define cubic @7 @8 @9 @theta @theta @theta
  crystal build cutoff @cutoff noper 0
  crystal write card unit 6
 label doocta
 if shape ne octahedral goto donextl
  Crystal define octa @7 @8 @9 @theta @theta @theta    
  open read card unit 1 name "rhomba.xtl"
  crystal read card unit 1
  close unit 1
 label donextl

image byseg xcen 0.0 ycen 0.0 zcen 0.0 select .not. resname tip3 end
image byres xcen 0.0 ycen 0.0 zcen 0.0 select resname tip3 end

!*********************************************************
!  turn on faster options and set-up SHAKE
!*********************************************************
 faster on

!*********************************************************
!   For a the first pass at equilibration,
!       do a minimization with solute fixed
!       and another with solute harmonically restrained
!*********************************************************
!
 if status ne start goto dodyn
  cons fix  select .not. resn tip3   end

 minimize conj nstep 30 tolenr 1.00 nprint 10 -
         cutnb 10.0 rdie eps 2.0 vswitch shift bycb
 cons fix  select none   end
 cons harm force 30.0  select .not. resn tip3   end
 minimize conj nstep 30 tolenr 1.00 nprint 10 -
         cutnb 10.0 rdie eps 2.0 vswitch shift bycb
 cons harm clear


  set begin 1
  set end 50
  set current @begin

!  Set-up directories for dynamics output if system allows. 
!            (broken for some parallel machines)
! system "mkdir kun res crd coor size"
!
!  Harmonically restrain the backbone to permit the solvent to become
!  "equilibrated" to the solute. First pass in equilibration only,
!   then turn off restraints.
!
cons harm force 5.0 mass select -
           ( type c .or. type o .or. type ca .or. type n ) end

!*********************************************************
!                DYNAMICS
!*********************************************************
label dodyn
 shake bonh tol 1.0e-8 para

title
*  Initial dynamics for @pdbfile.
*  Dynamics initiated from cooled structure and heated to 298K over
*  50K steps (100ps) of dynamics.  Production for @current 20ps window.
*  Harmonic force restraints applied during the first 20ps using a force
*  constant of 5*mass.
*

if current eq @begin goto skip
goto skip

decr current by 1
open unit 4 read form name "res/@PDBFILE_e@CURRENT.res"
incr current by 1

label skip
open unit 3 write form name "res/@PDBFILE_e@CURRENT.res"
open unit 20 write form name "kun/@PDBFILE_e@CURRENT.kun"
open unit 21 write unform name "crd/@PDBFILE_e@CURRENT.crd"

prnlev 4 node 0
! Run dynamics in periodic box
dynamics cpt leap @status timestep 0.002 nstep 10000 nprint 50 iprfrq 1000 -
     firstt 298 finalt 298 twindl -5.0 twindh 5.0 -
     ichecw 1 teminc 30.0 ihtfrq 50 ieqfrq 200 -
     iasors 1 iasvel 1 iscvel 0 isvfrq 1000 -
     iunwri 3 nsavc 100 nsavv 0 iunvel 0 -
     iunread 4 iuncrd 21 kunit 20   - !{* Nonbond options *}
     inbfrq -1 imgfrq -1 ilbfrq 0 bycb - !use bycube image nonbond list generator
     eps 1.0 cutnb @cutoff cutim @cutoff ctofnb @ctofnb ctonnb @ctonnb vswi -
     Ewald kappa 0.320 pmEwald order 4 fftx 64 ffty 64 fftz 64 ntrfq 200 - !PME
     pconstant pmass 100 pref 1 pgamma 10 !  Constant pressure


open unit 3 write form name "size/@PDBFILE_e@CURRENT.size"
write title unit 3
** Box size stream file
**
* set xdim ?xtla
* set ydim ?xtlb
* set zdim ?xtlc
* set alpha ?xtlalpha
* set beta ?xtlbeta
* set gamma ?xtlgamma
* end
*
prnlev 5 node 0

open unit 1 write form name -
    "coor/@PDBFILE_e@CURRENT.chr"
write coor card unit 1
*  Coordinates for @pdbfile
*  Equilibration after @current iterations of 10K steps.
*  Current boxsize is ?xtla
*

if current ne @begin goto increment
set status restart
cons harm force 0.0 select -
           ( type c .or. type o .or. type ca .or. type n ) end

label increment
incr current by 1
if current le @end goto dodyn

stop

label format
!*********************************************************
!                FORMAT
!*********************************************************
!
calc nunit = @end - @begin +1
  open unit 30 unform read name @pdbfile_e@begin.crd
  traj iread 30
  traj read
  close unit 30
set start ?start
set skip  ?skip
set nframes ?nfile
calc nfile = ?nfile * @nunit
calc stop =  @skip * @nframes * @nunit + @start - @skip

set unit 30
set current @begin
label dynopen
  open unit @unit unform read name @pdbfile_e@current.crd
  incr unit
  incr current
if @current le @end goto dynopen 
open unit 20 write unform name @pdbfile_e@begin-@end.crd
merge firstu 30 nunit @nunit skip @skip output 20 nfile @nfile
close unit 20
close unit 30

open unit 30 read unform name @pdbfile_e@begin-@end.crd
open unit 20 write form name @pdbfile_e@begin-@end.fcrd

dyna format first 30 nunit 1 output 20 -
            skip @skip begin @start stop @stop

stop
!
label unformat
!*********************************************************
!                  UNFORMAT
!*********************************************************

set unit 30
  open unit @unit unform write name @pdbfile_e@begin-@end.crd
!  open unit @unit unform write name a12

open unit 20 read form name @pdbfile_e@begin-@end.fcrd
dyna unformat input 20  output 30
stop
label unmerge
!*********************************************************
!                unmerge
!*********************************************************
!

calc nunit = @end - @begin +1
  open unit 30 unform read name @pdbfile_e@begin-@end.crd
  traj iread 30
  traj read
  close unit 30
set start ?start
set skip  ?skip
calc nframes ?nfile 
calc stop =  @skip * @nframes  + @start - @skip
calc nfile = ?nfile / @nunit

set unit 30
set current @begin
label mergopen
  open unit @unit unform write name @pdbfile_e@current.crd
  incr unit
  incr current
if @current le @end goto mergopen 
open unit 20 unform read name @pdbfile_e@begin-@end.crd

merge firstu 20 nunit 1 output 30 nfile @nfile -
            skip @skip begin @start stop @stop

stop
!===========================================================
!*********************************************************
!*******    END OF SCRIPT  subroutines follow     ********
!*********************************************************
!===========================================================

label iamdead
!system "banner Must specify coordinate/psf filename"

stop


!*********************************************************
!     get parameters set for this run
!

label getparams
if @?toppar eq 0 set toppar = "$toppar"
if @?T eq 0 set T = 298
if @?begin eq 0 set begin = 1
if @?end eq 0 set end = 10
if @?current eq 0 set current = 1
if @?shape eq 0 set shape = cubic
if @?crystal eq 0 set crystal = false
if @?cutoff eq 0 set cutoff = 9
if @?status eq 0 set status start

! Dimension of the solvent volume
Calc 7 = @size * 1.005 ! Initially set boxsize slightly larger for CPT equilibration
  if @shape ne octahedral goto sizeok
   if @crystal eq true calc 7 = @7 * sqrt(3) / 2.
  label sizeok
set 8 @7
set 9 @8

!  Set-up cutoffs for non-bonded list and energy calculations
Calc cutoff = @cutoff + 1  ! Add additional 1A to cutoff from build
set ctofnb = @cutoff
Calc ctofnb = @ctofnb - 1.5
set ctonnb = @ctofnb
Calc ctonnb = @ctonnb - 1.0

goto doneparams