! Module file: rmspairwise
!
! CNS MODULE
! **********
!
! Authors: Gregory L. Warren and Axel T. Brunger
!
! copyright Yale University
!
! version 03/05/98
!
! Function:
!          Calculates the pairwise rmsd values for structure
!          ensembles
!
! Requirements:
!       Use symbols initilized in the module acceptinit 
!
!

module { rmspairwise }
(
   &ave=ave;                 {OUPUT: pairwise average rmsd values}
   &ave2=ave2;               {OUPUT: pairwise standard deviation rmsd values}
   &count=count;             {INPUT: number of accepted or trial structures}
   &output=output;           {INPUT: file name prefix symbols}
   &pdb=pdb;                 {INPUT: atom selections}
)

checkversion 1.3

set message ? end
evaluate ($message_old=$result)
set echo ? end
evaluate ($echo_old=$result)
if ( $log_level = verbose ) then
  set echo=on message=normal end
else
  set echo=off message=off end
end if

if (&count > 0) then
   evaluate ($pair_count=0.)
   evaluate ($count1 = 0)
   while ($count1 < &count) loop par1
      evaluate ($count1 = $count1 +1)
      evaluate ($filename1=&output.prefix+encode($count1)+".pdb")
      set remarks=reset end
      coor disp=comp @@$filename1
      
      evaluate ($count2 = 0)
      while ($count2 < &count) loop par2
         evaluate ($count2 = $count2 + 1)
         if ($count1 < $count2) then
            evaluate ($filename2=&output.prefix+encode($count2)+".pdb")
            set remarks=reset end
            coor @@$filename2

            coor fit selection=( &pdb.atom.select ) end
            coor rms selection=( &pdb.atom.select ) end
            evaluate (&ave.pair.atom.slct=&ave.pair.atom.slct + $rms)
            evaluate (&ave2.pair.atom.slct=&ave2.pair.atom.slct + $rms^2)
            evaluate (&ave.min.pair.slct=min($rms,&ave.min.pair.slct))
            evaluate (&ave.max.pair.slct=max($rms,&ave.max.pair.slct))
         
            coor rms selection=( not chemical h* ) end
            evaluate (&ave.pair.atom.all=&ave.pair.atom.all + $rms)
            evaluate (&ave2.pair.atom.all=&ave2.pair.atom.all + $rms^2)
            evaluate (&ave.min.pair.all=min($rms,&ave.min.pair.all))
            evaluate (&ave.max.pair.all=max($rms,&ave.max.pair.all))
         
            evaluate ($pair_count=$pair_count + 1)
         end if
      end loop par2
      
   end loop par1
   
   if ($pair_count > 1) then
      evaluate (&ave2.pair.atom.slct=sqrt
                   ((&ave2.pair.atom.slct-(&ave.pair.atom.slct^2/$pair_count))
                                 /($pair_count-1)))
      evaluate (&ave2.pair.atom.all=sqrt
                   ((&ave2.pair.atom.all-(&ave.pair.atom.all^2/$pair_count))
                                 /($pair_count-1)))
   else
      evaluate (&ave2.pair.atom.slct=0.)
      evaluate (&ave2.pair.atom.all=0.)
   end if
   
   evaluate (&ave.pair.atom.slct=&ave.pair.atom.slct/$pair_count)
   evaluate (&ave.pair.atom.all=&ave.pair.atom.all/$pair_count)
else
   evaluate (&ave2.pair.atom.slct=0.)
   evaluate (&ave2.pair.atom.all=0.)
   evaluate (&ave.pair.atom.slct=0.)
   evaluate (&ave.pair.atom.all=0.)
end if   
set message=$message_old echo=$echo_old end