* Calculate accessible surface
* for each residu and calculate the
* difference bare-binding.
*
bomblev -1 

! Code for this job

$SET$

! Topology, parameters and PSF
$LIB$

!Initial psf 
open unit 1 card read name @path/scr/@name_complex.psf
read psf card unit 1
close unit 1

! Initial coordinates
open unit 1 card read name @path/scr/@name@ext_cx-cha.crd.his
read coor card unit 1
close unit 1


! Result file
open unit 11 write formatted name @path/scr/@name.sas.surf
write title unit 11
*# JOBNAME: @out
*# Analytical Accessible Surface Calculation
*# Resnum  Resname  Acc. Surf.     TOT.
*
close unit 11


! let first put the total surf to zero
set tot 0
! Now calculate the Acc. Surf. of the "bare" protein
! and "bare" first min
! Radii from param file, use weight to read them in .crd
! probe radius of 1.4 A.; accu=0 -> analyt. surf.

coor surf access accu 0.0 rprob 1.4 sele @fp end
scalar wmain mult -1
scalar wmain store 1
set a = ?AREA
coor surf access accu 0.0 rprob 1.4 select @sp end
scalar wmain mult -1 
scalar wmain +store 1
scalar wmain recall 1
incr a by ?AREA


!Allows to know the first and the last atom number of the fp
define first sele @fp end
set firstatomfp  ?SELATOM
set lastatomfp  ?SELATOM
incr lastatomfp by ?NSEL
decr lastatomfp by 1

!Allows to know the first and the last atom number of the sp
define second sele @sp end
set firstatomsp ?SELATOM
set lastatomsp ?SELATOM
incr lastatomsp by ?NSEL
decr lastatomsp by 1

! Allow to know the number of residue of the first part
set res = 0

!loop begins at the first atom of the fp
set j = 0
incr j by @firstatomfp
label loop10

!allows to know the segid and the residue of the atom j and the number of atom in the residue
define atom sele bynum @j end
set residue ?SELRESI
set segid ?SELSEGI
define resi sele resid @residue .and. segid @segid end
set resiatom ?NSEL

incr res by 1
incr j by @resiatom
if j le @lastatomfp goto loop10 

!loop begins at the first atom of the sp
set j = 0
incr j by @firstatomsp
label loop11

!allows to know the segid and the residue of the atom j and the number of atom in the residue
define atom sele bynum @j end
set residue ?SELRESI
set segid ?SELSEGI
define resi sele resid @residue .and. segid @segid end
set resiatom ?NSEL

incr res by 1
incr j by @resiatom
if j le @lastatomsp goto loop11 



! Now let's loop over residues; this will calculate
! the accessible surface for all protein residues 
! and one bare minimum

set i  1
label loop1
scalar wmain stat sele ires @i end
incr i by 1
incr tot by ?STOT

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@name.sas.surf
write title unit 12
*  ?SELRESI  ?SELRESN  ?STOT  @tot
*
close unit 12
FORMAT

if i le @res goto loop1 


! Compare the 2 values of total surf, just to check
open unit 12 append formatted name @path/scr/@name.sas.surf
write title unit 12
* sum of individual = @tot   TOTAL surf = @a
*
close unit 12



set tot 0

!let's do some bookeeping
scalar wmain recall 1
scalar wmain store 2
! Now calculate the Acc. Surf.
! Radii from param file, use weight to read them in .crd
! probe radius of 1.4 A.; accu=0 -> analyt. surf.
coor surf access accu 0.0 rprob 1.4
print coor
scalar wmain  +store 2 
scalar wmain recall 2 

set a = ?AREA
set k  1
label loop3
scalar wmain stat sele ires @k end
incr k by 1
incr tot by ?STOT

! let's print all that out
FORMAT (F9.3)
open unit 12 append formatted name @path/scr/@name.sas.surf
write title unit 12
*  ?SELRESI  ?SELRESN  ?STOT  @tot
*
close unit 12
FORMAT

if k le @res goto loop3 

! Compare the 2 values of total surf, just to check
open unit 12 append formatted name @path/scr/@name.sas.surf
write title unit 12
* sum of individual = @tot   TOTAL surf = @a
*
close unit 12

stop