diff --git a/modules/qa/pymod/export_reduced.cc b/modules/qa/pymod/export_reduced.cc index 0c230fb7d74aa88be2531ff122cadae94aeacd83..3384c45290d7237f8a9908f2e4b5ae3586014464 100644 --- a/modules/qa/pymod/export_reduced.cc +++ b/modules/qa/pymod/export_reduced.cc @@ -23,10 +23,17 @@ #include <ost/qa/reduced_potential.hh> using namespace ost::qa; - +using namespace ost::mol; using namespace boost::python; -uint64_t (ReducedStatistics::*get_count)(AminoAcid,AminoAcid,int,int)=&ReducedStatistics::GetCount; +uint64_t (ReducedStatistics::*get_count)(AminoAcid, + AminoAcid, + int,int)=&ReducedStatistics::GetCount; +void (ReducedStatistics::*ex_h)(EntityHandle)=&ReducedStatistics::Extract; +void (ReducedStatistics::*ex_v)(EntityView)=&ReducedStatistics::Extract; + +Real (ReducedPotential::*get_h)(EntityHandle,bool)=&ReducedPotential::GetTotalEnergy; +Real (ReducedPotential::*get_v)(EntityView,bool)=&ReducedPotential::GetTotalEnergy; void export_Reduced() { class_<ReducedStatOptions>("ReducedStatOptions", init<>()) @@ -43,7 +50,8 @@ void export_Reduced() class_<ReducedStatistics, ReducedStatisticsPtr>("ReducedStatistics", no_init) .def(init<Real,Real,uint,uint,uint>((arg("l_cutoff"), arg("u_cutoff"), arg("num_ang_bins"), arg("num_dist_bins"), arg("sequence_sep")))) - .def("Extract", &ReducedStatistics::Extract) + .def("Extract", ex_h, arg("ent")) + .def("Extract", ex_v, arg("ent")) .add_property("options", make_function(&ReducedStatistics::GetOptions, return_value_policy<copy_const_reference>())) .def("Save", &ReducedStatistics::Save) @@ -55,7 +63,9 @@ void export_Reduced() .def("Create", &ReducedPotential::Create).staticmethod("Create") .def("Save", &ReducedPotential::Save) .def("GetEnergy", &ReducedPotential::GetEnergy) - .def("GetTotalEnergy", &ReducedPotential::GetTotalEnergy, + .def("GetTotalEnergy", get_h, + (arg("ent"), arg("normalize")=true)) + .def("GetTotalEnergy", get_v, (arg("ent"), arg("normalize")=true)) .add_property("options", make_function(&ReducedPotential::GetOptions, return_value_policy<copy_const_reference>())) diff --git a/modules/qa/src/impl/reduced_impl.cc b/modules/qa/src/impl/reduced_impl.cc index 8994fada4c0d00b6bc074eb5f331a54440427b59..a3778cb39f4f2377e65bbed8795c3f61aa0d59ae 100644 --- a/modules/qa/src/impl/reduced_impl.cc +++ b/modules/qa/src/impl/reduced_impl.cc @@ -20,39 +20,65 @@ bool ReducedPotentialImpl::VisitResidue(const mol::ResidueHandle& res) if (!this->GetCAlphaCBetaPos(res, ca_pos_one, cb_pos_one)) { return false; } + + if (ent_) { + // we got a full entity handle. + mol::AtomHandleList within=ent_.FindWithin(ca_pos_one, + opts_.upper_cutoff-0.00001); + for (mol::AtomHandleList::const_iterator + i=within.begin(), e=within.end(); i!=e; ++i) { + if (i->GetName()=="CA") { + this->HandleResidue(aa_one, ca_pos_one, cb_pos_one, index, *i); + } - mol::AtomHandleList within=ent_.FindWithin(ca_pos_one, + } + return false; + } + mol::AtomViewList within=view_.FindWithin(ca_pos_one, opts_.upper_cutoff-0.00001); - for (mol::AtomHandleList::const_iterator + for (mol::AtomViewList::const_iterator i=within.begin(), e=within.end(); i!=e; ++i) { - mol::ResidueHandle res_two=i->GetResidue(); - if (res_two.GetIndex()-index<opts_.sequence_sep) { - continue; - } - if (i->GetName()!="CA" || !res_two.IsPeptideLinking()) { - continue; + if (i->GetName()=="CA") { + this->HandleResidue(aa_one, ca_pos_one, cb_pos_one, index, i->GetHandle()); } + } + return false; +} - AminoAcid aa_two=OneLetterCodeToAminoAcid(res_two.GetOneLetterCode()); - if (aa_two==Xxx) { - continue; - } - geom::Vec3 ca_pos_two; - geom::Vec3 cb_pos_two; - if (!this->GetCAlphaCBetaPos(res_two, ca_pos_two, cb_pos_two)) { - continue; - } +void ReducedPotentialImpl::HandleResidue(AminoAcid aa_one, + const geom::Vec3& ca_pos_one, + const geom::Vec3& cb_pos_one, + uint index_one, + const mol::AtomHandle& ca_two) +{ + + mol::ResidueHandle res_two=ca_two.GetResidue(); + if (res_two.GetIndex()-index_one<opts_.sequence_sep) { + return; + } + if (!res_two.IsPeptideLinking()) { + return; + } - Real dist=geom::Length(ca_pos_one-ca_pos_two); - if (dist<opts_.lower_cutoff) { - continue; - } - Real angle=geom::Angle(cb_pos_one-ca_pos_one, cb_pos_two-ca_pos_two); - this->OnInteraction(aa_one, aa_two, dist, angle); + AminoAcid aa_two=OneLetterCodeToAminoAcid(res_two.GetOneLetterCode()); + if (aa_two==Xxx) { + return; } - return false; + geom::Vec3 ca_pos_two; + geom::Vec3 cb_pos_two; + if (!this->GetCAlphaCBetaPos(res_two, ca_pos_two, cb_pos_two)) { + return; + } + + Real dist=geom::Length(ca_pos_one-ca_pos_two); + if (dist<opts_.lower_cutoff) { + return; + } + Real angle=geom::Angle(cb_pos_one-ca_pos_one, cb_pos_two-ca_pos_two); + + this->OnInteraction(aa_one, aa_two, dist, angle); } bool ReducedPotentialImpl::GetCAlphaCBetaPos(const mol::ResidueHandle& res, diff --git a/modules/qa/src/impl/reduced_impl.hh b/modules/qa/src/impl/reduced_impl.hh index 98651360bfdd1ef82994a7a2ce44ed15f7772b35..43c312ea2d792cc9e1c0dae77668904b3d7c01c1 100644 --- a/modules/qa/src/impl/reduced_impl.hh +++ b/modules/qa/src/impl/reduced_impl.hh @@ -35,6 +35,10 @@ public: ReducedPotentialImpl(const ReducedStatOptions& opts, mol::EntityHandle ent): opts_(opts), ent_(ent) { } + + ReducedPotentialImpl(const ReducedStatOptions& opts, mol::EntityView ent): + opts_(opts), view_(ent) + { } virtual bool VisitResidue(const mol::ResidueHandle& res); @@ -44,9 +48,13 @@ public: virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two, Real dist, Real angle)=0; -private: +private: + void HandleResidue(AminoAcid aa_one, const geom::Vec3& ca_pos_one, + const geom::Vec3& cb_pos_one, + uint index_one, const mol::AtomHandle& ca_two); ReducedStatOptions opts_; mol::EntityHandle ent_; + mol::EntityView view_; }; }}} diff --git a/modules/qa/src/reduced_potential.cc b/modules/qa/src/reduced_potential.cc index 78df5c9022077f1f0a25108b56143c412dd6a785..1799c18e4ee805411a11cd3f9a0f2ee702914d49 100644 --- a/modules/qa/src/reduced_potential.cc +++ b/modules/qa/src/reduced_potential.cc @@ -35,7 +35,11 @@ public: impl::ReducedPotentialImpl(opts, ent), energies_(energies), energy_(0.0), norm_(norm), count_(0) { } - + ReducedEnergiesCalc(const ReducedStatOptions& opts, ReducedEnergies& energies, + mol::EntityView ent, bool norm): + impl::ReducedPotentialImpl(opts, ent), + energies_(energies), energy_(0.0), norm_(norm), count_(0) + { } virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two, Real dist, Real angle) { @@ -138,5 +142,12 @@ Real ReducedPotential::GetTotalEnergy(ost::mol::EntityHandle ent, bool norm) return calc.GetEnergy(); } +Real ReducedPotential::GetTotalEnergy(ost::mol::EntityView ent, bool norm) +{ + ReducedEnergiesCalc calc(opts_, energies_, ent, norm); + ent.Apply(calc); + return calc.GetEnergy(); +} + }} diff --git a/modules/qa/src/reduced_potential.hh b/modules/qa/src/reduced_potential.hh index 534892678bed5b84d10c36b43688ce7ce71479f5..0f90debce307e003657877c1d9caec0c56b75f2a 100644 --- a/modules/qa/src/reduced_potential.hh +++ b/modules/qa/src/reduced_potential.hh @@ -49,7 +49,7 @@ public: Real GetTotalEnergy(ost::mol::EntityHandle ent, bool normalize=true); - + Real GetTotalEnergy(ost::mol::EntityView ent, bool normalize=true); static bool GetCAlphaCBetaPos(const mol::ResidueHandle& res, geom::Vec3& ca_pos, geom::Vec3& cb_pos); diff --git a/modules/qa/src/reduced_statistics.cc b/modules/qa/src/reduced_statistics.cc index 9c20a9cd20d3aa9f0c026417eae46479c4b4cdb5..723fe29f9ec6a4ddb65d260dd54ee2087e643b66 100644 --- a/modules/qa/src/reduced_statistics.cc +++ b/modules/qa/src/reduced_statistics.cc @@ -41,6 +41,12 @@ public: histo_(histo) { } + ReducedStatExtractor(const ReducedStatOptions& opts, ReducedHistogram& histo, + mol::EntityView ent): + impl::ReducedPotentialImpl(opts, ent), + histo_(histo) + { } + virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two, Real dist, Real angle) { @@ -130,4 +136,17 @@ void ReducedStatistics::Extract(mol::EntityHandle ent) ent.Apply(extractor); } +void ReducedStatistics::Extract(mol::EntityView ent) +{ + if (ent.GetChainCount()!=1) { + std::stringstream ss; + ss << "Expected exactly one chain, but entity has " + << ent.GetChainCount() << "chains"; + throw std::runtime_error(ss.str()); + } + ReducedStatExtractor extractor(opts_, histo_, ent); + ent.Apply(extractor); +} + + }} \ No newline at end of file diff --git a/modules/qa/src/reduced_statistics.hh b/modules/qa/src/reduced_statistics.hh index 419ac550f4b00c6ac216def800ba9860da295dbe..7d1a24800cc85f0363181705b1dabd2c579d7d7e 100644 --- a/modules/qa/src/reduced_statistics.hh +++ b/modules/qa/src/reduced_statistics.hh @@ -82,9 +82,12 @@ public: const ReducedStatOptions& GetOptions() const { return opts_; } - /// \brief extract the statistics from the given entity + /// \brief extract the statistics from the given entity (handle) void Extract(mol::EntityHandle ent); + /// \brief extract statistics from given entity (view) + void Extract(mol::EntityView ent); + void Save(const String& filename);