{+ 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(

{============================ coordinates ============================}

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

{==================== molecular information ==========================}

{* topology files *}
{===>} topology_infile_1="CNS_TOPPAR:protein.top";
{===>} topology_infile_2="CNS_TOPPAR:dna-rna.top";
{===>} topology_infile_3="CNS_TOPPAR:water.top";
{===>} topology_infile_4="CNS_TOPPAR:ion.top";
{===>} topology_infile_5="CNS_TOPPAR:carbohydrate.top";
{===>} topology_infile_6="";
{===>} topology_infile_7="";
{===>} topology_infile_8="";

{* linkage files for linear, continuous polymers (protein, DNA, RNA) *}
{===>} link_infile_1="CNS_TOPPAR:protein.link";
{===>} link_infile_2="CNS_TOPPAR:dna-rna-pho.link";
{===>} link_infile_3="";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param";
{===>} parameter_infile_2="CNS_TOPPAR:dna-rna_rep.param";
{===>} parameter_infile_3="CNS_TOPPAR:water_rep.param";
{===>} parameter_infile_4="CNS_TOPPAR:ion.param";
{===>} parameter_infile_5="CNS_TOPPAR:carbohydrate.param";
{===>} parameter_infile_6="";
{===>} parameter_infile_7="";
{===>} parameter_infile_8="";

{* molecular topology file: optional (leave blank for auto generation) *}
{* 
   Auto generation of the molecular topology from the coordinates should only 
   be used if:
   (1) Each distinct protein, DNA, or RNA chain must have a separate segid 
       (or chainid if the chainid is non-blank). 
   (2) Each contiguous protein, RNA, or RNA chain must not be disrupted by 
       other types of residues or ligands.  Rather, these other residues 
       should be listed after protein, RNA/DNA chains. 
   (3) Disulphides are automatically detected based on distances between the sulfur atoms
      (must be less than 3 A apart).
   (4) Broken protein/RNA/DNA chains without terminii must be more than 2.5 A apart to be recognized as such.
   (5) N-linked glycan links are automatically recognized if the bonded atoms are less than 2.5 A apart.
   (6) Automatic generation cannot be used with alternate conformations. 
   For ligands, the user must make suitable topology and parameter files.
   For non-standard covalent linkages, the custom patch file should be used.
   Alternatively, the generate.inp or generate_easy.inp task files
   can be used to generated the mtf prior to running this task file.
    *}
{===>} structure_infile="eg1_dimer.mtf";

{* for auto generation: extra linkages and modifications by custom patches *}
{===>} patch_infile="";

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

{* output file *}
{===>} outfile="contact.list";

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

 ) {- end block parameter definition -}

 checkversion 1.3

 evaluate ($log_level=quiet)

  if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 if ( &BLANK%structure_infile = true ) then
 
    {- read topology files -}
    topology
     evaluate ($counter=1)
     evaluate ($done=false)
     while ( $done = false ) loop read
      if ( &exist_topology_infile_$counter = true ) then
         if ( &BLANK%topology_infile_$counter = false ) then
            @@&topology_infile_$counter
         end if
      else
        evaluate ($done=true)
      end if
      evaluate ($counter=$counter+1)
     end loop read
    end
    
    @CNS_XTALMODULE:mtfautogenerate (
                                  coordinate_infile=&coordinate_infile;
                                  convert=true;
                                  separate=true;
                                  atom_delete=(not known);
                                  hydrogen_flag=true;
                                  break_cutoff=2.5;
                                  disulphide_dist=3.0;
                                  carbo_dist=2.5;
                                  patch_infile=&patch_infile;
                                  O5_becomes="O";
                                 )

 else

   structure @&structure_infile end
   coordinates @&coordinate_infile

 end if

 {- read parameter files -}
 parameter
  evaluate ($counter=1)
  evaluate ($done=false)
  while ( $done = false ) loop read
   if ( &exist_parameter_infile_$counter = true ) then
      if ( &BLANK%parameter_infile_$counter = false ) then
         @@&parameter_infile_$counter
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
  end loop read
 end

 set message=normal echo=on 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
  
 buffer layer
    to=file=&outfile
    dump
 end

 stop