diff --git a/actions/doc/index.rst b/actions/doc/index.rst
index 450ffa57d28ef1f491876e3f468880a55741713d..39c68baf7c8c46a2be6a8409d57b058c4060ce90 100644
--- a/actions/doc/index.rst
+++ b/actions/doc/index.rst
@@ -190,3 +190,8 @@ Several flags control the modelling behaviour:
Dont do subrotamer optimization if flexible rotamer model is used
+.. option:: -f, --energy_function
+
+ The energy function to be used. Default is SCWRL4, can be any function
+ supported by :meth:`promod3.modelling.ReconstructSidechains`.
+
diff --git a/actions/pm-build-sidechains b/actions/pm-build-sidechains
index 72a0ddfc6a57d9bdca6fd07d70269d47f570abce..2a57c5165c110c4ea5bdcfbf9bc9b6d1992d58df 100644
--- a/actions/pm-build-sidechains
+++ b/actions/pm-build-sidechains
@@ -54,6 +54,8 @@ parser.add_argument('-i', '--backbone-independent', action='store_true',
parser.add_argument('-s', '--no-subrotamer-optimization', action='store_true',
help='Dont do subrotamer optimization if flexible ' +
'rotamer model is used.')
+parser.add_argument('-f', '--energy_function', action = 'store',
+ dest = "energy_function", default = "SCWRL4")
opts = parser.Parse()
@@ -80,7 +82,8 @@ modelling.ReconstructSidechains(prot, keep_sidechains=opts.keep_sidechains,
build_disulfids=opts.no_disulfids==False,
rotamer_model=rotamer_model,
rotamer_library=lib,
- optimize_subrotamers=opt_subrot)
+ optimize_subrotamers=opt_subrot,
+ energy_function = opts.energy_function)
# output
ost.PopVerbosityLevel()
diff --git a/modelling/pymod/_reconstruct_sidechains.py b/modelling/pymod/_reconstruct_sidechains.py
index b94c37c93c389fb54b3d8661777e7f720c58dced..2586ac8f7199b1ba46784b0d6b71851f346a944d 100644
--- a/modelling/pymod/_reconstruct_sidechains.py
+++ b/modelling/pymod/_reconstruct_sidechains.py
@@ -402,7 +402,8 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
rotamer_model="frm", consider_ligands=True,
rotamer_library=None, optimize_subrotamers=True,
graph_max_complexity=100000000,
- graph_initial_epsilon=0.02):
+ graph_initial_epsilon=0.02,
+ energy_function = "SCWRL4"):
'''Reconstruct sidechains for the given structure.
:param ent: Structure for sidechain reconstruction. Note, that the sidechain
@@ -445,6 +446,9 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
:param graph_intial_epsilon: Initial epsilon for
:meth:`RotamerGraph.TreeSolve`.
:type graph_intial_epsilon: :class:`float`
+ :param energy_function: What energy function to use can be any in
+ ["SCWRL4", "SCWRL3"]
+ :type energy_function: :class:`str`
'''
prof_name = 'modelling::ReconstructSidechains'
prof = core.StaticRuntimeProfiler.StartScoped(prof_name)
@@ -463,8 +467,16 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
if type(rotamer_library) is sidechain.BBDepRotamerLib:
bbdep = True
- rotamer_constructor = sidechain.SCWRL4RotamerConstructor(False)
-
+ rotamer_constructor = None
+
+ if energy_function == "SCWRL4":
+ rotamer_constructor = sidechain.SCWRL4RotamerConstructor(False)
+ if energy_function == "SCWRL3":
+ rotamer_constructor = sidechain.SCWRL3RotamerConstructor(False)
+
+ if rotamer_constructor == None:
+ raise RuntimeError("Invalid Energy function to reconstruct sidechains!")
+
# take out ligand chain and any non-peptides
prot = ent.Select("peptide=true and cname!='_'")
diff --git a/sidechain/pymod/export_rotamer_constructor.cc b/sidechain/pymod/export_rotamer_constructor.cc
index cf28f603746dd64ab52ecb8a238057ee7081f05f..fdb91bb4ac3a745b6027141ff4c974e781c9ccf4 100644
--- a/sidechain/pymod/export_rotamer_constructor.cc
+++ b/sidechain/pymod/export_rotamer_constructor.cc
@@ -21,6 +21,7 @@
#include <promod3/sidechain/rotamer.hh>
#include <promod3/sidechain/rotamer_group.hh>
#include <promod3/sidechain/scwrl4_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl3_rotamer_constructor.hh>
using namespace boost::python;
using namespace promod3::sidechain;
@@ -34,6 +35,11 @@ SCWRL4RotamerConstructorPtr WrapSCWRL4RotamerConstructorInit(bool cb_in_sidechai
return p;
}
+SCWRL3RotamerConstructorPtr WrapSCWRL3RotamerConstructorInit(bool cb_in_sidechain) {
+ SCWRL3RotamerConstructorPtr p(new SCWRL3RotamerConstructor(cb_in_sidechain));
+ return p;
+}
+
RRMRotamerGroupPtr WrapRRMGroup_res(RotamerConstructorPtr constructor,
const ost::mol::ResidueHandle& res,
RotamerID id,
@@ -444,4 +450,10 @@ void export_RotamerConstructor(){
arg("comp_lib")))
;
+
+ class_<SCWRL3RotamerConstructor, SCWRL3RotamerConstructorPtr,
+ bases<RotamerConstructor> >("SCWRL3RotamerConstructor", no_init)
+ .def("__init__", boost::python::make_constructor(&WrapSCWRL3RotamerConstructorInit))
+ ;
+
}
diff --git a/sidechain/src/CMakeLists.txt b/sidechain/src/CMakeLists.txt
index 953b809ad3108aa9ae82d9431483ecdc8e6eaedc..1c142e51635c10fc9c47042dd65b0bf78639a1e4 100644
--- a/sidechain/src/CMakeLists.txt
+++ b/sidechain/src/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SIDECHAIN_HEADERS
particle_scoring_base.hh
particle_scoring.hh
scwrl4_particle_scoring.hh
+ scwrl3_particle_scoring.hh
rotamer.hh
rotamer_cruncher.hh
rotamer_density.hh
@@ -20,6 +21,7 @@ set(SIDECHAIN_HEADERS
rotamer_constructor.hh
rotamer_lookup.hh
scwrl4_rotamer_constructor.hh
+ scwrl3_rotamer_constructor.hh
subrotamer_optimizer.hh
)
@@ -29,6 +31,7 @@ set(SIDECHAIN_SOURCES
frame.cc
particle_scoring.cc
scwrl4_particle_scoring.cc
+ scwrl3_particle_scoring.cc
rotamer.cc
rotamer_cruncher.cc
rotamer_density.cc
@@ -42,6 +45,7 @@ set(SIDECHAIN_SOURCES
rotamer_constructor.cc
rotamer_lookup.cc
scwrl4_rotamer_constructor.cc
+ scwrl3_rotamer_constructor.cc
subrotamer_optimizer.cc
)
diff --git a/sidechain/src/particle_scoring.cc b/sidechain/src/particle_scoring.cc
index c95fc241162e7acf48234422e7bf71537cb484ee..de558f3f89248257b506561b5e4ad329b9ea9a21 100644
--- a/sidechain/src/particle_scoring.cc
+++ b/sidechain/src/particle_scoring.cc
@@ -16,6 +16,7 @@
#include <promod3/core/message.hh>
#include <promod3/sidechain/particle_scoring.hh>
#include <promod3/sidechain/scwrl4_particle_scoring.hh>
+#include <promod3/sidechain/scwrl3_particle_scoring.hh>
namespace promod3{ namespace sidechain{
@@ -31,6 +32,10 @@ Real PairwiseParticleScore(PScoringParam* p_one, PScoringParam* p_two) {
return SCWRL4PairwiseScore(reinterpret_cast<SCWRL4Param*>(p_one),
reinterpret_cast<SCWRL4Param*>(p_two));
}
+ case SCWRL3: {
+ return SCWRL3PairwiseScore(reinterpret_cast<SCWRL3Param*>(p_one),
+ reinterpret_cast<SCWRL3Param*>(p_two));
+ }
default: {
throw promod3::Error("Don't know what scoring function to call between "
"particles... maybe a lazy developer forgot to "
@@ -85,6 +90,16 @@ void PScoringParamContainer::PairwiseScores(const PScoringParamContainer& other,
}
break;
}
+ case SCWRL3: {
+ scores.resize(idx.size());
+ for(uint i = 0; i < idx.size(); ++i) {
+ scores[i] =
+ SCWRL3PairwiseScore(reinterpret_cast<SCWRL3Param*>(params_[idx[i].first]),
+ reinterpret_cast<SCWRL3Param*>(other.params_[idx[i].second]));
+ }
+ break;
+
+ }
default: {
throw promod3::Error("Don't know what scoring function to call between "
"particles... maybe a lazy developer forgot to "
diff --git a/sidechain/src/particle_scoring_base.hh b/sidechain/src/particle_scoring_base.hh
index 96cdb6d97b80353e3a5331d394d6e707db9ae501..e5cbdd167589f0a72df43d01dd9abd27ff8f9926 100644
--- a/sidechain/src/particle_scoring_base.hh
+++ b/sidechain/src/particle_scoring_base.hh
@@ -25,6 +25,7 @@ namespace promod3{ namespace sidechain{
// the available particle scoring functions
enum PScoringFunction{
SCWRL4,
+SCWRL3,
INVALID_PSCORING_FUNCTION
};
diff --git a/sidechain/src/scwrl3_particle_scoring.cc b/sidechain/src/scwrl3_particle_scoring.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9cf9b04e25e3b179af274489d2b32c8e0b9edb06
--- /dev/null
+++ b/sidechain/src/scwrl3_particle_scoring.cc
@@ -0,0 +1,65 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and
+// Biozentrum - University of Basel
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+
+#include <promod3/sidechain/scwrl3_particle_scoring.hh>
+
+
+namespace promod3 { namespace sidechain {
+
+void SCWRL3Param::ApplyTransform(const geom::Mat4& t) {
+
+ // get transformed pos
+ Real a = t(0,0)*pos_[0] + t(0,1)*pos_[1] + t(0,2)*pos_[2] + t(0,3);
+ Real b = t(1,0)*pos_[0] + t(1,1)*pos_[1] + t(1,2)*pos_[2] + t(1,3);
+ Real c = t(2,0)*pos_[0] + t(2,1)*pos_[1] + t(2,2)*pos_[2] + t(2,3);
+ pos_ = geom::Vec3(a, b, c);
+}
+
+
+SCWRL3Param* SCWRL3Param::Copy() const {
+ SCWRL3Param* return_ptr = new SCWRL3Param(pos_, radius_);
+ return return_ptr;
+}
+
+
+bool SCWRL3Param::EqualTo(PScoringParam* other) const {
+
+ if(other == NULL) return false;
+ if(other->GetScoringFunction() != SCWRL3) return false;
+
+ // as the other scoring function is also SCWRL3, we can assume that
+ // the following dynamic cast is fine...
+ SCWRL3Param* p = dynamic_cast<SCWRL3Param*>(other);
+
+ return pos_ == p->pos_ &&
+ radius_ == p->radius_;
+}
+
+
+Real SCWRL3PairwiseScore(SCWRL3Param* p_one, SCWRL3Param* p_two) {
+
+ Real d = geom::Distance(p_one->pos_, p_two->pos_);
+ Real R = p_one->radius_ + p_two->radius_;
+ Real e = 0.0;
+ if(d < Real(0.8254) * R) {
+ e = 10.0;
+ } else if(d <= R) {
+ e = 57.273 * (1 - d/R);
+ }
+ return e;
+}
+
+}}//ns
diff --git a/sidechain/src/scwrl3_particle_scoring.hh b/sidechain/src/scwrl3_particle_scoring.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1208dd43931b7922ab5e1cb91538fb02f8bb3cd8
--- /dev/null
+++ b/sidechain/src/scwrl3_particle_scoring.hh
@@ -0,0 +1,52 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and
+// Biozentrum - University of Basel
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef PROMOD3_SCWRL3_PARTICLE_SCORING_HH
+#define PROMOD3_SCWRL3_PARTICLE_SCORING_HH
+
+#include <ost/geom/mat4.hh>
+#include <ost/geom/vecmat3_op.hh>
+#include <promod3/sidechain/particle_scoring_base.hh>
+
+namespace promod3{ namespace sidechain{
+
+struct SCWRL3Param : public PScoringParam {
+
+ SCWRL3Param(const geom::Vec3& pos, Real r): pos_(pos), radius_(r) { }
+
+ virtual ~SCWRL3Param() { }
+
+ virtual const geom::Vec3& GetPos() const { return pos_; }
+
+ virtual Real GetCollisionDistance() const { return radius_; }
+
+ virtual void ApplyTransform(const geom::Mat4& t);
+
+ virtual SCWRL3Param* Copy() const;
+
+ virtual bool EqualTo(PScoringParam* other) const;
+
+ virtual PScoringFunction GetScoringFunction() const { return SCWRL3; }
+
+ geom::Vec3 pos_;
+ Real radius_;
+};
+
+
+Real SCWRL3PairwiseScore(SCWRL3Param* p_one, SCWRL3Param* p_two);
+
+}} //ns
+
+#endif
diff --git a/sidechain/src/scwrl3_rotamer_constructor.cc b/sidechain/src/scwrl3_rotamer_constructor.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d5f9048bbb5cbcf1b220bff635983cdbfe94c694
--- /dev/null
+++ b/sidechain/src/scwrl3_rotamer_constructor.cc
@@ -0,0 +1,140 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and
+// Biozentrum - University of Basel
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+
+#include <promod3/sidechain/scwrl3_rotamer_constructor.hh>
+#include <promod3/sidechain/scwrl3_particle_scoring.hh>
+#include <promod3/core/runtime_profiling.hh>
+#include <promod3/core/message.hh>
+#include <promod3/core/runtime_profiling.hh>
+#include <ost/conop/conop.hh>
+
+namespace promod3 { namespace sidechain {
+
+struct SCWRL3RotamerParam : public RotamerLookupParam {
+
+ SCWRL3RotamerParam() {
+
+ // overwrite default internal_e_prefactor from base class
+ internal_e_prefactor[ost::conop::ARG] = 3.0;
+ internal_e_prefactor[ost::conop::ASN] = 3.0;
+ internal_e_prefactor[ost::conop::ASP] = 3.0;
+ internal_e_prefactor[ost::conop::GLN] = 3.0;
+ internal_e_prefactor[ost::conop::GLU] = 3.0;
+ internal_e_prefactor[ost::conop::LYS] = 3.0;
+ internal_e_prefactor[ost::conop::SER] = 3.0;
+ internal_e_prefactor[ost::conop::CYS] = 3.0;
+ internal_e_prefactor[ost::conop::MET] = 3.0;
+ internal_e_prefactor[ost::conop::TRP] = 3.0;
+ internal_e_prefactor[ost::conop::TYR] = 3.0;
+ internal_e_prefactor[ost::conop::THR] = 3.0;
+ internal_e_prefactor[ost::conop::VAL] = 3.0;
+ internal_e_prefactor[ost::conop::ILE] = 3.0;
+ internal_e_prefactor[ost::conop::LEU] = 3.0;
+ internal_e_prefactor[ost::conop::PRO] = 3.0;
+ internal_e_prefactor[ost::conop::HIS] = 3.0;
+ internal_e_prefactor[ost::conop::PHE] = 3.0;
+ }
+
+ static const SCWRL3RotamerParam& Instance() {
+ static SCWRL3RotamerParam scwrl3_param;
+ return scwrl3_param;
+ }
+};
+
+
+SCWRL3RotamerConstructor::SCWRL3RotamerConstructor(bool cb_in_sidechain):
+ RotamerConstructor(cb_in_sidechain, HEAVY_ATOM_MODE,
+ SCWRL3RotamerParam::Instance()) { }
+
+void SCWRL3RotamerConstructor::AssignInternalEnergies(RRMRotamerGroupPtr group,
+ RotamerID id,
+ uint residue_index,
+ Real phi, Real psi,
+ bool n_ter, bool c_ter) {
+
+ core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
+ "SCWRL3RotamerConstructor::AssignInternalEnergies", 2);
+
+ Real max_p = group->GetMaxP();
+ for(uint i = 0; i < group->size(); ++i){
+ RRMRotamerPtr r = (*group)[i];
+ Real internal_e_prefactor = r->GetInternalEnergyPrefactor();
+ Real probability = r->GetProbability();
+ Real e = -internal_e_prefactor * std::log(probability/max_p);
+ r->SetInternalEnergy(e);
+ }
+}
+
+
+void SCWRL3RotamerConstructor::AssignInternalEnergies(FRMRotamerGroupPtr group,
+ RotamerID id,
+ uint residue_index,
+ Real phi, Real psi,
+ bool n_ter, bool c_ter) {
+
+ core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
+ "SCWRL3RotamerConstructor::AssignInternalEnergies", 2);
+
+ Real max_p = group->GetMaxP();
+ for(uint i = 0; i < group->size(); ++i){
+ FRMRotamerPtr r = (*group)[i];
+ Real internal_e_prefactor = r->GetInternalEnergyPrefactor();
+ Real probability = r->GetProbability();
+ Real e = -internal_e_prefactor * std::log(probability/max_p);
+ r->SetInternalEnergy(e);
+ }
+}
+
+
+void SCWRL3RotamerConstructor::ParametrizeParticle(int at_idx,
+ bool is_hydrogen,
+ Particle& particle) {
+
+ if(is_hydrogen) {
+ throw promod3::Error("Expect no hydrogen in SCWRL3RotamerConstructor!");
+ }
+
+ geom::Vec3 pos;
+ Real radius = 0.0;
+
+ if(at_idx == -2) {
+ pos = terminal_o_pos_;
+ radius = 1.3;
+ } else if (at_idx == -1){
+ pos = terminal_oxt_pos_;
+ radius = 1.3;
+ } else {
+ pos = pos_buffer_->GetPos(id_, at_idx);
+ ost::conop::AminoAcid aa = RotIDToAA(id_);
+ promod3::loop::AminoAcidAtom aaa =
+ promod3::loop::AminoAcidLookup::GetInstance().GetAAA(aa, at_idx);
+ String ele = promod3::loop::AminoAcidLookup::GetInstance().GetElement(aaa);
+ if(ele == "C") {
+ radius = 1.6;
+ } else if(ele == "N") {
+ radius = 1.3;
+ } else if(ele == "O") {
+ radius = 1.3;
+ } else if (ele == "S") {
+ radius = 1.7;
+ }
+ }
+
+ SCWRL3Param* p = new SCWRL3Param(pos, radius);
+ particle.SetSParam(p);
+}
+
+}} //ns
diff --git a/sidechain/src/scwrl3_rotamer_constructor.hh b/sidechain/src/scwrl3_rotamer_constructor.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4036ff20455bbb9d42e9b570d4467844ca145655
--- /dev/null
+++ b/sidechain/src/scwrl3_rotamer_constructor.hh
@@ -0,0 +1,59 @@
+// Copyright (c) 2013-2018, SIB - Swiss Institute of Bioinformatics and
+// Biozentrum - University of Basel
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+
+#ifndef PROMOD3_SCWRL3_ROTAMER_CONSTRUCTOR_HH
+#define PROMOD3_SCWRL3_ROTAMER_CONSTRUCTOR_HH
+
+#include <promod3/sidechain/rotamer_constructor.hh>
+
+namespace promod3 { namespace sidechain {
+
+class SCWRL3RotamerConstructor;
+typedef boost::shared_ptr<SCWRL3RotamerConstructor> SCWRL3RotamerConstructorPtr;
+
+class SCWRL3RotamerConstructor : public RotamerConstructor{
+
+public:
+
+ SCWRL3RotamerConstructor(bool cb_in_sidechain = false);
+
+ virtual ~SCWRL3RotamerConstructor() { }
+
+ // Assign internal energies to rotamer groups
+ virtual void AssignInternalEnergies(RRMRotamerGroupPtr group,
+ RotamerID id,
+ uint residue_index,
+ Real phi = -1.0472,
+ Real psi = -0.7854,
+ bool n_ter = false,
+ bool c_ter = false);
+
+ virtual void AssignInternalEnergies(FRMRotamerGroupPtr group,
+ RotamerID id,
+ uint residue_index,
+ Real phi = -1.0472,
+ Real psi = -0.7854,
+ bool n_ter = false,
+ bool c_ter = false);
+
+private:
+
+ void ParametrizeParticle(int atom_idx, bool is_hydrogen, Particle& p);
+};
+
+}} // ns
+
+#endif