diff --git a/modules/qa/pymod/export_interaction.cc b/modules/qa/pymod/export_interaction.cc
index dfc2636cb3cb84503e7b01864f36fdbfe4c54cbb..fe61394c0b819f8d32e96dbdc1c025cc9f174ac2 100644
--- a/modules/qa/pymod/export_interaction.cc
+++ b/modules/qa/pymod/export_interaction.cc
@@ -236,11 +236,12 @@ void export_Interaction()
.def("GetCount", f2, args("partner1", "partner2", "distance"))
.def("SaveToFile", &InteractionStatistics::SaveToFile, args("file_name"))
.def("LoadFromFile", &InteractionStatistics::LoadFromFile).staticmethod("LoadFromFile")
- .def("RepairCbetaStatistics", &InteractionStatistics::RepairCbetaStatistics)
- .def("IsCBetaOnly", &InteractionStatistics::IsCBetaOnly)
;
class_<AllAtomPotentialOpts>("AllAtomPotentialOpts", init<>())
+ .def(init<float,float,float,int,float>((arg("lower_cutoff"),
+ arg("upper_cutoff"), arg("dist_bin_size"),
+ arg("sequence_sep"), arg("sigma")=0.02)))
.def_readwrite("LowerCutoff", &AllAtomPotentialOpts::lower_cutoff)
.def_readwrite("UpperCutoff", &AllAtomPotentialOpts::upper_cutoff)
.def_readwrite("DistanceBucketSize",
@@ -257,6 +258,7 @@ void export_Interaction()
.def("GetEnergyCounts", &AllAtomPotential::GetEnergyCounts)
.def("GetEnergy", three_args, args("atom_type_1", "atom_type_1", "distance"))
.def("SaveToFile", &AllAtomPotential::SaveToFile)
+ .add_property("options", make_function(&AllAtomPotential::GetOptions, return_value_policy<copy_const_reference>()))
.def("SetSequenceSeparation", &AllAtomPotential::SetSequenceSeparation)
;
register_ptr_to_python<AllAtomPotentialPtr>();
diff --git a/modules/qa/src/all_atom_potential.cc b/modules/qa/src/all_atom_potential.cc
index 606d74d2aecafa12d05e41e354b5ef863fd5251b..18a4ce5f01488a74a5d62d24553e5804f7938e10 100644
--- a/modules/qa/src/all_atom_potential.cc
+++ b/modules/qa/src/all_atom_potential.cc
@@ -36,7 +36,8 @@ namespace {
class AllAtomPotentialCalculator : public mol::EntityVisitor {
public:
AllAtomPotentialCalculator(AllAtomPotential::AllAtomEnergies& energies,
- AllAtomPotentialOpts& opts, mol::EntityView target_view):
+ AllAtomPotentialOpts& opts,
+ mol::EntityView target_view):
energies_(energies),
options_(opts),
energy_(0.0),
@@ -53,12 +54,6 @@ public:
virtual bool VisitAtom(const mol::AtomHandle& atom_handle)
{
- //TODO: The handle contains all atoms of the original loded PDB file.
- //i.e. the selection on C-beta atoms for the residue-based potential
- //is lost. Therefore a first view contains all atoms of which the energy
- //should be alculated and an additional (target-)view is provided which
- //specifies the set of atoms against which the potential is calculated
-
atom::ChemType type_a=GetAtomTypeByName(curr_aa_, atom_handle.GetName());
if (type_a!=atom::UNKNOWN) {
mol::AtomViewList nearby_atoms=target_view_.FindWithin(atom_handle.GetPos(),
@@ -75,24 +70,24 @@ public:
float d=geom::Distance(atom_handle.GetPos(),
a.GetPos());
- // avoid rounding problems: distances very close to the cut-off may
- // be rounded up the to cut-off and lead to an overflow in container energies_
+ // avoid rounding problems: distances very close to the cut-off may
+ // be rounded up the to cut-off and lead to an overflow in container
+ // energies_
d=d-0.0001; //used to prevent overflow if distance == upper_cutoff
if (d<0) d=0.0;
if (d>options_.lower_cutoff) {
- // GetHandle() necessary to retriev original index of the view (as compared to index in the view):
- if (abs(atom_handle.GetResidue().GetIndex()-a.GetResidue().GetHandle().GetIndex())<options_.sequence_sep) {
+ // GetHandle() necessary to retriev original index of the view (as
+ // compared to index in the view):
+ int abs_sep=abs(atom_handle.GetResidue().GetIndex()-
+ a.GetResidue().GetHandle().GetIndex());
+ if (abs_sep<options_.sequence_sep) {
continue;
}
energy_+=energies_.Get(type_a, type_b, d);
interaction_counts_++;
- //std::cout << type_a << "(" << atom_handle.GetQualifiedName() << ")" << " " << type_b << "(" << a.GetQualifiedName() << ")" << " : " << d << " " << energies_.Get(type_a, type_b, d) << std::endl;
}
}
-// else {
-// std::cout << "Unknown type: " << type_b << std::endl;
-// }
}
}
return false;
@@ -143,8 +138,13 @@ void AllAtomPotential::SaveToFile(const String& filename)
AllAtomPotentialPtr AllAtomPotential::LoadFromFile(const String& filename)
{
- if(!boost::filesystem::exists(filename))
- throw io::IOException("Could not open interaction potential data file.\nFile does not exist at: "+filename);
+ if (!boost::filesystem::exists(filename)) {
+ std::stringstream ss;
+ ss << "Could not open interaction potential. File '"
+ << filename << "' does not exist";
+ throw io::IOException(ss.str());
+ }
+
std::ifstream stream(filename.c_str(), std::ios_base::binary);
io::BinaryDataSource ds(stream);
@@ -222,7 +222,6 @@ void AllAtomPotential::Fill(const InteractionStatisticsPtr& stats)
options_.lower_cutoff+
options_.distance_bucket_size*0.5);
}
-// std::cout << "total_counts: " << total_counts << std::endl;
for (int i=0; i<atom::UNKNOWN; ++i) {
for (int j=0; j<atom::UNKNOWN; ++j) {
@@ -243,26 +242,6 @@ void AllAtomPotential::Fill(const InteractionStatisticsPtr& stats)
// prop = (Nd_xy / Nxy) / (Nd / Ntot)
float e=0.582*log(1+options_.sigma*t2)-0.582*log(1+options_.sigma*t2*d);
energies_.Set(Index(i, j, k), e);
-
- // potential is only calculated for Cbeta atoms:
- // copy the potential for Calpha (Calphas are selected if
- // Cbetas are not present)
- if (stats->IsCBetaOnly() == true) {
- //check if Cbeta (has counts) and not Glycin-Calpha
- if (t3 != 0) {
- if (i==3) {
- energies_.Set(Index(i, j-1, k), e);
- }
- else if (j==3) {
- energies_.Set(Index(i-1, j, k), e);
- }
- else {
- energies_.Set(Index(i-1, j, k), e);
- energies_.Set(Index(i, j-1, k), e);
- energies_.Set(Index(i-1, j-1, k), e);
- }
- }
- }
}
}
}
diff --git a/modules/qa/src/all_atom_potential.hh b/modules/qa/src/all_atom_potential.hh
index b32e8d078e461f60d3c66c382066cc35bfa6fdc3..0a1a99b7aeec1fa296645b6f7aeba81fd7a42ab5 100644
--- a/modules/qa/src/all_atom_potential.hh
+++ b/modules/qa/src/all_atom_potential.hh
@@ -37,7 +37,9 @@ struct DLLEXPORT_OST_QA AllAtomPotentialOpts {
AllAtomPotentialOpts();
- AllAtomPotentialOpts(float lower_cutoff, float upper_cutoff, float distance_bucket_size, int sequence_sep, float sigma=0.02);
+ AllAtomPotentialOpts(float lower_cutoff, float upper_cutoff,
+ float distance_bucket_size, int sequence_sep,
+ float sigma=0.02);
public:
/// \brief weight factor
diff --git a/modules/qa/src/interaction_statistics.cc b/modules/qa/src/interaction_statistics.cc
index 27422beaca060cecb54667c0005c927c088221c8..81399b27d725cd6626a82f3d28990e59cf5590f0 100644
--- a/modules/qa/src/interaction_statistics.cc
+++ b/modules/qa/src/interaction_statistics.cc
@@ -28,7 +28,6 @@ using namespace ost::mol::impl;
namespace ost { namespace qa {
InteractionStatistics::InteractionStatistics()
-: isCbetaStatisticsFlag_(false)
{
}
@@ -87,8 +86,7 @@ InteractionStatistics::InteractionStatistics(Real lower_cutoff,
histogram_(IntegralClassifier(atom::UNKNOWN, 0),
IntegralClassifier(atom::UNKNOWN, 0),
ContinuousClassifier(int((upper_cutoff_-lower_cutoff_)/bucket_size),
- lower_cutoff_, upper_cutoff_)),
- isCbetaStatisticsFlag_(false)
+ lower_cutoff_, upper_cutoff_))
{
upper_sqr_=upper_cutoff_*upper_cutoff_;
lower_sqr_=lower_cutoff_*lower_cutoff_;
@@ -170,49 +168,6 @@ void InteractionStatistics::Set(atom::ChemType a, atom::ChemType b,
}
-void InteractionStatistics::RepairCbetaStatistics()
-{
- int num=int((upper_cutoff_-
- lower_cutoff_)/bucket_size_);
-
- for (int k=0; k<num; ++k) {
- for (int i=1; i<atom::UNKNOWN; ++i) {
- for (int j=1; j<atom::UNKNOWN; ++j) {
- //all C_alpha counts of non-glycine atome need to be asigned to Cbeta
- if (this->GetCount(atom::ChemType(i), atom::ChemType(j),k) != 0 &&
- this->GetCount(atom::ChemType(i-1), atom::ChemType(j-1),k) != 0) {
- this->Set(atom::ChemType(i), atom::ChemType(j),k,
- this->GetCount(atom::ChemType(i), atom::ChemType(j), k) +
- this->GetCount(atom::ChemType(i-1), atom::ChemType(j-1),k));
- this->Set(atom::ChemType(i-1), atom::ChemType(j-1),k,0);
- }
-
- if (this->GetCount(atom::ChemType(i), atom::ChemType(j),k) != 0 &&
- this->GetCount(atom::ChemType(i), atom::ChemType(j-1),k) != 0) {
- this->Set(atom::ChemType(i), atom::ChemType(j), k,
- this->GetCount(atom::ChemType(i), atom::ChemType(j),k) +
- this->GetCount(atom::ChemType(i), atom::ChemType(j-1),k));
- this->Set(atom::ChemType(i), atom::ChemType(j-1),k,0);
- }
-
- if (this->GetCount(atom::ChemType(i), atom::ChemType(j),k) != 0 &&
- this->GetCount(atom::ChemType(i-1), atom::ChemType(j),k) != 0) {
- this->Set(atom::ChemType(i), atom::ChemType(j),k,
- this->GetCount(atom::ChemType(i), atom::ChemType(j), k) +
- this->GetCount(atom::ChemType(i-1), atom::ChemType(j),k));
- this->Set(atom::ChemType(i-1), atom::ChemType(j),k,0);
- }
-
- }
- }
- }
-
- isCbetaStatisticsFlag_ = true;
-
-
-}
-
-
template <typename DS>
void InteractionStatistics::Serialize(DS& ds)
{
diff --git a/modules/qa/src/interaction_statistics.hh b/modules/qa/src/interaction_statistics.hh
index 070346c7d76b4a83701c3b6cdebb862b743adb86..ca66774bdfd4f1ae47c490fdb169d70cb791062c 100644
--- a/modules/qa/src/interaction_statistics.hh
+++ b/modules/qa/src/interaction_statistics.hh
@@ -95,15 +95,6 @@ public:
/// account.
void Set(atom::ChemType a, atom::ChemType b, int distance_bin, int counts);
-
- /// for Cbeta potentials:
- /// in the initial selection Cbetas are counted and, if a residue
- /// has no Cbeta, the Calpha is taken instead.
- /// Afterwards the statistics needs to be corrected:
- /// Calpha counts are Cbeta counts ...
- void RepairCbetaStatistics();
-
-
/// \internal
template <typename DS>
void Serialize(DS& ds);
@@ -111,9 +102,6 @@ public:
virtual bool VisitResidue(const mol::ResidueHandle& r);
virtual bool VisitAtom(const mol::AtomHandle& a);
- inline bool IsCBetaOnly() {return isCbetaStatisticsFlag_;}
-
-
private:
InteractionStatistics();
@@ -128,8 +116,6 @@ private:
mol::EntityView view_b_;
AminoAcid amino_acid_;
InteractionHistogram histogram_;
-
- bool isCbetaStatisticsFlag_;
};
}}