* Membrane System : step3.inp
* For randomly picking the lipids from the 2000 member DMPC library
*

stream system.str

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

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

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

!Read big-spheres (POL)
stream pol_toppar.str

read sequence POL @Nlipid
generate HEAD warn

open read card unit 1 name step2.crd
read coor card unit 1
close unit 1

!==============================================================================
! Randomly choosing the lipids and hydration water

!Random number generation for picking lipids
RAND UNIF  ISEED 1384716341

set LipidNum = 1

label loop_gen

label repeatrand

set a ?RAND
if a le 0.5 set Lipdir conf1   ! conf1 and conf2 are two directories 
if a gt 0.5 set Lipdir conf2   ! containing lipid structures

if a le 0.001 goto repeatrand
calc N = int( @a * 1000 )   

open read card unit 1 name dmpc/@Lipdir/dmpc_@N.crd
read sequence coor unit 1 
generate L@LipidNum setup noangle nodihedrals first none last none
close unit 1

open read card unit 1 name dmpc/@Lipdir/dmpc_@N.crd
read coor card unit 1 append
close unit 1

! check to delete any lipids that have tail near head groups
coor stat sele segid L@LipidNum .and. type C214 end
if ?ZAVE ge -10.0 then
   delete atoms sele segid L@LipidNum end 
   goto repeatrand
endif
coor stat sele segid L@LipidNum .and. type C314 end
if ?ZAVE ge -10.0 then
   delete atoms sele segid L@LipidNum end 
   goto repeatrand
endif

! put the X-Y center of mass of the lipid head at the origin 
coor stat mass select segid L@LipidNum end
coor trans xdir -?XAVE ydir -?YAVE zdir 0.0 select segid L@LipidNum end

! to put the P-N of the polar head at the z=0 point
coor stat sele segid L@LipidNum .and. ( type P .or. type N ) end
coor trans zdir -?ZAVE select segid L@LipidNum end 

! flip the orientation of the lipids in the lower monolayer if head Z le 0
coor stat sele segid HEAD .and. resid @LipidNum end 
if ?ZMIN le 0  coor rotate xdir 1.0  phi 180.0 select segid L@LipidNum end

! put the head group at the location of the spheres
coor stat sele segid HEAD .and. resid @LipidNum end
!wp:coor trans xdir ?XMIN ydir ?YMIN zdir ?ZMIN select segid L@LipidNum end
if ?ZMIN le 0.0 -
   coor trans xdir ?XMIN ydir ?YMIN zdir @ZPosBOT select segid L@LipidNum end
if ?ZMIN gt 0.0 -
   coor trans xdir ?XMIN ydir ?YMIN zdir @ZPosTOP select segid L@LipidNum end

! check whether the lipid penetrates the peptide
coor dist cut 1.5 sele segid L@LipidNum .and. (type c2* .or. type c3*) end  -
                  sele .not. ( segid L* .or. type POLO ) end

set b  ?npair
if b ne 0 delete atoms sele segid L@LipidNum end 
if b ne 0 goto repeatrand

! continue for total number of lipids that are needed
incr LipidNum by 1
if LipidNum le @NLipid goto loop_gen
!==============================================================================

! delete the vdW spheres 
delete atom sele type POLO end

ic param
ic build
hbuild

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

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

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

stop