diff --git a/modules/qa/pymod/export_interaction.cc b/modules/qa/pymod/export_interaction.cc index 6e6e119f0a6b2a2d2a9fad3662b6815516b37935..dfc2636cb3cb84503e7b01864f36fdbfe4c54cbb 100644 --- a/modules/qa/pymod/export_interaction.cc +++ b/modules/qa/pymod/export_interaction.cc @@ -237,7 +237,7 @@ void export_Interaction() .def("SaveToFile", &InteractionStatistics::SaveToFile, args("file_name")) .def("LoadFromFile", &InteractionStatistics::LoadFromFile).staticmethod("LoadFromFile") .def("RepairCbetaStatistics", &InteractionStatistics::RepairCbetaStatistics) - .def("isCbetaStaistics", &InteractionStatistics::isCbetaStaistics) + .def("IsCBetaOnly", &InteractionStatistics::IsCBetaOnly) ; class_<AllAtomPotentialOpts>("AllAtomPotentialOpts", init<>()) diff --git a/modules/qa/src/all_atom_potential.cc b/modules/qa/src/all_atom_potential.cc index 348880447a93c9f8530bc1a4b65a146f809e2d72..606d74d2aecafa12d05e41e354b5ef863fd5251b 100644 --- a/modules/qa/src/all_atom_potential.cc +++ b/modules/qa/src/all_atom_potential.cc @@ -153,36 +153,25 @@ AllAtomPotentialPtr AllAtomPotential::LoadFromFile(const String& filename) return all_atom_p; } - -float AllAtomPotential::GetTotalEnergy(mol::EntityView view, mol::EntityView target_view) +float AllAtomPotential::GetTotalEnergy(mol::EntityView view, + mol::EntityView target_view) { AllAtomPotentialCalculator c(energies_, options_, target_view); mol::EntityHandle e=view.GetHandle(); view.Apply(c); - interaction_counts_member_=c.GetEnergyCounts(); + interaction_counts_=c.GetEnergyCounts(); return c.GetEnergy(); } -float AllAtomPotential::GetEnergy(atom::ChemType type_a, atom::ChemType type_b, float distance) -{ - return energies_.Get(type_a, type_b, distance); -} - float AllAtomPotential::GetTotalEnergy(mol::EntityView view) { AllAtomPotentialCalculator c(energies_, options_, view); mol::EntityHandle e=view.GetHandle(); view.Apply(c); - interaction_counts_member_=c.GetEnergyCounts(); + interaction_counts_=c.GetEnergyCounts(); return c.GetEnergy(); } - -int AllAtomPotential::GetEnergyCounts() const { - return interaction_counts_member_; -} - - void AllAtomPotential::SetSequenceSeparation(int seq_sep) { options_.sequence_sep=seq_sep; } @@ -258,7 +247,7 @@ void AllAtomPotential::Fill(const InteractionStatisticsPtr& stats) // potential is only calculated for Cbeta atoms: // copy the potential for Calpha (Calphas are selected if // Cbetas are not present) - if (stats->isCbetaStaistics() == true) { + if (stats->IsCBetaOnly() == true) { //check if Cbeta (has counts) and not Glycin-Calpha if (t3 != 0) { if (i==3) { diff --git a/modules/qa/src/all_atom_potential.hh b/modules/qa/src/all_atom_potential.hh index 81beef95d7a1bf50168aa309c8865789035a38bd..5f0e15b4e4691090396088435a17f11d8e6802b1 100644 --- a/modules/qa/src/all_atom_potential.hh +++ b/modules/qa/src/all_atom_potential.hh @@ -76,8 +76,11 @@ public: /// \brief extract energy of a specific interaction /// (for plotting pseudo Lennard-Jones potential). - float GetEnergy(atom::ChemType type_a, atom::ChemType type_b, float distance); - + float GetEnergy(atom::ChemType type_a, atom::ChemType type_b, + float distance) + { + return energies_.Get(type_a, type_b, distance); + } /// \brief calculate all-atom interaction between two entities. /// Two entities need to be provided: /// the atoms for which the energy should be derived and @@ -85,7 +88,7 @@ public: float GetTotalEnergy(mol::EntityView view, mol::EntityView target_view); /// \brief retrieve total number of interactions (for normalisation) - int GetEnergyCounts() const; + int GetEnergyCounts() const { return interaction_counts_; } /// \brief set different seqeunce separation than used for training void SetSequenceSeparation(int seq_sep); @@ -102,7 +105,7 @@ private: AllAtomPotentialOpts options_; AllAtomEnergies energies_; mol::EntityView target_view_; - int interaction_counts_member_; + int interaction_counts_; }; }} diff --git a/modules/qa/src/interaction_statistics.hh b/modules/qa/src/interaction_statistics.hh index a1e8eeb48e1c2eee1efb969eb74a17987968335f..070346c7d76b4a83701c3b6cdebb862b743adb86 100644 --- a/modules/qa/src/interaction_statistics.hh +++ b/modules/qa/src/interaction_statistics.hh @@ -111,7 +111,7 @@ public: virtual bool VisitResidue(const mol::ResidueHandle& r); virtual bool VisitAtom(const mol::AtomHandle& a); - inline bool isCbetaStaistics() {return isCbetaStatisticsFlag_;} + inline bool IsCBetaOnly() {return isCbetaStatisticsFlag_;} private: diff --git a/modules/qa/src/packing_potential.cc b/modules/qa/src/packing_potential.cc index e7a473fd3ac496291c2c459da974ae930fec9548..10a7de4aeb20f1642f743915bc3cab7719d7062c 100644 --- a/modules/qa/src/packing_potential.cc +++ b/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; } diff --git a/modules/qa/src/packing_potential.hh b/modules/qa/src/packing_potential.hh index dec0054f672319ab71d3ab108a907b18a46b1c9b..6d11d2e6c6ff43535e725c1f4267910edbc4a4cf 100644 --- a/modules/qa/src/packing_potential.hh +++ b/modules/qa/src/packing_potential.hh @@ -84,10 +84,17 @@ public: /// /// See documentation for PackingStatistics::Extract() Real GetTotalEnergy(const mol::EntityView& view, - const mol::EntityViewList& views); + const mol::EntityViewList& views); 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); @@ -109,7 +116,7 @@ private: // used to calculate total energy... mol::EntityView view_; mol::EntityViewList views_; - Real energy_; + Real energy_; int energy_counts_; }; diff --git a/modules/qa/src/torsion_potential.cc b/modules/qa/src/torsion_potential.cc index 651718a4f04bc0d53faf6fadf3dd1ab8eb328e0c..2ecf8c655f24de49220a00be1b09b24d1e14b856 100644 --- a/modules/qa/src/torsion_potential.cc +++ b/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(); diff --git a/modules/qa/src/torsion_potential.hh b/modules/qa/src/torsion_potential.hh index 68faeb24fc2ab71cc4c261edf6d719452dd25201..d0048e84ccce0ada40901e93dd70fea536d57d04 100644 --- a/modules/qa/src/torsion_potential.hh +++ b/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); diff --git a/modules/qa/src/torsion_statistics.hh b/modules/qa/src/torsion_statistics.hh index 15f20089a5d8e7950621d79c347c612bbd0005d0..74ad15562e4a504ebd20dfa83645a7f67002ce44 100644 --- a/modules/qa/src/torsion_statistics.hh +++ b/modules/qa/src/torsion_statistics.hh @@ -70,8 +70,8 @@ public: Real next_phi_angle, Real next_psi_angle) const; uint64_t GetCount(Real prev_phi_angle, Real prev_psi_angle, - Real central_phi_angle, Real central_psi_angle, - Real next_phi_angle, Real next_psi_angle) const; + Real central_phi_angle, Real central_psi_angle, + Real next_phi_angle, Real next_psi_angle) const; uint64_t GetCount(const AminoAcidSet& central_aa) const;