expose additional methods to calculate mean force energies

We now grant low-level access to the energy calculation, but without
exposing the implementation details. These methods are all completely
reentrant and can thus be used in a multi-threaded environment.

modules/qa/src/packing_potential.cc

@@ -187,10 +187,7 @@ bool PackingPotential::VisitAtom(const mol::AtomHandle& atom)
        e=views_.end(); i!=e; ++i) {
     count+=i->FindWithin(atom.GetPos(), options_.cutoff).size();
   }
-
-  count=count>options_.max_counts ? options_.max_counts : count;
-
-  energy_+=energies_.Get(aa, count/options_.bucket_size);
+  energy_+=this->GetPackingEnergy(aa, count);
   energy_counts_++;
   return false;
 }

modules/qa/src/packing_potential.hh

@@ -88,6 +88,13 @@ public:
 
   int GetEnergyCounts();
   
+  
+  Real GetPackingEnergy(AminoAcid aa, int count) const
+  {
+    count=count>options_.max_counts ? options_.max_counts : count;
+    return energies_.Get(aa, count/options_.bucket_size);
+  }
+  
   template <typename DS>
   void Serialize(DS& ds);
   

modules/qa/src/torsion_potential.cc

@@ -16,6 +16,7 @@
 // along with this library; if not, write to the Free Software Foundation, Inc.,
 // 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 //------------------------------------------------------------------------------
+
 #include "torsion_potential.hh"
 #include "amino_acids.hh"
 
@@ -39,10 +40,9 @@ namespace {
 
 class TorsionEnergyCalc : public mol::EntityVisitor {
 public:
-  TorsionEnergyCalc(TorsionPotential::TorsionEnergies& energies, 
+  TorsionEnergyCalc(const TorsionPotentialPtr& pot, 
                     TorsionPotentialOpts opts):
-    energies_(energies), energy_(0.0),
-    opts_(opts), num_torsions_(0)
+    pot_(pot), energy_(0.0), num_torsions_(0)
   {
   }
   
@@ -75,24 +75,16 @@ public:
     Real next_phi=next_phit.GetAngle()*180/M_PI - 0.001; if (next_phi<-180) next_phi=-180;
     Real next_psi=next_psit.GetAngle()*180/M_PI - 0.001; if (next_psi<-180) next_psi=-180;
 
+    
     // calculate position of the amino acid in the alphabet
-    int icenter=this->GetAAIndex(ca);
-
-    if (icenter >-1) {
-      energy_+=energies_.Get(icenter, prev_phi, prev_psi,
-                             central_phi, central_psi,
-                             next_phi, next_psi);
-      num_torsions_++;
-    } 
-    else {
-      LOG_INFO("Amino acid not found in alphabets...");
-    }
+    energy_+=pot_->GetTorsionEnergy(ca, prev_phi, prev_psi, central_phi, 
+                                    central_psi, next_phi, next_psi);
+    ++num_torsions_;
     return false;
   }
 
   Real GetEnergy() const {
     return energy_;
-    //return num_torsions_==0?0:energy_/num_torsions_;
   }
 
   int GetEnergyCounts() const {
@@ -100,29 +92,11 @@ public:
   }
 
 private:
-  int GetAAIndex(AminoAcid aa)
-  {
-    return this->GetIndex(aa, opts_.alphabet);
-  }
-
-  
-  int GetIndex(AminoAcid aa, AminoAcidAlphabet& alpha) {
-    AminoAcidAlphabet::iterator i=alpha.begin();
-    for (int index=0; i!=alpha.end(); ++i, ++index) {
-      if ((*i).Contains(aa)) {
-        return index;        
-      }
-
-    }
-    return -1;
-  }
-private:
-  TorsionPotential::TorsionEnergies& energies_;
+  TorsionPotentialPtr pot_;
   AminoAcid prev_;
   AminoAcid center_;
   mol::ResidueHandle cr_;
   Real energy_;  
-  TorsionPotentialOpts opts_;  
   int    num_torsions_;  
 };
 
@@ -153,6 +127,18 @@ TorsionPotentialOpts::TorsionPotentialOpts():
 }
 
 
+int TorsionPotential::GetAAIndex(AminoAcid aa) const 
+{
+  AminoAcidAlphabet::const_iterator i=options_.alphabet.begin();
+  for (int index=0; i!=options_.alphabet.end(); ++i, ++index) {
+    if ((*i).Contains(aa)) {
+      return index;        
+    }
+  }
+  return -1;
+}
+
+
 TorsionPotentialPtr TorsionPotential::Create(TorsionStatisticsPtr statistics,
                                              const TorsionPotentialOpts& opts,
                                              bool calculate_average_energy_flag)
@@ -225,7 +211,7 @@ void TorsionPotential::SaveToFile(const String& path)
 
 Real TorsionPotential::GetTotalEnergy(mol::EntityHandle entity) 
 {
-  TorsionEnergyCalc c(energies_, options_);  
+  TorsionEnergyCalc c(this->shared_from_this(), options_);  
   entity.Apply(c);  
   num_torsions_ = c.GetEnergyCounts();
   return c.GetEnergy();
@@ -233,7 +219,7 @@ Real TorsionPotential::GetTotalEnergy(mol::EntityHandle entity)
 
 Real TorsionPotential::GetTotalEnergy(mol::EntityView entity) 
 {
-  TorsionEnergyCalc c(energies_, options_);
+  TorsionEnergyCalc c(this->shared_from_this(), options_);
   entity.Apply(c);
   num_torsions_ = c.GetEnergyCounts();
   return c.GetEnergy();

modules/qa/src/torsion_potential.hh

@@ -30,6 +30,7 @@
 */
 #include <ost/qa/torsion_statistics.hh>
 #include <boost/shared_ptr.hpp>
+#include <boost/enable_shared_from_this.hpp>
 #include <vector>
 #include <boost/scoped_ptr.hpp>
 
@@ -67,7 +68,8 @@ typedef boost::shared_ptr<TorsionPotential> TorsionPotentialPtr;
 /// \brief Torsion potential
 /// 
 /// The torsion potential class is parametrisable by TorsionPotentialOpts.
-class DLLEXPORT_OST_QA TorsionPotential {
+class DLLEXPORT_OST_QA TorsionPotential : 
+   public boost::enable_shared_from_this<TorsionPotential> {
 public:
   /// \brief  create new torsion potential with the given torsion statistics 
   ///         and options
@@ -89,6 +91,18 @@ public:
 
   /// \brief retrieve total number of energy local (i.e. valid residues)
   int GetEnergyCounts() const;
+  
+  
+  
+  Real GetTorsionEnergy(AminoAcid central_aa, Real prev_phi, Real prev_psi,
+                        Real central_phi, Real central_psi,  
+                        Real next_phi, Real next_psi) const {
+    int icenter=this->GetAAIndex(central_aa);
+    return energies_.Get(icenter, prev_phi, prev_psi,
+                         central_phi, central_psi,
+                         next_phi, next_psi);
+  }
+                        
 
   /// \brief save torsion potential
   /// 
@@ -106,6 +120,9 @@ public:
   typedef MultiClassifier<float, int, Real, Real,
                           Real, Real, Real, Real> TorsionEnergies;
 private:
+  
+  int GetAAIndex(AminoAcid aa) const;
+  
   void Fill(const TorsionStatisticsPtr& stat,
             bool calculate_average_energy_flag);