SCHWED-1740: test multi-loop in MmSystemCreator

loop/tests/test_mm_system_creator.cc

@@ -275,11 +275,34 @@ void SetupLoop(const AllAtomPositions& all_pos, uint start_res_idx,
   // get disulfid bridges
   std::vector< std::pair<uint,uint> > disulfid_bridges;
   disulfid_bridges = mm_sim.GetDisulfidBridges(all_pos, res_indices);
-  // build system and check it
+  // build system
   mm_sim.SetupSystem(all_pos, res_indices, loop_length, is_n_ter, is_c_ter,
                      disulfid_bridges);
 }
 
+// same as above but passing loop_start_indices and loop_lengths
+void SetupLoopShift(const AllAtomPositions& all_pos, uint start_res_idx,
+                    uint end_res_idx, MmSystemCreator& mm_sim) {
+  // setup
+  const uint num_residues = all_pos.GetNumResidues();
+  std::vector<uint> loop_start_indices(1, start_res_idx);
+  std::vector<uint> loop_lengths(1, end_res_idx - start_res_idx + 1);
+  // get indices
+  std::vector<uint> res_indices;
+  for (uint i = 0; i < num_residues; ++i) res_indices.push_back(i);
+  // get terminal info
+  std::vector<bool> is_n_ter(num_residues, false);
+  std::vector<bool> is_c_ter(num_residues, false);
+  is_n_ter[0] = true;
+  is_c_ter[num_residues-1] = true;
+  // get disulfid bridges
+  std::vector< std::pair<uint,uint> > disulfid_bridges;
+  disulfid_bridges = mm_sim.GetDisulfidBridges(all_pos, res_indices);
+  // build system
+  mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices, loop_lengths,
+                     is_n_ter, is_c_ter, disulfid_bridges);
+}
+
 void CheckLoops(const AllAtomPositions& all_pos, bool fix_surrounding_hydrogens,
                 bool kill_es, Real cutoff) {
   // generate MM sim
@@ -333,23 +356,26 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
   for (uint i = start_res_idx; i <= end_res_idx; ++i) res_indices.push_back(i);
   const uint first_idx = all_pos.GetFirstIndex(start_res_idx);
   const uint last_idx = all_pos.GetLastIndex(end_res_idx);
-  // run full chain for a bit
+  // run it for a bit
   SetupLoop(all_pos, start_res_idx, end_res_idx, mm_sim);
   ost::mol::mm::SimulationPtr sim = mm_sim.GetSimulation();
+  Real pot_ref_loop = sim->GetPotentialEnergy();
   sim->ApplySD(0.01, 25);
+  BOOST_CHECK(sim->GetPotentialEnergy() < pot_ref_loop);
   // check result
   AllAtomPositionsPtr new_pos = all_pos.Copy();
   mm_sim.ExtractLoopPositions(*new_pos, res_indices);
   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));
     if (i < first_idx || i > last_idx) {
-      const Real dist = geom::Distance(all_pos.GetPos(i), new_pos->GetPos(i));
-      BOOST_CHECK(dist == 0);
-    } else {
+      BOOST_CHECK_EQUAL(dist, Real(0));
+    } else if (dist > 0) {
       ++num_moves;
     }
   }
   BOOST_CHECK(num_moves > 4*loop_length);
+  
   // repeat with inaccurate one
   SetupLoop(all_pos, start_res_idx, end_res_idx, mm_sim_inaccurate);
   sim = mm_sim_inaccurate.GetSimulation();
@@ -359,14 +385,31 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
   mm_sim_inaccurate.ExtractLoopPositions(*new_pos2, res_indices);
   for (uint i = 0; i < all_pos.GetNumAtoms(); ++i) {
     const Real dist = geom::Distance(new_pos->GetPos(i), new_pos2->GetPos(i));
-    BOOST_CHECK(dist == 0);
+    BOOST_CHECK_EQUAL(dist, Real(0));
+  }
+
+  // repeat without res_indices reshuffling (boring if idx = 0)
+  if (start_res_idx > 0) {
+    SetupLoopShift(all_pos, start_res_idx, end_res_idx, mm_sim);
+    sim = mm_sim.GetSimulation();
+    BOOST_CHECK_EQUAL(sim->GetPotentialEnergy(), pot_ref_loop);
+    sim->ApplySD(0.01, 25);
+    // check that result is identical
+    new_pos2 = all_pos.Copy();
+    std::vector<uint> res_indices2;
+    for (uint i = 0; i <= end_res_idx; ++i) res_indices2.push_back(i);
+    mm_sim.ExtractLoopPositions(*new_pos2, res_indices2);
+    for (uint i = 0; i < all_pos.GetNumAtoms(); ++i) {
+      const Real dist = geom::Distance(new_pos->GetPos(i), new_pos2->GetPos(i));
+      BOOST_CHECK_EQUAL(dist, Real(0));
+    }
   }
 }
 
 void CheckLoopMinimize(const AllAtomPositions& all_pos,
                        bool fix_surrounding_hydrogens,
                        bool kill_es, Real cutoff) {
-  // run all
+  // run single loops
   const uint num_residues = all_pos.GetNumResidues();
   CheckLoopMinimize(all_pos, fix_surrounding_hydrogens, kill_es, cutoff, 0,
                     num_residues-1);
@@ -374,6 +417,74 @@ void CheckLoopMinimize(const AllAtomPositions& all_pos,
                     15);
   CheckLoopMinimize(all_pos, fix_surrounding_hydrogens, kill_es, cutoff,
                     num_residues-5, num_residues-1);
+  
+  // set up multi-loop problem (follows CheckLoopMinimize/SetupLoopShift)
+  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,
+                                    kill_es, cutoff, true);
+  mm_sim.SetCpuPlatformSupport(false);
+  mm_sim_inaccurate.SetCpuPlatformSupport(false);
+  std::vector<uint> res_indices;
+  for (uint i = 0; i < num_residues; ++i) res_indices.push_back(i);
+  // 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);
+  // get terminal info
+  std::vector<bool> is_n_ter(num_residues, false);
+  std::vector<bool> is_c_ter(num_residues, false);
+  is_n_ter[0] = true;
+  is_c_ter[num_residues-1] = true;
+  // get disulfid bridges
+  std::vector< std::pair<uint,uint> > disulfid_bridges;
+  disulfid_bridges = mm_sim.GetDisulfidBridges(all_pos, res_indices);
+  // build system
+  mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices, loop_lengths,
+                     is_n_ter, is_c_ter, disulfid_bridges);
+  // check it
+  BOOST_CHECK_EQUAL(mm_sim.GetNumResidues(), num_residues);
+  BOOST_CHECK_EQUAL(mm_sim.GetNumLoopResidues(), 6u + 5u);
+  BOOST_CHECK(mm_sim.GetLoopStartIndices() == loop_start_indices);
+  BOOST_CHECK(mm_sim.GetLoopLengths() == loop_lengths);
+  // run it for a bit
+  ost::mol::mm::SimulationPtr sim = mm_sim.GetSimulation();
+  Real pot_ref_loop = sim->GetPotentialEnergy();
+  sim->ApplySD(0.01, 25);
+  BOOST_CHECK(sim->GetPotentialEnergy() < pot_ref_loop);
+  // 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;
+    }
+  }
+  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,
+                                disulfid_bridges);
+  sim = mm_sim_inaccurate.GetSimulation();
+  sim->ApplySD(0.01, 25);
+  // check that result is identical
+  AllAtomPositionsPtr new_pos2 = all_pos.Copy();
+  mm_sim_inaccurate.ExtractLoopPositions(*new_pos2, res_indices);
+  for (uint i = 0; i < all_pos.GetNumAtoms(); ++i) {
+    const Real dist = geom::Distance(new_pos->GetPos(i), new_pos2->GetPos(i));
+    BOOST_CHECK_EQUAL(dist, Real(0));
+  }
 }
 
 } // anon ns
@@ -533,6 +644,30 @@ BOOST_AUTO_TEST_CASE(test_mm_system_creator_loop_exceptions) {
   BOOST_CHECK_THROW(mm_sim.ExtractLoopPositions(*loop_pos), promod3::Error);
   loop_pos = all_pos.Extract(res_indices_short);
   BOOST_CHECK_THROW(mm_sim.ExtractLoopPositions(*loop_pos), promod3::Error);
+
+  // check/break multi-loop
+  std::vector<uint> loop_start_indices;
+  std::vector<uint> loop_lengths;
+  BOOST_CHECK_THROW(mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices,
+                                       loop_lengths, is_n_ter, is_c_ter,
+                                       dis_bridges),
+                    promod3::Error);
+  loop_start_indices.push_back(1);
+  loop_lengths.push_back(5);
+  BOOST_CHECK_THROW(mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices,
+                                       loop_lengths, is_n_ter, is_c_ter,
+                                       dis_bridges),
+                    promod3::Error);
+  loop_lengths[0] = 0;
+  BOOST_CHECK_THROW(mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices,
+                                       loop_lengths, is_n_ter, is_c_ter,
+                                       dis_bridges),
+                    promod3::Error);
+  loop_start_indices.push_back(1);
+  BOOST_CHECK_THROW(mm_sim.SetupSystem(all_pos, res_indices, loop_start_indices,
+                                       loop_lengths, is_n_ter, is_c_ter,
+                                       dis_bridges),
+                    promod3::Error);
 }
 
 BOOST_AUTO_TEST_SUITE_END();