* Membrane System : step4.inp
* Add BULK water boxes
*

stream system.str

!Read topologie
open read card unit 10 name @charmmdir/top_all27_prot_lipid.rtf
read  rtf card unit 10

!Read parameters
open read card unit 20 name @charmmdir/par_all27_prot_lipid.prm
read para card unit 20

!==============================================================================
! The height should be changed based up the size of the protein

calc BoxSizeX  = @BoxsizeXY
calc BoxSizeY  = @BoxsizeXY
calc BoxSizeZ  = 22

stream waterbox.str
delete atom sele resname TIP3 end
!==============================================================================

!Read PSF
open read card unit 10 name step3.psf
read psf  card unit 10

open read card unit 10 name step3.crd
read coor card unit 10 
close unit 10

!Read Water box
open read card unit 10 name waterbox.crd
read sequence coor unit 10
generate WAT1 setup noangle nodihedrals

open read card unit 10 name waterbox.crd
read coor card unit 10 append

!Duplicate the water box
generate WAT2 setup duplicate WAT1
coor duplicate select segid WAT1 end select segid WAT2 end

!Translate one box to +Z
coor stat sele segid WAT1 .and. type OH2 end
coor trans zdir -?zmin sele segid WAT1 end
coor trans zdir   12.0 sele segid WAT1 end

!Translate one box to -Z
coor stat sele segid WAT2 .and. type OH2 end
coor trans zdir -?zmax sele segid WAT2 end
coor trans zdir  -12.0 sele segid WAT2 end

!Remove all the water overlapped with protein + lipids
define WBOX sele segid WAT1 .or. segid WAT2 end
define TOTO sele .not. WBOX .and. .not. hydrogen end
delete atom sele .byres. ( WBOX .and. ( TOTO .around. 2.6 ) ) end

open write card unit 1 name step4.psf
write  psf card unit 1

open write card unit 1 name step4.pdb
write coor pdb unit 1

open write card unit 1 name step4.crd
write coor card unit 1

calc BoxsizeZ = ( 12.0 + @BoxsizeZ ) * 2.0

open write card  unit 90 name system.str
write title unit 90
* set  charmmdir = @charmmdir
* calc Nlipid    = @Nlipid
* calc BoxsizeXY = @BoxsizeXY
* calc BoxsizeZ  = @BoxsizeZ
* calc ZPosTop   =  17      ! Position of the top layer (for DPPC use 18)
* calc ZPosBot   = -17      ! Position of the bottom layer (for DPPC use -18)
*

!========================================================================
! Separate Lipid and Water

delete atom sele .not. resname TIP3 end

open write card unit 1 name step4_tip3.crd
write coor card unit 1

delete atom sele all end
open read card unit 1 name step4.psf
read  psf card unit 1

open read card unit 1 name step4.crd
read coor card unit 1

delete atom sele .not. resname DMPC end

open write card unit 1 name step4_dmpc.crd
write coor card unit 1

stop