From ee8e3193b9647774495fd846c7a8f7696ee46fb3 Mon Sep 17 00:00:00 2001
From: Gerardo Tauriello <gerardo.tauriello@unibas.ch>
Date: Tue, 19 Sep 2017 18:23:36 +0200
Subject: [PATCH] SCHWED-1740: handle loops of length 1 in MM system creator
---
loop/src/mm_system_creator.cc | 10 +-
loop/src/mm_system_creator.hh | 19 ++--
loop/tests/test_mm_system_creator.cc | 133 +++++++++++++++++++++++----
3 files changed, 136 insertions(+), 26 deletions(-)
diff --git a/loop/src/mm_system_creator.cc b/loop/src/mm_system_creator.cc
index 5ffe7a97..2395048d 100644
--- a/loop/src/mm_system_creator.cc
+++ b/loop/src/mm_system_creator.cc
@@ -557,8 +557,8 @@ MmSystemCreator::SetupTopology_(const DisulfidBridgeVector& disulfid_bridges) {
masses = ff_lookup_->GetMasses(ff_aa, is_nter, is_cter);
std::copy(masses.begin(), masses.end(), &atom_masses[first_idx]);
// loop check
- const IndexLocation_ i_loc = GetIndexLocation_(i_res);
- if (i_loc == ON_N_STEM && !is_nter) {
+ const uint i_loc = GetIndexLocation_(i_res);
+ if ((i_loc & ON_N_STEM) && !is_nter) {
// fix N-stem (only fix N, CA, CB)
const uint idx_N = ff_lookup_->GetHeavyIndex(ff_aa, BB_N_INDEX);
atom_masses[first_idx + idx_N] = 0;
@@ -568,7 +568,8 @@ MmSystemCreator::SetupTopology_(const DisulfidBridgeVector& disulfid_bridges) {
const uint idx_CB = ff_lookup_->GetHeavyIndex(ff_aa, BB_CB_INDEX);
atom_masses[first_idx + idx_CB] = 0;
}
- } else if (i_loc == ON_C_STEM && !is_cter) {
+ }
+ if ((i_loc & ON_C_STEM) && !is_cter) {
// fix C-stem (only fix CA, CB, C, O)
const uint idx_CA = ff_lookup_->GetHeavyIndex(ff_aa, BB_CA_INDEX);
atom_masses[first_idx + idx_CA] = 0;
@@ -580,7 +581,8 @@ MmSystemCreator::SetupTopology_(const DisulfidBridgeVector& disulfid_bridges) {
const uint idx_CB = ff_lookup_->GetHeavyIndex(ff_aa, BB_CB_INDEX);
atom_masses[first_idx + idx_CB] = 0;
}
- } else if (i_loc == OUT_OF_LOOP) {
+ }
+ if (i_loc == OUT_OF_LOOP) {
// fix surrounding (either all or only heavy atoms)
if (fix_surrounding_hydrogens_) {
std::fill_n(&atom_masses[first_idx], masses.size(), Real(0));
diff --git a/loop/src/mm_system_creator.hh b/loop/src/mm_system_creator.hh
index 9ee72384..2602a27e 100644
--- a/loop/src/mm_system_creator.hh
+++ b/loop/src/mm_system_creator.hh
@@ -112,16 +112,23 @@ protected:
// dealing with loops
enum IndexLocation_ {
- ON_N_STEM, ON_C_STEM, IN_LOOP, OUT_OF_LOOP
+ // to be or-ed
+ OUT_OF_LOOP = 0,
+ ON_N_STEM = 1,
+ ON_C_STEM = 2,
+ IN_LOOP = 4
};
- IndexLocation_ GetIndexLocation_(uint i_res) {
- // check all loops
+ uint GetIndexLocation_(uint i_res) {
+ // check all loops (no loops of length 0 possible!)
const uint num_loops = loop_start_indices_.size();
for (uint i_loop = 0; i_loop < num_loops; ++i_loop) {
const uint tst = i_res - loop_start_indices_[i_loop];
- if (tst == 0) return ON_N_STEM;
- if (tst == loop_lengths_[i_loop]-1) return ON_C_STEM;
- if (tst < loop_lengths_[i_loop]-1) return IN_LOOP;
+ if (tst < loop_lengths_[i_loop]) {
+ uint res = IN_LOOP;
+ if (tst == 0) res |= ON_N_STEM;
+ if (tst == loop_lengths_[i_loop]-1) res |= ON_C_STEM;
+ return res;
+ }
}
return OUT_OF_LOOP;
}
diff --git a/loop/tests/test_mm_system_creator.cc b/loop/tests/test_mm_system_creator.cc
index 48f14fcf..743cac93 100644
--- a/loop/tests/test_mm_system_creator.cc
+++ b/loop/tests/test_mm_system_creator.cc
@@ -365,6 +365,36 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
// check result
AllAtomPositionsPtr new_pos = all_pos.Copy();
mm_sim.ExtractLoopPositions(*new_pos, res_indices);
+ const uint num_residues = all_pos.GetNumResidues();
+ for (uint res_idx = 0; res_idx < num_residues; ++res_idx) {
+ // expect stems out-of-loop and parts of stems to be fixed and rest moving
+ const bool in_loop = (res_idx >= start_res_idx && res_idx <= end_res_idx);
+ const bool is_nstem = (start_res_idx > 0 && res_idx == start_res_idx);
+ const bool is_cstem = ( end_res_idx < num_residues - 1
+ && res_idx == end_res_idx);
+ // check atoms
+ const uint first_idx = all_pos.GetFirstIndex(res_idx);
+ const uint last_idx = all_pos.GetLastIndex(res_idx);
+ for (uint i = first_idx; i <= last_idx; ++i) {
+ const uint atom_idx = i - first_idx;
+ const Real dist = geom::Distance(all_pos.GetPos(i), new_pos->GetPos(i));
+ if (!in_loop) {
+ BOOST_CHECK_EQUAL(dist, Real(0));
+ } else if (is_nstem && ( atom_idx == BB_N_INDEX
+ || atom_idx == BB_CA_INDEX
+ || atom_idx == BB_CB_INDEX)) {
+ BOOST_CHECK(dist < Real(1e-6));
+ } else if (is_cstem && ( atom_idx == BB_CA_INDEX
+ || atom_idx == BB_CB_INDEX
+ || atom_idx == BB_C_INDEX
+ || atom_idx == BB_O_INDEX)) {
+ BOOST_CHECK(dist < Real(1e-6));
+ } else {
+ BOOST_CHECK(dist > 1e-6);
+ }
+ }
+ }
+ // TEST
uint num_moves = 0;
for (uint i = 0; i < all_pos.GetNumAtoms(); ++i) {
const Real dist = geom::Distance(all_pos.GetPos(i), new_pos->GetPos(i));
@@ -409,7 +439,9 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
void CheckLoopMinimize(const AllAtomPositions& all_pos,
bool fix_surrounding_hydrogens,
bool kill_es, Real cutoff) {
+ ////////////////////////////////////////////////////////////////////////////
// run single loops
+ ////////////////////////////////////////////////////////////////////////////
const uint num_residues = all_pos.GetNumResidues();
CheckLoopMinimize(all_pos, fix_surrounding_hydrogens, kill_es, cutoff, 0,
num_residues-1);
@@ -418,7 +450,17 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
CheckLoopMinimize(all_pos, fix_surrounding_hydrogens, kill_es, cutoff,
num_residues-5, num_residues-1);
+ ////////////////////////////////////////////////////////////////////////////
// set up multi-loop problem (follows CheckLoopMinimize/SetupLoopShift)
+ ////////////////////////////////////////////////////////////////////////////
+ // fix problem
+ const uint start_res_idx1 = 10;
+ const uint end_res_idx1 = 15;
+ const uint loop_length1 = end_res_idx1 - start_res_idx1 + 1;
+ const uint start_res_idx2 = num_residues-5;
+ const uint end_res_idx2 = num_residues-1;
+ const uint loop_length2 = end_res_idx2 - start_res_idx2 + 1;
+ // setup mm_sim
ForcefieldLookupPtr ff_lookup = ForcefieldLookup::GetDefault();
MmSystemCreator mm_sim(ff_lookup, fix_surrounding_hydrogens, kill_es, cutoff);
MmSystemCreator mm_sim_inaccurate(ff_lookup, fix_surrounding_hydrogens,
@@ -430,10 +472,10 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
// same loops as above
std::vector<uint> loop_start_indices;
std::vector<uint> loop_lengths;
- loop_start_indices.push_back(10);
- loop_lengths.push_back(6);
- loop_start_indices.push_back(num_residues-5);
- loop_lengths.push_back(5);
+ loop_start_indices.push_back(start_res_idx1);
+ loop_lengths.push_back(loop_length1);
+ loop_start_indices.push_back(start_res_idx2);
+ loop_lengths.push_back(loop_length2);
// get terminal info
std::vector<bool> is_n_ter(num_residues, false);
std::vector<bool> is_c_ter(num_residues, false);
@@ -458,20 +500,38 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
// check result
AllAtomPositionsPtr new_pos = all_pos.Copy();
mm_sim.ExtractLoopPositions(*new_pos, res_indices);
- uint num_moves = 0;
- const uint first_idx1 = all_pos.GetFirstIndex(10);
- const uint last_idx1 = all_pos.GetLastIndex(15);
- const uint first_idx2 = all_pos.GetFirstIndex(num_residues-5);
- const uint last_idx2 = all_pos.GetLastIndex(num_residues-1);
- for (uint i = 0; i < all_pos.GetNumAtoms(); ++i) {
- const Real dist = geom::Distance(all_pos.GetPos(i), new_pos->GetPos(i));
- if ((i < first_idx1 || i > last_idx1) && (i < first_idx2 || i > last_idx2)) {
- BOOST_CHECK_EQUAL(dist, Real(0));
- } else if (dist > 0) {
- ++num_moves;
+ for (uint res_idx = 0; res_idx < num_residues; ++res_idx) {
+ // expect stems out-of-loop and parts of stems to be fixed and rest moving
+ const bool in_loop = ( ( res_idx >= start_res_idx1
+ && res_idx <= end_res_idx1)
+ || (res_idx >= start_res_idx2
+ && res_idx <= end_res_idx2));
+ const bool is_nstem = ( res_idx == start_res_idx1
+ || res_idx == start_res_idx2);
+ // HACK: end_res_idx2 is C-terminal!
+ const bool is_cstem = (res_idx == end_res_idx1);
+ // check atoms
+ const uint first_idx = all_pos.GetFirstIndex(res_idx);
+ const uint last_idx = all_pos.GetLastIndex(res_idx);
+ for (uint i = first_idx; i <= last_idx; ++i) {
+ const uint atom_idx = i - first_idx;
+ const Real dist = geom::Distance(all_pos.GetPos(i), new_pos->GetPos(i));
+ if (!in_loop) {
+ BOOST_CHECK_EQUAL(dist, Real(0));
+ } else if (is_nstem && ( atom_idx == BB_N_INDEX
+ || atom_idx == BB_CA_INDEX
+ || atom_idx == BB_CB_INDEX)) {
+ BOOST_CHECK(dist < Real(1e-6));
+ } else if (is_cstem && ( atom_idx == BB_CA_INDEX
+ || atom_idx == BB_CB_INDEX
+ || atom_idx == BB_C_INDEX
+ || atom_idx == BB_O_INDEX)) {
+ BOOST_CHECK(dist < Real(1e-6));
+ } else {
+ BOOST_CHECK(dist > 1e-6);
+ }
}
}
- BOOST_CHECK(num_moves > 4*mm_sim.GetNumLoopResidues());
// repeat with inaccurate one
mm_sim_inaccurate.SetupSystem(all_pos, res_indices, loop_start_indices,
loop_lengths, is_n_ter, is_c_ter,
@@ -485,6 +545,47 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
const Real dist = geom::Distance(new_pos->GetPos(i), new_pos2->GetPos(i));
BOOST_CHECK_EQUAL(dist, Real(0));
}
+
+ ////////////////////////////////////////////////////////////////////////////
+ // try loop of length 1 (only sc should move!)////////////////////////////////////////////////////////////////////////////
+ // find interesting case
+ uint loop_idx = 1;
+ for (; loop_idx < num_residues-1; ++loop_idx) {
+ const uint first_idx = all_pos.GetFirstIndex(loop_idx);
+ const uint last_idx = all_pos.GetLastIndex(loop_idx);
+ if (last_idx - first_idx > 4) break;
+ }
+ // build system
+ loop_start_indices.assign(1, loop_idx);
+ loop_lengths.assign(1, 1);
+ mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices, loop_lengths,
+ is_n_ter, is_c_ter, disulfid_bridges);
+ // run it for a bit
+ sim = mm_sim.GetSimulation();
+ pot_ref_loop = sim->GetPotentialEnergy();
+ sim->ApplySD(0.01, 25);
+ BOOST_CHECK(sim->GetPotentialEnergy() < pot_ref_loop);
+ // check result
+ new_pos = all_pos.Copy();
+ mm_sim.ExtractLoopPositions(*new_pos, res_indices);
+ const uint first_idx = all_pos.GetFirstIndex(loop_idx);
+ const uint last_idx = all_pos.GetLastIndex(loop_idx);
+ for (uint i = 0; i < all_pos.GetNumAtoms(); ++i) {
+ const Real dist = geom::Distance(all_pos.GetPos(i), new_pos->GetPos(i));
+ if (i < first_idx || i > last_idx) {
+ BOOST_CHECK_EQUAL(dist, Real(0));
+ } else {
+ // expect backbone to be fixed
+ const uint atom_idx = i - first_idx;
+ if ( atom_idx == BB_N_INDEX || atom_idx == BB_CA_INDEX
+ || atom_idx == BB_C_INDEX || atom_idx == BB_O_INDEX
+ || atom_idx == BB_CB_INDEX) {
+ BOOST_CHECK(dist < Real(1e-6));
+ } else {
+ BOOST_CHECK(dist > 1e-6);
+ }
+ }
+ }
}
} // anon ns
--
GitLab