From 31d4bd9f9934506d230d5ba3a221d7e97ee4380a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Conrad=20H=C3=BCbler?= Date: Fri, 19 Apr 2024 14:43:58 +0200 Subject: [PATCH] save changes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Conrad Hübler --- .../optimiser/OptimiseDipoleScaling.h | 19 +- src/capabilities/rmsd.cpp | 381 ++++++++++++------ src/capabilities/rmsd.h | 99 +++-- src/core/molecule.cpp | 13 + 4 files changed, 361 insertions(+), 151 deletions(-) diff --git a/src/capabilities/optimiser/OptimiseDipoleScaling.h b/src/capabilities/optimiser/OptimiseDipoleScaling.h index d7929bd..9e9f2c9 100644 --- a/src/capabilities/optimiser/OptimiseDipoleScaling.h +++ b/src/capabilities/optimiser/OptimiseDipoleScaling.h @@ -113,16 +113,25 @@ inline std::pair> OptimiseScaling(const Molecule* mo // std::cout << parameter.transpose() << std::endl; double fx; int converged = solver.InitializeSingleSteps(fun, parameter, fx); - - for (iteration = 1; iteration <= maxiter && fx > threshold; ++iteration) { - solver.SingleStep(fun, parameter, fx); + bool loop = true; + for (iteration = 1; iteration <= maxiter && fx > threshold && loop; ++iteration) { + try { + solver.SingleStep(fun, parameter, fx); + + } catch (const std::logic_error& error_result) { + loop = false; + } catch (const std::runtime_error& error_result) { + loop = false; + } std::cout << "Iteration: " << iteration << " Difference: " << fx << std::endl; } std::vector scaling(parameter.size()); - for (int i = 0; i < scaling.size(); ++i) + for (int i = 0; i < scaling.size(); ++i) { scaling[i] = parameter(i); + std::cout << scaling[i] << " "; + } auto dipole_vector = molecule->CalculateDipoleMoment(scaling); dipole = sqrt(dipole_vector[0] * dipole_vector[0] + dipole_vector[1] * dipole_vector[1] + dipole_vector[2] * dipole_vector[2]); - + std::cout << std::endl; return std::pair>(dipole, scaling); } diff --git a/src/capabilities/rmsd.cpp b/src/capabilities/rmsd.cpp index be77d9f..3599cd8 100644 --- a/src/capabilities/rmsd.cpp +++ b/src/capabilities/rmsd.cpp @@ -16,7 +16,6 @@ * along with this program. If not, see . * */ - #include "rmsd_functions.h" #include "munkres.h" @@ -30,6 +29,7 @@ #include "external/CxxThreadPool/include/CxxThreadPool.h" #include +#include #include #include #include @@ -48,6 +48,8 @@ #include #include +#include + #include "json.hpp" using json = nlohmann::json; @@ -172,6 +174,7 @@ void RMSDDriver::LoadControlJson() m_split = Json2KeyWord(m_defaults, "split"); m_nofree = Json2KeyWord(m_defaults, "nofree"); m_dmix = Json2KeyWord(m_defaults, "dmix"); + m_maxtrial = Json2KeyWord(m_defaults, "maxtrial"); #pragma message("these hacks to overcome the json stuff are not nice, TODO!") try { std::string element = m_defaults["Element"].get(); @@ -216,6 +219,7 @@ void RMSDDriver::LoadControlJson() } else m_method = 1; + m_costmatrix = Json2KeyWord(m_defaults, "costmatrix"); std::string order = Json2KeyWord(m_defaults, "order"); int cycles = Json2KeyWord(m_defaults, "cycles"); if (cycles != -1) @@ -261,8 +265,8 @@ void RMSDDriver::start() m_reference_centered = m_reference; m_target_centered = m_target; - m_reference_centered.Center(); - m_target_centered.Center(); + m_reference.Center(); + m_target.Center(); m_target_aligned = m_target; m_reference_aligned.LoadMolecule(m_reference); @@ -273,6 +277,7 @@ void RMSDDriver::start() m_target_reordered = ApplyOrder(m_reorder_rules, m_target); m_target = m_target_reordered; m_rmsd = m_results.begin()->first; + // std::cout << m_rmsd << std::endl; rmsd_calculated = true; } } @@ -564,10 +569,11 @@ StructComp RMSDDriver::Rule2RMSD(const std::vector rules, int fragment) void RMSDDriver::clear() { m_results.clear(); - m_connectivity.clear(); m_stored_rules.clear(); m_prepared_cost_matrices.clear(); m_intermediate_cost_matrices.clear(); + m_intermedia_rules.clear(); + m_rmsd = 0.0; } @@ -672,9 +678,7 @@ std::vector RMSDDriver::IndivRMSD(const Molecule& reference_mol, const M Geometry reference = CenterMolecule(reference_mol, m_fragment_reference); Geometry target = CenterMolecule(target_mol, m_fragment_target); - Eigen::MatrixXd tar = target.transpose(); - - Geometry rotated = tar.transpose() * R; + Geometry rotated = target * R; for (int i = 0; i < rotated.rows(); ++i) { terms.push_back((rotated(i, 0) - reference(i, 0)) * (rotated(i, 0) - reference(i, 0)) + (rotated(i, 1) - reference(i, 1)) * (rotated(i, 1) - reference(i, 1)) + (rotated(i, 2) - reference(i, 2)) * (rotated(i, 2) - reference(i, 2))); @@ -763,9 +767,14 @@ std::pair RMSDDriver::InitialisePair() void RMSDDriver::ReorderMolecule() { - if (!m_nofree) - TemplateFree(); + auto R = GetOperateVectors(m_reference, m_target); + Eigen::Matrix3d rotation = R.first; + + auto blob = MakeCostMatrix(OptimiseRotation(rotation)); + if (m_method == 4) + blob.first = 1; + InsertRotation(blob); if (m_method == 1) ReorderIncremental(); else if (m_method == 2) @@ -774,6 +783,8 @@ void RMSDDriver::ReorderMolecule() HeavyTemplate(); else if (m_method == 4) AtomTemplate(); + else if (m_method == 5) + TemplateFree(); else if (m_method == 6) { if (!MolAlignLib()) { return; @@ -787,6 +798,40 @@ void RMSDDriver::ReorderMolecule() m_target_aligned = m_target; } +void RMSDDriver::FinaliseTemplate() +{ + Molecule target = m_target; + std::vector> rules = m_stored_rules; + double rmsd_prev = 10; + int eq_counter = 0, incr_counter = 0; + int iter = 0; + RunTimer time; + double prev_cost = -1; + for (auto permutation : m_prepared_cost_matrices) { + iter++; + if (eq_counter > m_maxtrial || incr_counter > m_maxtrial || iter > m_maxtrial) + break; + if (prev_cost != -1 && prev_cost < permutation.first / 10) { + std::cout << permutation.first << std::endl; + break; + } + prev_cost = permutation.first; + auto result = SolveCostMatrix(permutation.second); + m_target_reordered = ApplyOrder(result, target); + double rmsd = Rules2RMSD(result); + if (!m_silent) + std::cout << permutation.first << " " << rmsd << " " << eq_counter << " " << time.Elapsed() << std::endl; + m_results.insert(std::pair>(rmsd, result)); + time.Reset(); + + eq_counter += std::abs(rmsd - rmsd_prev) < 1e-3; + incr_counter += rmsd > rmsd_prev; + rmsd_prev = rmsd; + } + + m_reorder_rules = m_results.begin()->second; +} + void RMSDDriver::DistanceTemplate() { if (!m_silent) @@ -811,50 +856,158 @@ void RMSDDriver::HeavyTemplate() PrepareHeavyTemplate(); } +Matrix RMSDDriver::OptimiseRotation(const Eigen::Matrix3d& rotation) +{ + Eigen::Vector3d euler = rotation.eulerAngles(0, 1, 2); + + LBFGSRotation opt(3); + opt.m_reference = m_reference.getGeometry(); + opt.m_target = m_target.getGeometry(); + opt.m_reference_atoms = m_reference.Atoms(); + opt.m_target_atoms = m_target.Atoms(); + opt.m_costmatrix = m_costmatrix; + + int iteration; + int maxiter = 1000; + LBFGSParam param; + LBFGSSolver solver(param); + Vector vector(3); + for (int i = 0; i < 3; ++i) + vector(i) = euler(i); + + double fx; + int converged = solver.InitializeSingleSteps(opt, vector, fx); + Vector v_old = vector; + for (iteration = 1; iteration <= maxiter; ++iteration) { + try { + solver.SingleStep(opt, vector, fx); + } catch (const std::logic_error& error_result) { + break; + } catch (const std::runtime_error& error_result) { + break; + } + + if ((vector - v_old).norm() < 1e-5) + break; + v_old = vector; + } + Eigen::Matrix3d n; + n = Eigen::AngleAxisd(vector[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxisd(vector[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxisd(vector[2], Eigen::Vector3d::UnitZ()); + + return n; +} void RMSDDriver::TemplateFree() { - /* Eigen::Matrix3d indentity = Eigen::MatrixXd::Zero(3, 3); - indentity(0, 0) = -1; - indentity(1, 1) = -1; - indentity(2, 2) = -1; */ auto R = GetOperateVectors(m_reference, m_target); - auto blob = MakeCostMatrix(R.first); + Eigen::Matrix3d rotation = R.first; - m_cost_limit = blob.first; - if (!m_silent) - std::cout << blob.first << std::endl; + // Eigen::Matrix3d n = OptimiseRotation(R.first); + for (int i = 0; i < 360; i += 60) { + Eigen::Matrix3d m = Eigen::Matrix3d::Zero(); + m(0, 0) = 1; + m(1, 1) = 1; + m(2, 2) = 1; + m = m * Eigen::AngleAxisd(i, Eigen::Vector3d::UnitX()); + // std::cout << m << std::endl; - if (m_method == 4) - blob.first = 0; - m_prepared_cost_matrices.insert(blob); + auto blob = MakeCostMatrix(OptimiseRotation(m)); + // std::cout << blob.first << " " << std::endl; + m_cost_limit = blob.first; + + if (m_method == 4) + blob.first = 1; + InsertRotation(blob); + } + for (int j = 0; j < 360; j += 60) { + Eigen::Matrix3d m = Eigen::Matrix3d::Zero(); + m(0, 0) = 1; + m(1, 1) = 1; + m(2, 2) = 1; + m = m * Eigen::AngleAxisd(j, Eigen::Vector3d::UnitY()); + // std::cout << m << std::endl; + + auto blob = MakeCostMatrix(OptimiseRotation(m)); + // std::cout << blob.first << " " << std::endl ; + + m_cost_limit = blob.first; + + if (m_method == 4) + blob.first = 1; + InsertRotation(blob); + } + for (int k = 0; k < 360; k += 60) { + Eigen::Matrix3d m = Eigen::Matrix3d::Zero(); + m(0, 0) = 1; + m(1, 1) = 1; + m(2, 2) = 1; + m = m * Eigen::AngleAxisd(k, Eigen::Vector3d::UnitZ()); + // std::cout << m << std::endl; + + auto blob = MakeCostMatrix(OptimiseRotation(m)); + // std::cout << blob.first << " " << std::endl ; + + m_cost_limit = blob.first; + + if (m_method == 4) + blob.first = 1; + InsertRotation(blob); + } } -void RMSDDriver::FinaliseTemplate() +void RMSDDriver::InsertRotation(std::pair& rotation) { - Molecule target = m_target; - std::vector> rules = m_stored_rules; - double rmsd_prev = 10; - int eq_counter = 0, incr_counter = 0; - int iter = 0; - RunTimer time; - for (auto permutation : m_prepared_cost_matrices) { - iter++; - auto result = SolveCostMatrix(permutation.second); - m_target_reordered = ApplyOrder(result, target); - double rmsd = Rules2RMSD(result); - if (!m_silent) - std::cout << permutation.first << " " << rmsd << " " << eq_counter << " " << time.Elapsed() << std::endl; - m_results.insert(std::pair>(rmsd, result)); - time.Reset(); - if (eq_counter > 5 || incr_counter > 5 || iter > 5) - break; - eq_counter += std::abs(rmsd - rmsd_prev) < 1e-3; - incr_counter += rmsd > rmsd_prev; - rmsd_prev = rmsd; + for (auto i : m_prepared_cost_matrices) { + // std::cout << (rotation.second - i.second).cwiseAbs().sum()/double(i.second.rows()*i.second.rows()) << std::endl; + if ((rotation.second - i.second).cwiseAbs().sum() / double(i.second.rows() * i.second.rows()) < 10) + return; } + m_prepared_cost_matrices.insert(rotation); +} - m_reorder_rules = m_results.begin()->second; +bool RMSDDriver::TemplateReorder() +{ + auto fragments = CheckFragments(); + + if (fragments.first == -1 || fragments.second == -1) + return false; + + m_fragment = -1; + m_fragment_target = -1; + m_fragment_reference = -1; + + Molecule reference_mol = m_reference.getFragmentMolecule(fragments.first); + Molecule target_mol = m_target.getFragmentMolecule(fragments.second); + + auto operators = GetOperateVectors(fragments.first, fragments.second); + Eigen::Matrix3d R = operators.first; + + Geometry cached_reference = m_reference.getGeometryByFragment(fragments.first, m_protons); + Geometry cached_target = m_target.getGeometryByFragment(fragments.second, m_protons); + + Geometry ref = GeometryTools::TranslateMolecule(m_reference, GeometryTools::Centroid(cached_reference), Position{ 0, 0, 0 }); + Geometry tget = GeometryTools::TranslateMolecule(m_target, GeometryTools::Centroid(cached_target), Position{ 0, 0, 0 }); + + // Eigen::MatrixXd tar = tget.transpose(); + + Geometry rotated = tget * R; + + Molecule ref_mol = m_reference; + ref_mol.setGeometry(ref); + + Molecule tar_mol = m_target; + tar_mol.setGeometry(rotated); + + auto blob = MakeCostMatrix(cached_reference, rotated); + m_cost_limit = blob.first; + if (!m_silent) + + if (m_method == 4) + blob.first = 1; + InsertRotation(blob); + return true; } std::pair, std::vector> RMSDDriver::PrepareHeavyTemplate() @@ -893,7 +1046,6 @@ std::pair, std::vector> RMSDDriver::PrepareHeavyTemplate() for (auto index : m_stored_rules[i]) tmp.push_back(target_indicies[index]); transformed_rules.push_back(tmp); - m_intermedia_rules.push_back(tmp); } m_stored_rules = transformed_rules; std::vector target_indices = m_reorder_rules; @@ -943,18 +1095,18 @@ std::pair, std::vector> RMSDDriver::PrepareAtomTemplate(in std::vector tmp; for (auto index : m_stored_rules[i]) tmp.push_back(target_indicies[index]); - m_intermedia_rules.push_back(tmp); + auto OperateVectors = GetOperateVectors(reference_indicies, tmp); auto result = MakeCostMatrix(OperateVectors.first); result.first = m_prepared_cost_matrices.size(); - m_prepared_cost_matrices.insert(result); + InsertRotation(result); } for (const auto& indices : m_intermedia_rules) { auto result = MakeCostMatrix(reference_indicies, indices); result.first = m_prepared_cost_matrices.size(); - m_prepared_cost_matrices.insert(result); + InsertRotation(result); } return std::pair, std::vector>(reference_indicies, reference_indicies); @@ -1003,15 +1155,15 @@ std::pair, std::vector> RMSDDriver::PrepareAtomTemplate(co m_intermedia_rules.push_back(tmp); auto OperateVectors = GetOperateVectors(reference_indicies, tmp); auto result = MakeCostMatrix(OperateVectors.first); - result.first = m_prepared_cost_matrices.size(); + result.first = m_prepared_cost_matrices.size() + 1; - m_prepared_cost_matrices.insert(result); + InsertRotation(result); } for (const auto& indices : m_intermedia_rules) { auto result = MakeCostMatrix(reference_indicies, indices); result.first = m_prepared_cost_matrices.size(); - m_prepared_cost_matrices.insert(result); + InsertRotation(result); } m_stored_rules = transformed_rules; std::vector target_indices = m_reorder_rules; @@ -1024,6 +1176,8 @@ std::pair, std::vector> RMSDDriver::PrepareDistanceTemplat if (!m_silent) std::cout << "Start Prepare Template" << std::endl; RunTimer time; + std::map> m_distance_reference, m_distance_target; + Position ref_centroid = m_reference.Centroid(true); Position tar_centroid = m_target.Centroid(true); @@ -1100,7 +1254,7 @@ std::pair, std::vector> RMSDDriver::PrepareDistanceTemplat m_target = cached_target_mol; for (const auto& indices : m_intermedia_rules) { auto result = MakeCostMatrix(reference_indicies, indices); - m_prepared_cost_matrices.insert(result); + InsertRotation(result); } m_intermedia_rules.clear(); @@ -1113,7 +1267,7 @@ std::pair, std::vector> RMSDDriver::PrepareDistanceTemplat auto OperateVectors = GetOperateVectors(reference_indicies, tmp); auto result = MakeCostMatrix(OperateVectors.first); - m_prepared_cost_matrices.insert(result); + InsertRotation(result); } m_stored_rules.clear(); @@ -1137,9 +1291,9 @@ std::pair RMSDDriver::MakeCostMatrix(const std::vector& ref Geometry cached_reference = m_reference_centered.getGeometry(); Geometry cached_target = m_target_centered.getGeometry(); - Eigen::MatrixXd tar = cached_target.transpose(); + // Eigen::MatrixXd tar = cached_target.transpose(); - Geometry rotated = tar.transpose() * R; + Geometry rotated = cached_target * R; Molecule ref_mol = m_reference; ref_mol.setGeometry(cached_reference); @@ -1154,21 +1308,41 @@ std::pair RMSDDriver::MakeCostMatrix(const std::pair RMSDDriver::MakeCostMatrix(const Matrix& rotation) +{ + // Geometry target = m_target.getGeometry().transpose(); + Geometry rotated = m_target.getGeometry() * rotation; + Geometry reference = m_reference.getGeometry(); + return MakeCostMatrix(reference, rotated); +} + std::pair RMSDDriver::MakeCostMatrix(const Geometry& reference, const Geometry& target /*, const std::vector reference_atoms, const std::vector target_atoms*/) { - double penalty = 100; + double penalty = 1e23; Eigen::MatrixXd distance = Eigen::MatrixXd::Zero(m_reference.AtomCount(), m_reference.AtomCount()); double sum = 0; for (int i = 0; i < m_reference.AtomCount(); ++i) { double min = penalty; int row_count = 0; for (int j = 0; j < m_reference.AtomCount(); ++j) { - double d = (target.row(j) - reference.row(i)).norm() - * (target.row(j) - reference.row(i)).norm(); - double norm = (target.row(j).norm() - reference.row(i).norm()) - * (target.row(j).norm() - reference.row(i).norm()); - distance(i, j) = d // + norm - + penalty * (m_reference.Atom(j).first != m_target.Atom(i).first); + double d = (target.row(j) - reference.row(i)).norm(); + double norm = (target.row(j).norm() - reference.row(i).norm()); + if (m_costmatrix == 2) + distance(i, j) = d; + else if (m_costmatrix == 1) + distance(i, j) = d * d; + else if (m_costmatrix == 3) + distance(i, j) = d + norm; + else if (m_costmatrix == 4) + distance(i, j) = d * d + norm * norm; + else if (m_costmatrix == 5) + distance(i, j) = d * norm; + else if (m_costmatrix == 5) + distance(i, j) = d * d * norm * norm; + else + distance(i, j) = d * d; + // + norm + distance(i, j) += penalty * (m_reference.Atom(i).first != m_target.Atom(j).first); min = std::min(min, distance(i, j)); } sum += min; @@ -1176,31 +1350,39 @@ std::pair RMSDDriver::MakeCostMatrix(const Geometry& reference, return std::pair(sum, distance); } -std::pair RMSDDriver::MakeCostMatrix(const Matrix& rotation) +std::pair RMSDDriver::MakeCostMatrix(const Geometry& reference, const Geometry& target, const std::vector& reference_atoms, const std::vector& target_atoms, int costmatrix) { - Geometry target = m_target.getGeometry().transpose(); - Geometry rotated = target.transpose() * rotation; - Geometry reference = m_reference.getGeometry(); - return MakeCostMatrix(reference, rotated); - /* - double penalty = 100; - Eigen::MatrixXd distance = Eigen::MatrixXd::Zero(m_reference.AtomCount(), m_reference.AtomCount()); + double penalty = 1e23; + Eigen::MatrixXd distance = Eigen::MatrixXd::Zero(reference_atoms.size(), reference_atoms.size()); double sum = 0; - for (int i = 0; i < m_reference.AtomCount(); ++i) { + for (int i = 0; i < reference_atoms.size(); ++i) { double min = penalty; int row_count = 0; - for (int j = 0; j < m_reference.AtomCount(); ++j) { - double d = (rotated.row(j) - reference.row(i)).norm() - * (rotated.row(j) - reference.row(i)).norm(); - double norm = (rotated.row(j).norm() - reference.row(i).norm()) - * (rotated.row(j).norm() - reference.row(i).norm()); - distance(i, j) = d - + penalty * (m_reference.Atom(j).first != m_target.Atom(i).first); - min = std::min(min, d); + for (int j = 0; j < reference_atoms.size(); ++j) { + double d = (target.row(j) - reference.row(i)).norm(); + double norm = (target.row(j).norm() - reference.row(i).norm()); + if (costmatrix == 2) + distance(i, j) = d; + else if (costmatrix == 1) + distance(i, j) = d * d; + else if (costmatrix == 3) + distance(i, j) = d + norm; + else if (costmatrix == 4) + distance(i, j) = d * d + norm * norm; + else if (costmatrix == 5) + distance(i, j) = d * norm; + else if (costmatrix == 5) + distance(i, j) = d * d * norm * norm; + else + distance(i, j) = d * d; + + // + norm + distance(i, j) += penalty * (reference_atoms[i] != target_atoms[j]); + min = std::min(min, distance(i, j)); } sum += min; } - return std::pair(sum, distance);*/ + return std::pair(sum, distance); } std::vector RMSDDriver::SolveCostMatrix(Matrix& distance) @@ -1293,49 +1475,6 @@ std::pair RMSDDriver::GetOperateVectors(const Molecule& refere return std::pair(R, translate); } -bool RMSDDriver::TemplateReorder() -{ - auto fragments = CheckFragments(); - - if (fragments.first == -1 || fragments.second == -1) - return false; - - m_fragment = -1; - m_fragment_target = -1; - m_fragment_reference = -1; - - Molecule reference_mol = m_reference.getFragmentMolecule(fragments.first); - Molecule target_mol = m_target.getFragmentMolecule(fragments.second); - - auto operators = GetOperateVectors(fragments.first, fragments.second); - Eigen::Matrix3d R = operators.first; - - Geometry cached_reference = m_reference.getGeometryByFragment(fragments.first, m_protons); - Geometry cached_target = m_target.getGeometryByFragment(fragments.second, m_protons); - - Geometry ref = GeometryTools::TranslateMolecule(m_reference, GeometryTools::Centroid(cached_reference), Position{ 0, 0, 0 }); - Geometry tget = GeometryTools::TranslateMolecule(m_target, GeometryTools::Centroid(cached_target), Position{ 0, 0, 0 }); - - Eigen::MatrixXd tar = tget.transpose(); - - Geometry rotated = tar.transpose() * R; - - Molecule ref_mol = m_reference; - ref_mol.setGeometry(ref); - - Molecule tar_mol = m_target; - tar_mol.setGeometry(rotated); - - auto blob = MakeCostMatrix(cached_reference, rotated); - m_cost_limit = blob.first; - if (!m_silent) - - if (m_method == 4) - blob.first = 0; - m_prepared_cost_matrices.insert(blob); - return true; -} - bool RMSDDriver::MolAlignLib() { m_reference.writeXYZFile("molaign_ref.xyz"); diff --git a/src/capabilities/rmsd.h b/src/capabilities/rmsd.h index 1890c20..54d5e6b 100644 --- a/src/capabilities/rmsd.h +++ b/src/capabilities/rmsd.h @@ -31,6 +31,8 @@ #include #include +#include + #include "json.hpp" using json = nlohmann::json; @@ -67,28 +69,7 @@ class RMSDThread : public CxxThread { int m_calculations = 0; std::function m_evaluator; }; -/* -class IntermediateStorage { -public: - inline IntermediateStorage(unsigned int size) - : m_size(size) - { - } - - inline void addItem(const std::vector& vector, double rmsd) - { - m_shelf.insert(std::pair>(rmsd, vector)); - if (m_shelf.size() >= m_size) - m_shelf.erase(--m_shelf.end()); - } - - const std::map>* data() const { return &m_shelf; } -private: - unsigned int m_size; - std::map> m_shelf; -}; -*/ static const json RMSDJson = { { "reorder", false }, { "check", false }, @@ -100,7 +81,7 @@ static const json RMSDJson = { { "pt", 0 }, { "silent", false }, { "storage", 1.0 }, - { "method", "incr" }, + { "method", "free" }, { "noreorder", false }, { "threads", 1 }, { "Element", 7 }, @@ -119,7 +100,9 @@ static const json RMSDJson = { { "molaligntol", 10 }, { "cycles", -1 }, { "nofree", false }, - { "limit", 10 } + { "limit", 10 }, + { "costmatrix", 1 }, + { "maxtrial", 3 } }; class RMSDDriver : public CurcumaMethod { @@ -242,6 +225,7 @@ class RMSDDriver : public CurcumaMethod { void setThreads(int threads) { m_threads = threads; } bool MolAlignLib(); + static std::pair MakeCostMatrix(const Geometry& reference, const Geometry& target, const std::vector& reference_atoms, const std::vector& target_atoms, int costmatrix); private: /* Read Controller has to be implemented for all */ @@ -270,6 +254,8 @@ class RMSDDriver : public CurcumaMethod { void CheckTopology(); + Matrix OptimiseRotation(const Eigen::Matrix3d& rotation); + std::pair, std::vector> PrepareHeavyTemplate(); std::pair, std::vector> PrepareDistanceTemplate(); @@ -290,6 +276,7 @@ class RMSDDriver : public CurcumaMethod { } std::vector FillMissing(const Molecule& molecule, const std::vector& order); + void InsertRotation(std::pair& rotation); void InitialiseOrder(); std::pair InitialisePair(); @@ -324,7 +311,6 @@ class RMSDDriver : public CurcumaMethod { std::vector m_reorder_rules; std::vector> m_stored_rules, m_intermedia_rules; std::vector m_tmp_rmsd; - std::map> m_connectivity; double m_rmsd = 0, m_rmsd_raw = 0, m_scaling = 1.5, m_intermedia_storage = 1, m_threshold = 99, m_damping = 0.8, m_dmix = -1; bool m_check = false; bool m_check_connections = false, m_postprocess = true, m_noreorder = false, m_swap = false, m_dynamic_center = false; @@ -334,10 +320,73 @@ class RMSDDriver : public CurcumaMethod { int m_munkress_cycle = 1; int m_molaligntol = 10; int m_limit = 10; + int m_costmatrix = 1; + int m_maxtrial = 2; double m_cost_limit = 0; mutable int m_fragment = -1, m_fragment_reference = -1, m_fragment_target = -1; std::vector m_initial, m_element_templates; std::string m_molalign = "molalign"; - std::map> m_distance_reference, m_distance_target; std::map m_prepared_cost_matrices; }; + +using namespace LBFGSpp; + +class LBFGSRotation { +public: + LBFGSRotation(int n_) + { + } + double operator()(const Eigen::VectorXd& x, Eigen::VectorXd& grad) + { + Eigen::Matrix3d n; + n = Eigen::AngleAxisd(x[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxisd(x[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxisd(x[2], Eigen::Vector3d::UnitZ()); + + // Eigen::MatrixXd tar = m_target.transpose(); + + // Geometry rotated = tar.transpose() * n; + double fx = 0.0; + double dx = 1e-5; + + auto result = RMSDDriver::MakeCostMatrix(m_reference, m_target * n, m_reference_atoms, m_target_atoms, m_costmatrix); + // std::cout << result.first << std::endl; + Eigen::VectorXd tmp = x; + for (int i = 0; i < 3; ++i) { + tmp[i] += dx; + // std::cout << tmp.transpose() << std::endl; + Eigen::Matrix3d n; + n = Eigen::AngleAxisd(tmp[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxisd(tmp[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxisd(tmp[2], Eigen::Vector3d::UnitZ()); + + // Geometry rotated = tar.transpose() * n; + + auto p = RMSDDriver::MakeCostMatrix(m_reference, m_target * n, m_reference_atoms, m_target_atoms, m_costmatrix).first; + + tmp[i] -= 2 * dx; + // std::cout << tmp.transpose() << std::endl; + + n = Eigen::AngleAxisd(tmp[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxisd(tmp[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxisd(tmp[2], Eigen::Vector3d::UnitZ()); + + // rotated = tar.transpose() * n; + + auto m = RMSDDriver::MakeCostMatrix(m_reference, m_target * n, m_reference_atoms, m_target_atoms, m_costmatrix).first; + // std::cout << p << " " << m << " " << (p-m)/(2*dx) << std::endl << std::endl; + grad[i] = (p - m) / (2 * dx); + } + + return result.first; + } + + Vector Parameter() const { return m_parameter; } + + Geometry m_reference, m_target; + std::vector m_reference_atoms, m_target_atoms; + int m_costmatrix; + +private: + Vector m_parameter; +}; diff --git a/src/core/molecule.cpp b/src/core/molecule.cpp index f1704b3..2635e44 100644 --- a/src/core/molecule.cpp +++ b/src/core/molecule.cpp @@ -957,6 +957,19 @@ Position Molecule::CalculateDipoleMoment(const std::vector& scaling) con dipole(0) += m_charges[i] * (m_geometry(i, 0) - pos(0)) * scale; dipole(1) += m_charges[i] * (m_geometry(i, 1) - pos(1)) * scale; dipole(2) += m_charges[i] * (m_geometry(i, 2) - pos(2)) * scale; + for (int j = 0; j < m_geometry.rows(); ++j) { + if (i == j) + continue; + double scale = 3; + if (scaling.size() > i) + scale = scaling[i]; + double revr = 1 / (m_geometry.row(i) - m_geometry.row(j)).norm(); + // std::cout << m_charges[i] *revr * scale << " "; + + dipole(0) -= m_charges[i] * revr * scale; + dipole(1) -= m_charges[i] * revr * scale; + dipole(2) -= m_charges[i] * revr * scale; + } // std::cout << scale << " "; } // std::cout << std::endl;