* read pdb 
* add protons and water molecules
*

bomblev -2

! { parameters passed }
set diam_w  2.8  ! { pseudo diameter of water model }
$SET$

! Read in Topology and  Parameter files
open unit 1 card read name @path/lib/top_all22_prot_na.inp
read RTF card unit 1
close unit 1
open unit 1 card read name @path/lib/par_all22_prot_na.inp
read PARA card unit 1
close unit 1

!-----------------------------------------------> define protein
! Initial psf 
open unit 1 card read name @path/scr/@pdb.psf.his
read psf card unit 1
close unit 1

! Initial coordinates
open unit 1 card read name @path/scr/@pdb.mini_@ext_abnr.crd
read coor card unit 1
close unit 1

! BEGIN { retranslate system to reaction region origin }
define origin select all end

! { find the centre of the reaction region }
coordinate statistics select ( origin ) end

! { make this point the origin }
coordinate orient sele all end
!coordinate translate xdirection -?XAVE ydirection -?YAVE zdirection -?ZAVE

! { verify }
coordinate statistics select ( origin ) end

! end { translate system to reaction region origin }

! { read box of equilibrated water }
open unit 1 read formatted name @path/lib/box90r.crd
read sequence coor unit 1 
rewind unit 1

! { generate new segment }
generate bulk noangle nodihedral 

! { get water box }
read coor card unit 1 append 
close unit 1

! BEGIN { retranslate water to origin }

define origin select  segid bulk end

! { find the centre of the reaction region }
coordinate statistics select ( origin ) end

! { make this point the origin }
coordinate rotate  ydir 1.0 phi -90. sele origin end
coordinate statistics select ( origin ) end
coordinate translate xdirection -?XAVE ydirection -?YAVE zdirection -?ZAVE -
    sele segid bulk end

! { verify }
coordinate statistics select ( origin ) end

! end { translate system to reaction region origin }
   
! verify location of protein and water sphere

COOR STATISTICS SELE ( SEGID bulk .or. SEGID W ) END
COOR STATISTICS SELE ( SEGID @fc .or SEGID @sc ) END

! { delete waters which overlap with protein or ca ion}
delete atom sort -
            select .byres. (segid bulk .AND. type oh2 .and. -
            ((.not. (segid bulk .OR. hydrogen)) .around. @{diam_w})) end 

!find out number of waters added
define count select (segid bulk .AND. type oh2) end 
set c ?nsel

! write out the psf and the coordinates
open unit 1 write formatted name @path/scr/@pdb.dyna_solv.psf 
write psf card unit 1
*  @job
*  PSF 
*  INITIAL SOLVENT OVERLAY: @c waters added to @pdb.
*

open unit 1 write formatted name @path/scr/@pdb.dyna_solv.crd 
write coordinate card unit 1
*  @job
*  COORdinates 
*  INITIAL SOLVENT OVERLAY: @c waters added to @pdb.
*

OPEN WRITE CARD NAME @path/scr/@pdb_solv.pdb UNIT 20
WRITE COOR PDB UNIT 20
*  @job
*  COORdinates 
*  INITIAL SOLVENT OVERLAY: @c waters added to @pdb.
*

energy

STOP