! Module file: ensavecoord
!
! 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 
!          multi-conformer structure ensembles
!
! Requirements:
!       Use symbols initilized in the module acceptinit or initave
!
!

module {ensavecoord}
(
   &ave=ave;                 {OUPUT: average rmsd values}
   &ave2=ave2;               {OUPUT: standard deviation rmsd values}
   &count=count;             {INPUT: number of accepted or trial structures}
   &nmr=nmr;                 {INPUT: nmr restraints parameters}
   &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
    
      @CNS_NMRMODULE:enssubavecoord ( ave=&ave;
                                      ave2=&ave2;
                                      nmr=&nmr;
                                      output=&output;
                                      pdb=&pdb; )
      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*&output.copies)))
                              /((&count-1)*&output.copies)))
      evaluate (&ave2.rms.atom.all=sqrt
                  ((&ave2.rms.atom.all-(&ave.rms.atom.all^2/(&count*&output.copies)))
                              /((&count-1)*&output.copies)))
   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*&output.copies))
   evaluate (&ave.rms.atom.all =&ave.rms.atom.all  / (&count*&output.copies))

   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