* generate DOPC lipid bilayer with water
*

!generation via modular approach where the head group (PCGL)
!is joined to two aliphatic tails (OLEO) via the patches
!EST1 and EST2

bomlev -1

faster on

set id1 dopc
set id2 1

! 1: topology
! 2: parameter
! 3: starting coordinates
! 4: cutnb
! 5: ctonnb
! 6: ctofnb 
! 7: electrostatic cutoff regime (shift or switch)
! 8: electrostatic pairing method (atom or group)
! 9: VDW pairing method (vatom or vgroup)
! VDW cutoff regime is vwsitch
set 1 top_all27_lipid.rtf
set 2 par_all27_lipid.prm
set 3 12.0  ! cutimg
set 4 12.0  ! cutnb
set 5  8.0  ! ctonnb
set 6 10.0  ! ctofnb
set 7 switch
set 8 atom
set 9 vatom

open unit 10 form read name @1
read rtf card unit 10
close unit 10

open unit 10 form read name @2
read para card unit 10
close unit 10

set lname L00
set count 1

label gen_loop

if count ge 10 set lname L0

read sequence card
* DPPC via modular approach
* 
3
OLEO PCGL OLEO

GENERATE @lname@count FIRST none LAST none SETUP WARN

patch est1 @lname@count 2 @lname@count 1
patch est2 @lname@count 2 @lname@count 3

incr count by 1
if count le 72 goto gen_loop

autogenerate angles dihedrals

!generate water
read sequence tip3 386

generate wat first none last none noangle nodihedral

!read in coordinates from previous simulation
open unit 20 read form name dopc_1500ps.pdb
read coor pdb unit 20

!coor print

coor stat
coor stat sele .not. hydrogen end
coor stat sele (resn palm .or. resn pcgl) .and. .not. hydrogen end
coor stat sele resn tip3 .and. .not. hydrogen end

!write psf and coordinates
open unit 20 write form name @id1_@id2.psf
write psf card unit 20
* 72 DPPC + 2094 WATER MOLECULES
* XDIM = 47.5857  YDIM = 47.5857  ZDIM = 66.8901
*

open unit 20 write form name @id1_@id2.crd
write coor card unit 20
* 72 DOPC + 386 WATER MOLECULES
* XDIM = 47.5857  YDIM = 47.5857  ZDIM = 66.8901
*

!create images and calculate energy

! cell parameters
set a 46.20000
set b 46.20000
set c 49.14077

crystal define tetrahedral @a @b @c 90. 90. 90.

crystal build noper 0 cutoff @3

!turn on image centering of solvent
image byres xcen 0.0 ycen 0.0 zcen 0.0 -
     sele resn tip3 end

!kappa should be optimized based on the value of
!kappa at which the 1st derivative of the energy
!with respect to kappa goes to 0.0
set kap 0.36  !NOT Optimized
set ord 6

!fftx dimensions should correspond to 1 A or less; must
!be power of 2, 3 or 5
energy inbfrq -1 imgfrq -1 ihbfrq 0 -
ewald pmewald kappa @kap fftx 64 ffty 64 fftz 64 order @ord -
@7 @8 @9 vfswitch cutimg @3 cutnb @4 ctofnb @6 ctonnb @5

stop