{+ file: rms_matrix.inp +}
{+ directory: general +}
{+ description: Compute rms differences between all pairs of an ensemble of coordinates +}
{+ comment: 
            The resulting file with the rms differences can be analyzed with the cluster
            analysis program cluster_struc +}
{+ 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 -}

{- begin block parameter definition -} define(

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

{* coordinate file *}
{* if an ensemble is read, only specify the "root" file name here,
   i.e., "refine" if the coordinate files are called 
   refine_1.pdb, refine_2.pdb, refine_3.pdb ... *}
{===>} coordinate_infile="amy_ensemble";

{* coordinate ensemble count *}
{* if greater than 1 a coordinate ensemble is read.  All coordinate files 
   must have the same topology (e.g., annealing refinements repeats with
   different random number seeds) *}
{===>} coordinate_ensemble_count=3;

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

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

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

{* select atom(s) for rms calculation *}
{===>} atom_select=( resid 2:30 );

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

{* output list file *}
{* contains triangular rmsd matrix.  Can be read processed by the program cluster_struc *}
{===>} list_outfile="rms_matrix.list";

{* temporary trajectory file *}
{===>} trajectory="rms_matrix.crd";

{===========================================================================}
{        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 (&coordinate_ensemble_count > 1) then
    {- use first coordinate file to auto generate MTF if requested -}
    evaluate ($input_coor = &coordinate_infile+"_"+encode(1)+".pdb")
 else
    {- use single coordinate file -}
    evaluate ($input_coor = &coordinate_infile)
 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=$input_coor;
                                  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
   coor @$input_coor

 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=off echo=off end

 set display=&list_outfile end
 
 evaluate ($nstruc_i = 0)

 eval ($nstruc_total=&coordinate_ensemble_count) 

 eval ($trajfile=&trajectory)

 while ($nstruc_i < $nstruc_total) loop traj

  evaluate ($nstruc_i = $nstruc_i + 1)

  evaluate ($filename = &coordinate_infile+"_"+encode($nstruc_i)+".pdb")
  
  fileexist $filename end
  
  if($result=true) then
      coor @@$filename
  else 
      do (x=0) ( all )
      do (y=0) ( all )
      do (z=0) ( all )
  end if

  if ($nstruc_i = 1) then
      write trajectory
         ascii true
         output= $trajfile
      end
   else
      write traj next end
   end if
 
 end loop traj

 close $trajfile end

 evaluate ($nstruc_i = 0)

 while ($nstruc_i < $nstruc_total) loop fiti

  evaluate ($nstruc_i = $nstruc_i + 1)

  read trajectory
     ascii true
     input= $trajfile
     begin $nstruc_i skip 1 stop $nstruc_i
  end
  coor copy end

  evaluate ($nstruc_j = $nstruc_i)
  while ($nstruc_j <  $nstruc_total) loop fitj
  
    if ($nstruc_j=$nstruc_i) then
    
      evaluate ($nstruc_j = $nstruc_j + 1)
      read trajectory
          ascii true
          input= $trajfile
          begin $nstruc_j skip 1 stop $nstruc_total
       end
    else
    
       evaluate ($nstruc_j = $nstruc_j + 1)
       read trajectory next end
       
    end if

    coor sele= ( &atom_select ) rms end
    
    display $nstruc_i $nstruc_j $result
  end loop fitj
  
 end loop fiti

 stop