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

module {avecoord}
(
   &ave=ave;                 {OUPUT: average rmsd values}
   &ave2=ave2;               {OUPUT: 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 ($count1 = 0)
   while ($count1 < &count ) loop aver
      evaluate ($count1=$count1+1)
      evaluate ($filename=&output.prefix+encode($count1)+".pdb")
      set remarks=reset end
      coor @@$filename
      coor fit selection=( &pdb.atom.select ) end
      coor rms selection=( &pdb.atom.select ) end
      evaluate (&ave.rms.atom.slct =&ave.rms.atom.slct  + $rms)
      evaluate (&ave2.rms.atom.slct=&ave2.rms.atom.slct + $rms^2)
      coor rms selection=( not chemical h* ) end
      evaluate (&ave.rms.atom.all  =&ave.rms.atom.all  + $rms)
      evaluate (&ave2.rms.atom.all =&ave2.rms.atom.all + $rms^2)
      do (store1=store1+x) (all)
      do (store2=store2+y) (all)
      do (store3=store3+z) (all)
      do (store4=store4+x^2+y^2+z^2) (all)
   end loop aver

   if (&count > 1) then
      evaluate (&ave2.rms.atom.slct=sqrt
                  ((&ave2.rms.atom.slct-(&ave.rms.atom.slct^2/&count))
                              /(&count-1)))
      evaluate (&ave2.rms.atom.all=sqrt
                  ((&ave2.rms.atom.all-(&ave.rms.atom.all^2/&count))
                              /(&count-1)))
   else
      evaluate (&ave2.rms.atom.slct=0.)
      evaluate (&ave2.rms.atom.all=0.)
   end if

   evaluate (&ave.rms.atom.slct=&ave.rms.atom.slct / &count)
   evaluate (&ave.rms.atom.all =&ave.rms.atom.all  / &count)

   do (x = store1 / &count) (all)
   do (y = store2 / &count) (all)
   do (z = store3 / &count) (all)
   do (bcomp =sqrt(max(0,store4/&count-(x^2+y^2+z^2)))) (all)
else
   evaluate (&ave2.rms.atom.slct=0.)
   evaluate (&ave2.rms.atom.all=0.)
   evaluate (&ave.rms.atom.slct=0.)
   evaluate (&ave.rms.atom.all =0.)
end if

set message=$message_old echo=$echo_old end