{+ file: difference_distance.inp +}
{+ directory: general +}
{+ description: Calculate a difference distance matrix between 2 structures +}
{+ comment: 
            outputs a Mathematica matrix of differences in
            distances between the 2 structures +}
{+ authors: Axel T. Brunger, and Paul D. Adams +}
{+ 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(

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

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

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

{* coordinate file *}
{===>} main_coordinate_infile="amy_anneal.pdb";

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

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

{* select atom type to calculate distances between *}
{* typically this will be CA for proteins or P for DNA/RNA *}
{===>} atom_type=(name CA);

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

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

{===========================================================================}
{        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 disp=comp @&comp_coordinate_infile

 coordinates @&main_coordinate_infile

 set display=&list_outfile end

 display (* Difference distance matrix between 2 structures: *)
 display (*   structure file= &STRIP%structure_infile *)
 display (*   main coordinate set= &STRIP%main_coordinate_infile *)
 display (*   reference coordinate set= &STRIP%comp_coordinate_infile *)
 display

 display matrix={

 set message=off end
 set echo=off end

 evaluate ($first=true)

 for $id in id ( &atom_select and &atom_type ) loop ca1

   if ( $first = true ) then
     display {
     evaluate ($first=false)
   else
     display ,{
   end if

   show (segid) (id $id)
   evaluate ($segid1=$result)
   show (resid) (id $id)
   evaluate ($resid1=$result)

   show (x) (id $id)
   evaluate ($xm1=$result)
   show (y) (id $id)
   evaluate ($ym1=$result)
   show (z) (id $id)
   evaluate ($zm1=$result)

   show (xcomp) (id $id)
   evaluate ($xc1=$result)
   show (ycomp) (id $id)
   evaluate ($yc1=$result)
   show (zcomp) (id $id)
   evaluate ($zc1=$result)
   
   evaluate ($flag=true)

   for $id in id ( &atom_select and &atom_type ) loop ca2

     show (segid) (id $id)
     evaluate ($segid2=$result)
     show (resid) (id $id)
     evaluate ($resid2=$result)

     show (x) (id $id)
     evaluate ($xm2=$result)
     show (y) (id $id)
     evaluate ($ym2=$result)
     show (z) (id $id)
     evaluate ($zm2=$result)

     show (xcomp) (id $id)
     evaluate ($xc2=$result)
     show (ycomp) (id $id)
     evaluate ($yc2=$result)
     show (zcomp) (id $id)
     evaluate ($zc2=$result)

     evaluate ($dist1=sqrt(($xm2-$xm1)^2+($ym2-$ym1)^2+($ym2-$ym1)^2))

     evaluate ($dist2=sqrt(($xc2-$xc1)^2+($yc2-$yc1)^2+($yc2-$yc1)^2))

     evaluate ($diff=$dist1-$dist2)
   
     if ( $flag = true ) then
       display $diff[f8.4]
       evaluate ($flag=false)
     else
       display , $diff[f8.4]
     end if

   end loop ca2

   display }

 end loop ca1

 display };

 stop