remarks file  nmr/average.inp 
remarks Computes the average structure, atomic rms differences from the 
remarks mean for a family of structures, and average overall rms 
remarks difference between the family and the mean structure. 

parameter @TOPPAR:protein.par @TOPPAR:nucleic.par 
          @TOPPAR:axes.par end
{====>} 
structure @axis_600.psf @axis_700.psf @axis_new.psf @sry_wt.psf @sry_dna.psf end                     {*Read the structure file.*}


{====>}       {*"Backbone" selection--excluding the more flexible ends of the *}
              {* sequence.*}
vector idend ( store9 ) ( (name ca or name n or name c) and resid 4:81)  
vector idend ( store8) (not hydrogen  and (resid 105:112 or resid 117:124))

   {* The first stage consists of computing the mean structure.*}

{====>}                 {*Loop through the family of 10 accepted structures.*}
evaluate ($end_count=3)   

eval ($nfile=0)
vector do (store1=0) (all)
vector do (store2=0) (all)
vector do (store3=0) (all)
vector do (store4=0) (all)

evaluate ($count = 0)
while ($count < $end_count ) loop main 
   evaluate ($count=$count+1)      

{====>}                     {*This is the name of the family of structures.*}
   evaluate ($filename="final_"+encode($count)+".sa")

   coor @@$filename  

   if ($count=1) then
      coor copy end              {*Store first structure in comparison set.*} 
   end if
   coor sele=( recall9 ) fit end  {*fit current to first (stored in copy)*}
   vector do (store1=store1+x) (all)
   vector do (store2=store2+y) (all)
   vector do (store3=store3+z) (all)
   vector do (store4=store4+x*x+y*y+z*z) (all)
   eval ($nfile=$nfile+1)
end loop main

vector do (x = store1 / $nfile) (all)
vector do (y = store2 / $nfile) (all)
vector do (z = store3 / $nfile) (all)
vector do (bcomp=sqrt(max(0,store4/$nfile-(x**2+y**2+z**2)))) (all)
{* <x-xave>^2 equals <x>^2-<xave>^2, so bcomp has RMS difference from average *}
{* structure for each atom. *}

        {*The second stage consists of computing an overall rms difference.*}

evaluate ($ave_rmsd_all=0.)
evaluate ($ave_rmsd_back_core=0.)
evaluate ($ave_rmsd_back_all=0.0)
evaluate ($ave_rmsd_ca=0.0)
evaluate ($ave_rmsd_dna=0.)

coor copy end

evaluate ($count = 0)
while ($count < $end_count ) loop main 
   evaluate ($count=$count+1)      
                          
{====>}                     {*This is the name of the family of structures.*}
   evaluate ($filename="final_"+encode($count)+".sa")

   coor @@$filename 
   coor fit sele=( recall9 ) end
   coor rms selection=( recall9 )end
   evaluate ($ave_rmsd_back_core=$ave_rmsd_back_core + $result)
   coor rms selection=( not hydrogen )end
   evaluate ($ave_rmsd_all =$ave_rmsd_all + $result)
   coor rms selection = (name c or name n or name ca) end
   evaluate ($ave_rmsd_back_all =$ave_rmsd_back_all + $result)
   coor rms selection= (recall8) end
   evaluate ($ave_rmsd_dna=$ave_rmsd_dna + $result)
end loop main

evaluate ($count = 0)
while ($count < $end_count ) loop main 
   evaluate ($count=$count+1)      
                          
{====>}                     {*This is the name of the family of structures.*}
   evaluate ($filename="final_"+encode($count)+".sa")

   coor @@$filename 
   coor fit sele=( name ca ) end
   coor rms selection=( name ca )end
   evaluate ($ave_rmsd_ca=$ave_rmsd_ca + $result)
end loop main



evaluate ($ave_rmsd_back_core=$ave_rmsd_back_core / $nfile)
evaluate ($ave_rmsd_all =$ave_rmsd_all  / $nfile)
evaluate ($ave_rmsd_back_all=$ave_rmsd_back_all /$nfile)
evaluate ($ave_rmsd_ca=$ave_rmsd_ca / $nfile)
evaluate ($ave_rmsd_dna=$ave_rmsd_dna / $nfile)
display ave. rms diff. to the mean struct. for non-h atoms= $ave_rmsd_all
display ave. rms diff. to the mean struct. for the backbone= $ave_rmsd_back_all

   {* Finally, the average structure and RMSDs are written to a file.*}
coor swap end
vector do (b=bcomp) ( all )

remarks unminimized average over $nfile files
remarks ave. rms diff. to the mean struct. for non-h atoms= $ave_rmsd_all
remarks ave. rms diff. to the mean struct. for the central backbone= $ave_rmsd_back_core
remarks ave. rms diff to mean struct for all backbone= $ave_rmsd_back_all
remarks ave rms diff for CA atoms = $ave_rmsd_ca
remarks ave rms diff to mean struct for DNA = $ave_rmsd_dna
remarks b array (last column) is the rms difference from the mean

{====>}            {*Write average coordinates and RMSDs to specified file.*}
write coordinates output=average.pdb end

delete selection = (resname ANI) end

flags
 exclude * include bonds angles impropers
end

mini powell nstep 1000 nprint 100 end
  {* Minimize with only the covalent constraints. *}

write coordinates output=average_min.pdb end


stop