{+ file: buried_surface.inp +}
{+ directory: general +}
{+ description: Lee & Richards buried surface accessibility calculation +}
{+ comment:
           Coordinate files will be written with:
             - interface residues only
             - all atoms, interface residues will have a B-factor=1
               all other residues will have B-factor=0). +}
{+ authors: Axel T. Brunger, and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: B. Lee, and F.M. Richards, The Interpretation of
              Protein Structures:  Estimation of Static Accessibility,
              J. Mol. Biol. 55, 379-400, 1971. +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) are not quoted
   - do not remove any evaluate statements from the file -}

{- begin block parameter definition -} define(

{======================= molecular structure =========================}

{* structure file *}
{===>} structure_infile="eg1_dimer.mtf";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param";
{===>} parameter_infile_2="CNS_TOPPAR:carbohydrate.param";
{===>} parameter_infile_3="";
{===>} parameter_infile_4="";
{===>} parameter_infile_5="";

{* coordinate file *}
{===>} coordinate_infile="eg1_dimer.pdb";

{========================== parameters ==============================}

{* probe radius in Angstroms *}
{===>} probe=1.4;

{========================== atom selection ===========================}

{* select atoms of first molecule *}
{===>} atom_select_1=(segid AAAA or segid BBBB);

{* select atoms of second molecule *}
{===>} atom_select_2=(segid CCCC or segid DDDD);

{=========================== output files ============================}

{* root name for output files *}
{+ list:
        buried area listing: <output_root>.list
        coordinates (B=1 for interface residues): <output_root>.pdb
        coordinates (interface residues only): <output_root>_interface.pdb +}
{===>} output_root="buried_surface";

{===========================================================================}
{        things below this line do not normally need to be changed          }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates @&coordinate_infile

 parameter
   if ( &BLANK%parameter_infile_1 = false ) then
      @@&parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      @@&parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      @@&parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      @@&parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      @@&parameter_infile_5
   end if
 end

 evaluate ($list_outfile=&output_root + ".list")

 set display=$list_outfile end

 do ( store1=0 ) ( all )

 surface
    mode=access
    rh2o=&probe
    selection=&atom_select_1
    radius=vdw
 end

 show sum ( rmsd ) ( &atom_select_1 )
 evaluate ($molecule_1=$result)

 do ( store1=rmsd ) ( &atom_select_1 )

 surface
    mode=access
    rh2o=&probe
    selection=&atom_select_2
    radius=vdw
 end

 show sum ( rmsd ) ( &atom_select_2 )
 evaluate ($molecule_2=$result)

 do ( store1=rmsd ) ( &atom_select_2 )

 surface
    mode=access
    rh2o=&probe
    selection=( &atom_select_1 or &atom_select_2 )
    radius=vdw
 end

 show sum ( rmsd ) ( &atom_select_1 or &atom_select_2 )
 evaluate ($molecule_both=$result)

 do (store2=rmsd) ( &atom_select_1 or &atom_select_2 )

 evaluate ($buried=$molecule_1 + $molecule_2 - $molecule_both)

 do ( b=0 ) ( all )

 do ( b=max(0,step(store1-store2-4) )) ( &atom_select_1 or &atom_select_2 )

 do ( b=1 ) ( ( &atom_select_1 or &atom_select_2 ) and
                byresidue ( attribute b > 0) )

 display >>>> Buried surface area calculation:
 display >>>>   structure file= &STRIP%structure_infile
 display >>>>   coordinate set= &STRIP%coordinate_infile
 display >>>>   molecule 1: &atom_select_1
 display >>>>   molecule 2: &atom_select_2
 display
 display  surface area for molecule 1: $molecule_1  A^2
 display  surface area for molecule 2: $molecule_2  A^2
 display  surface area for the complex of both molecules: $molecule_both  A^2
 display
 display  surface area buried at interface: $buried  A^2

 evaluate ($coord_outfile=&output_root + ".pdb")

 write coordinates 
   format=pdbo
   output=$coord_outfile
 end

 evaluate ($coord_outfile=&output_root + "_interface.pdb")

 write coordinates 
   format=pdbo 
   output=$coord_outfile 
   selection=( attribute b > 0 )
 end

 stop