{+ file: contact.inp +}
{+ directory: general +}
{+ description: determine contacting residues between two sets of atoms. +}
{+ comment:
           Lists atoms involved in close contacts and the distance
           between them +}
{+ authors: Axel T. Brunger +}
{+ copyright: Yale University +}

{- 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 
   - the selections store1 through store8 are available for general use -}

{- 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";

{=========================== selections =============================}

{* selection for first atom set *}
{===>} set_select.1=(segid AAAA);

{* selection for second atom set *}
{===>} set_select.2=(segid CCCC);

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

{* distance cutoff value in Angstroms *}
{* only interactions which are closers than the specified distance
   cutoff are included *}
{===>} cutoff=3.2;

{* selection for contacting atoms *}
{* only interactions between selected atoms are included *}
{===>} contact_select=(not(name c or name n or name ca or name o));

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

{* root name for output files *}
{* list of contacts will be in: <output_root>.list *}
{===>} output_root="contact";

{===========================================================================}
{        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

 {- delete everything except selected atoms -}
 delete selection=( not ( &set_select.1 or &set_select.2 ) or not known ) end
 
 {- save original coordinates -}
 coor copy end
 
 buffer layer
    display  >>> Contact analysis between specified sets of atoms
    display  ========================================================================
    display  >>> set1= &set_select.1
    display  >>> set2= &set_select.2
    display  >>> definition of contact atoms: &contact_select
    display  >>> maximum distance cutoff between contact atoms: &cutoff
    display  ========================================================================
    display
    display  List of contacting residue pairs between set1 and set2.  The atoms
    display  that form the closest contact between the particular pair of residues 
    display  and the corresponding distance are listed. 
    display 
    display    atom in set 1      atom in set 2      distance (A)
    display  ========================================================================
 end
  
 evaluate ($cutoff_on=&cutoff-1)
 evaluate ($cutoff_off=&cutoff-0.5)
 parameter
  nbonds
   cutnb=&cutoff   wmin=1.5
   repel = 1.0          
   rexponent = 4
   irexponent = 1
   rconst = 16.0
   tolerance = 0.5  nbxmod = 5
   {- the following for consistency only, not needed for repel -}
   ctonnb=$cutoff_on ctofnb=$cutoff_off
  end
 end
   
 {- loop over all residues in set 1 -}
 
 for $atom_id in id ( tag and &set_select.1 ) loop main    
                  
    {- determine atom in other set that is closest to the
       residue in the first set -}
            
    igroup 
       interaction=( byresidue ( id $atom_id ) and &contact_select ) 
                               ( &set_select.2 and &contact_select ) 
    end
            
    do (rmsd=0) ( all ) 
    distance 
       to=( byresidue ( id $atom_id ) and &contact_select ) 
       from=( &set_select.2 and &contact_select ) 
       cuton= 0.0
       cutoff= &cutoff 
       disposition=rmsd
    end 

    {- catch all zero entries ( not selected or no interactions
                                in selected range ) -}           
    do ( rmsd=9999 ) ( attribute rmsd < 0.01 ) 
            
    show min ( rmsd ) ( &set_select.2 and &contact_select )
            
    if ($select>0) then
       evaluate ($distance_value=$result) 
    else
       evaluate ($distance_value=9999)
    end if

    if ($distance_value < 9999 ) then
    
       {- loop over all contacting residues in second set -}
       for $1 in id ( tag and &set_select.2 and 
                      byresidue ( attribute rmsd < 9999 ) ) loop ggg
       
          show min ( rmsd ) ( byresidue (id $1) and &contact_select )
          
          if ($select>0) then
             evaluate ($distance_value=$result) 
          else
             evaluate ($distance_value=9999)
          end if

          if ($distance_value < 9999 ) then
             {- this kludge gets the id of the minimum distance atom, uhh... -}
             evaluate ($distance_minus=$distance_value-0.001)
             evaluate ($distance_plus=$distance_value+0.001) 
             for $2 in id (  attribute rmsd > $distance_minus and
                             attribute rmsd < $distance_plus  ) loop ppp
                evaluate ($min_atom_id=$2) 
             end loop ppp
              
             {- store resid and resname in symbols -}                 
             show element ( resid ) ( id $min_atom_id ) 
             evaluate ($resid2=$result) 
             show element ( resname ) ( id $min_atom_id ) 
             evaluate ($resname2=$result)
             show element ( name ) ( id $min_atom_id ) 
             evaluate ($atomname2=$result) 

             {- get atom in first set which corresponds to the
                minimum distance -}
             igroup 
                interaction=( byresidue ( id $atom_id ) and &contact_select ) 
                                        ( id $min_atom_id ) 
             end
             
             do (store1=rmsd) ( all )
                  
             do (rmsd=0) ( all )
             distance 
                to=( id $min_atom_id  ) 
                from=( byresidue ( id $atom_id ) and &contact_select ) 
                cuton= 0.0
                cutoff= &cutoff 
                disposition=rmsd
             end 
                  
             {- this kludge gets the id of the minimum distance atom, uhh... -}
             evaluate ($distance_minus=$distance_value-0.001)
             evaluate ($distance_plus=$distance_value+0.001) 
             for $1 in id (  attribute rmsd > $distance_minus and
                             attribute rmsd < $distance_plus  ) loop iii
                evaluate ($min2_atom_id=$1) 
             end loop iii
             
             do (rmsd=store1) ( all ) 

             {- store resid and resname in symbols -}                 
             show element ( resid ) ( id $min2_atom_id ) 
             evaluate ($resid1=$result) 
             show element ( resname ) ( id $min2_atom_id ) 
             evaluate ($resname1=$result)
             show element ( name ) ( id $min2_atom_id ) 
             evaluate ($atomname1=$result) 

                  
             buffer layer
             display   [ $resname1[a4] $resid1[a4] $atomname1[a4] ] \
[  $resname2[a4] $resid2[a4] $atomname2[a4] ]    $distance_value
             end
           end if
                  
        end loop ggg
    end if
 
 end loop main
 
 evaluate ($outfile=&output_root + ".list") 
 
 buffer layer
    to=file=$outfile
    dump
 end

 stop