{+ file: surface_plot.inp +}
{+ directory: general +}
{+ description: Lee & Richards surface accessibility calculation +}
{+ 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="amy.mtf";

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

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

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

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

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

{* select atoms to be included in calculation *}
{===>} atom_select=(known and not(hydrogen));

{* select main chain atoms *}
{* note: atoms outside this selection will be considered to be
         side chain atoms *}
{===>} atom_main=(name ca or name n or name c or name o or name ot+);

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

{* list average or total surface area per residue *}
{+ choice: "total" "average" +}
{===>} list_mode="average";

{* output list file *}
{===>} list_outfile="surface_plot.list";

{* output coordinate file *}
{* the B-factor column will contain the per atom surface accessibility *}
{===>} coordinate_outfile="surface_plot.pdb";

{===========================================================================}
{        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
 
 set display=&list_outfile end

 surface
    mode=access
    rh2o=&probe
    selection=&atom_select
    radius=vdw
 end
 
 display >>>> Surface plot for:
 display >>>>   structure file= &STRIP%structure_infile
 display >>>>   coordinate set= &STRIP%coordinate_infile
 display
 if ( &list_mode = "average" ) then
   display                                   average accessible area
 elseif ( &list_mode = "total" ) then
   display                                    total accessible area
 end if
 display  segid  resid  resname     residue      mainchain       sidechain
 display

 set display=OUTPUT end

 for $id in id ( &atom_select and tag ) loop bfac

   show (segid) (id $id)
   evaluate ($segid=$result)
   show (resid) (id $id)
   evaluate ($resid=$result)
   show (resname) (id $id)
   evaluate ($resname=$result)

   if ( &list_mode = "average" ) then

     show ave(rmsd) (segid $segid and resid $resid)
     evaluate ($access=$result)

     show ave(rmsd) (segid $segid and resid $resid and &atom_main)
     evaluate ($access_main=$result)

     show ave(rmsd) (segid $segid and resid $resid and not(&atom_main))
     evaluate ($access_side=$result)

   elseif ( &list_mode = "total" ) then

     show sum(rmsd) (segid $segid and resid $resid)
     evaluate ($access=$result)

     show sum(rmsd) (segid $segid and resid $resid and &atom_main)
     evaluate ($access_main=$result)

     show sum(rmsd) (segid $segid and resid $resid and not(&atom_main))
     evaluate ($access_side=$result)

   else

     display unknown list mode
     abort

   end if

   set display=&list_outfile end
   display  $segid[a4]    $resid[a4]    $resname[a4]   $access[f10.4] \
   $access_main[f10.4]      $access_side[f10.4]
   set display=OUTPUT end
   
 end loop bfac

 do (b=rmsd) ( all ) 

 set remarks=reset end
 set remarks=accumulate end

 remarks Lee & Richards surface accessibility calculation
 remarks Surface area for each atom is in the B-factor column

 write coordinates output=&coordinate_outfile end
 
 stop