* Membrane System : step1.inp
* Caculate the cross-sectional area of a molecule
*

set charmmdir = .

!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

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

!Read Coordinate
open read card unit 10 name pept.pdb
read coor pdb  unit 10 resid

!===== Cross Section ======

open write card  unit 51 name step1.plo
write title unit 51
* #Cross sections along z axis
*

calc dcel   = 0.2  ! this is an approximate one!
calc factor = 1.0 / @dcel

coor stat sele all end
calc x1 = ?xmin - 5
calc x2 = ?xmax + 5
calc y1 = ?ymin - 5
calc y2 = ?ymax + 5
calc xncel = int( ( @x2 - @x1 ) / @dcel )
calc yncel = int( ( @y2 - @y1 ) / @dcel )

calc a = int( ?zmin ) 
calc b = @a + @dcel

calc toparea = -1.0
calc botarea = -1.0

label DOAREA

scalar wmain = radius
scalar wmain mult 2.5   sele resname TIP3 end
scalar wmain add  1.529 sele .not. resname TIP3 end
scalar wmain mult 0.85  sele .not. resname TIP3 end

coor search -
     xmin @x1  xmax @x2  xgrip @xncel -
     ymin @y1  ymax @y2  ygrip @yncel -
     zmin @a   zmax @b   zgrip 5

calc c = ?volume * @factor

if a .lt. 0 if c .gt. @botarea set botarea = @c
if a .gt. 0 if c .gt. @toparea set toparea = @c

write title unit 51
* @b  @c
*

incr a by @dcel
incr b by @dcel
if b le ?ZMAX goto DOAREA

open write card  unit 51 name step1.str
write title unit 51
* calc peptarea = ( @toparea + @botarea ) / 2.0
*

stop