SCHWED-1740: handle loops of length 1 in MM system creator

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));

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;
   }

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