From 1055f2713bdd34bd5c5f722d8517c1752d830cac Mon Sep 17 00:00:00 2001 From: Gerardo Tauriello <gerardo.tauriello@unibas.ch> Date: Mon, 18 Sep 2017 19:00:28 +0200 Subject: [PATCH] SCHWED-1740: test multi-loop in MmSystemCreator --- loop/tests/test_mm_system_creator.cc | 149 +++++++++++++++++++++++++-- 1 file changed, 142 insertions(+), 7 deletions(-) diff --git a/loop/tests/test_mm_system_creator.cc b/loop/tests/test_mm_system_creator.cc index 0ba27a2c..48f14fcf 100644 --- a/loop/tests/test_mm_system_creator.cc +++ b/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(); -- GitLab