diff --git a/CMakeLists.txt b/CMakeLists.txt index 453bbc0..f975206 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -164,9 +164,9 @@ set(curcuma_core_SRC src/core/eigen_uff.cpp src/core/qmdff.cpp src/core/eht.cpp - #src/core/forcefield.h - #src/core/forcefield.cpp - #src/core/forcefieldfunctions.h + src/core/forcefieldthread.cpp + src/core/forcefield.cpp + src/core/forcefieldfunctions.h #src/core/forcefield_terms/qmdff_terms.h src/tools/formats.h src/tools/geometry.h @@ -277,12 +277,12 @@ if(USE_D3) if(NOT TARGET s-dftd3) add_subdirectory(external/simple-dftd3) endif() - #if(HELPERS) + if(HELPERS) add_executable(dftd3_helper src/helpers/dftd3_helper.cpp ) target_link_libraries(dftd3_helper PUBLIC s-dftd3 gfortran pthread) - #endif() + endif() set(curcuma_D3_SRC src/core/dftd3interface.cpp ) diff --git a/src/capabilities/rmsd.cpp b/src/capabilities/rmsd.cpp index f78fd66..2234bdd 100644 --- a/src/capabilities/rmsd.cpp +++ b/src/capabilities/rmsd.cpp @@ -1,6 +1,6 @@ /* * - * Copyright (C) 2019 - 2023 Conrad Hübler + * Copyright (C) 2019 - 2024 Conrad Hübler * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -803,32 +804,35 @@ void RMSDDriver::HeavyTemplate() void RMSDDriver::TemplateFree() { + Molecule ref_mol = m_reference; Molecule tar_mol = m_target; + for (int i = 0; i < 2; ++i) { - Geometry cached_reference = m_reference.getGeometry(); - Geometry cached_target = m_target.getGeometry(); - - Geometry tref = GeometryTools::TranslateMolecule(m_reference, m_reference.Centroid(true), Position{ 0, 0, 0 }); - Geometry tget = GeometryTools::TranslateMolecule(m_target, m_target.Centroid(true), Position{ 0, 0, 0 }); - ref_mol.setGeometry(tref); - tar_mol.setGeometry(tget); - - if (m_moi) { - ref_mol.CalculateRotationalConstants(); - tar_mol.CalculateRotationalConstants(); - - Eigen::MatrixXd tar = tget.transpose(); - Eigen::MatrixXd ref = tref.transpose(); - - Geometry rotated_reference = ref.transpose() * ref_mol.AlignmentAxes(); - Geometry rotated_target = tar.transpose() * tar_mol.AlignmentAxes(); + Geometry cached_reference = m_reference.getGeometry(); + Geometry cached_target = m_target.getGeometry(); - ref_mol.setGeometry(rotated_reference); - tar_mol.setGeometry(rotated_target); - ref_mol.writeXYZFile("reference.moi.xyz"); - tar_mol.writeXYZFile("target.moi.xyz"); - } else { + Geometry tref = GeometryTools::TranslateMolecule(m_reference, m_reference.Centroid(true), Position{ 0, 0, 0 }); + Geometry tget = GeometryTools::TranslateMolecule(m_target, m_target.Centroid(true), Position{ 0, 0, 0 }); + ref_mol.setGeometry(tref); + tar_mol.setGeometry(tget); + /* + if (m_moi) { + ref_mol.CalculateRotationalConstants(); + tar_mol.CalculateRotationalConstants(); + + Eigen::MatrixXd tar = tget.transpose(); + Eigen::MatrixXd ref = tref.transpose(); + + Geometry rotated_reference = ref.transpose() * ref_mol.AlignmentAxes(); + Geometry rotated_target = tar.transpose() * tar_mol.AlignmentAxes(); + + ref_mol.setGeometry(rotated_reference); + tar_mol.setGeometry(rotated_target); + ref_mol.writeXYZFile("reference.moi.xyz"); + tar_mol.writeXYZFile("target.moi.xyz"); + } else { + */ auto operators = GetOperateVectors(ref_mol, tar_mol); Eigen::Matrix3d R = operators.first; Eigen::MatrixXd tar = tget.transpose(); @@ -839,14 +843,17 @@ void RMSDDriver::TemplateFree() ref_mol.setGeometry(tref); Molecule tar_mol = m_target; tar_mol.setGeometry(rotated); - ref_mol.writeXYZFile("reference.nomoi.xyz"); - tar_mol.writeXYZFile("target.nomoi.xyz"); + // ref_mol.writeXYZFile("reference.nomoi.xyz"); + // tar_mol.writeXYZFile("target.nomoi.xyz"); + // } + + std::vector new_order = DistanceReorder(ref_mol, tar_mol); + auto rules = ApplyOrder(new_order, m_target); + double rmsd = Rules2RMSD(new_order); + m_target = rules; + m_target_aligned = m_target; + m_reorder_rules = new_order; } - std::vector new_order = DistanceReorder(ref_mol, tar_mol); - m_target_reordered = ApplyOrder(new_order, m_target); - m_target = m_target_reordered; - m_target_aligned = m_target; - m_reorder_rules = new_order; } void RMSDDriver::FinaliseTemplate(std::pair, std::vector> pairs) @@ -863,13 +870,15 @@ void RMSDDriver::FinaliseTemplate(std::pair, std::vector> auto result = AlignByVectorPair(pairs); m_reorder_rules = result; m_target_reordered = ApplyOrder(m_reorder_rules, target); - local_results.insert(std::pair>(Rules2RMSD(m_reorder_rules), m_reorder_rules)); + double rmsd = Rules2RMSD(m_reorder_rules); + // std::cout << rmsd << " "; + local_results.insert(std::pair>(rmsd, m_reorder_rules)); /*std::set s(result.second.begin(), result.second.end()); if (s.size() != result.second.size()) // make sure, that only results with non-duplicate vectors are accepted continue;*/ result = AlignByVectorPair(tmp, m_reorder_rules); - if (m_reorder_rules == result) - break; + // if (m_reorder_rules == result) + // break; m_reorder_rules = result; m_target_reordered = ApplyOrder(m_reorder_rules, target); StructComp structcomp = Rule2RMSD(m_reorder_rules); @@ -1158,7 +1167,7 @@ bool RMSDDriver::TemplateReorder() std::vector RMSDDriver::DistanceReorder(const Molecule& reference, const Molecule& target) { std::map> rules; - + /* std::vector orderV1 = DistanceReorderV1(reference, target); std::vector orderV2 = DistanceReorderV2(reference, target); @@ -1167,24 +1176,25 @@ std::vector RMSDDriver::DistanceReorder(const Molecule& reference, const Mo rules.insert(std::pair>(rmsdV1, orderV1)); rules.insert(std::pair>(rmsdV2, orderV2)); - +*/ if (m_nomunkres == false) { std::vector munkress = Munkress(reference, target); double rmsdM = Rules2RMSD(munkress); + // std::cout << rmsdM << ": " << std::endl; rules.insert(std::pair>(rmsdM, munkress)); } - - if (m_update_rotation) { - auto orderV3 = DistanceReorderV3(reference, target); - double rmsdV3 = Rules2RMSD(orderV3.first); - - rules.insert(std::pair>(rmsdV3, orderV3.first)); - if (m_split) { - double rmsdV4 = Rules2RMSD(orderV3.second); - rules.insert(std::pair>(rmsdV4, orderV3.first)); + /* + if (m_update_rotation) { + auto orderV3 = DistanceReorderV3(reference, target); + double rmsdV3 = Rules2RMSD(orderV3.first); + + rules.insert(std::pair>(rmsdV3, orderV3.first)); + if (m_split) { + double rmsdV4 = Rules2RMSD(orderV3.second); + rules.insert(std::pair>(rmsdV4, orderV3.first)); + } } - } - + */ for (auto iter : rules) m_stored_rules.push_back(iter.second); @@ -1391,19 +1401,25 @@ std::vector RMSDDriver::Munkress(const Molecule& reference, const Molecule& Eigen::MatrixXd distance = Eigen::MatrixXd::Zero(reference.AtomCount(), reference.AtomCount()); std::vector element_reference = reference.Atoms(); std::vector element_target = target.Atoms(); - + double min = penalty; + static std::random_device rd{}; + static std::mt19937 gen{ rd() }; + static std::normal_distribution<> d{ 0, 3 }; for (int i = 0; i < reference.AtomCount(); ++i) { for (int j = 0; j < target.AtomCount(); ++j) { - distance(i, j) = GeometryTools::Distance(target.Atom(j).second, reference.Atom(i).second) + penalty * (target.Atom(j).first != reference.Atom(i).first); + + distance(i, j) = GeometryTools::Distance(target.Atom(j).second, reference.Atom(i).second) * GeometryTools::Distance(target.Atom(j).second, reference.Atom(i).second) + // + (target.Atom(j).second.norm() - reference.Atom(i).second.norm())*(target.Atom(j).second.norm() - reference.Atom(i).second.norm()) + // + d(gen) + + penalty * (target.Atom(j).first != reference.Atom(i).first); } } if (m_dmix <= 1 && 0 < m_dmix) { Matrix d = target.DistanceMatrix().first; distance = (1 - m_dmix) * distance + m_dmix * d; } - // std::cout << distance << std::endl; + auto result = MunkressAssign(distance); - // std::cout << result << std::endl; for (int i = 0; i < result.cols(); ++i) { for (int j = 0; j < result.rows(); ++j) { @@ -1413,9 +1429,6 @@ std::vector RMSDDriver::Munkress(const Molecule& reference, const Molecule& } } } - // for(auto i : new_order) - // std::cout << i << " "; - // std::cout << std::endl; return new_order; } @@ -1425,7 +1438,7 @@ bool RMSDDriver::MolAlignLib() m_target.writeXYZFile("molalign_tar.xyz"); FILE* FileOpen; - std::string command = m_molalign + " molaign_ref.xyz " + " molalign_tar.xyz -reorder -fast -tol " + std::to_string(m_molaligntol) + " 2>&1"; + std::string command = m_molalign + " molaign_ref.xyz " + " molalign_tar.xyz -sort -fast -tol " + std::to_string(m_molaligntol) + " 2>&1"; if (!m_silent) std::cout << command << std::endl; FileOpen = popen(command.c_str(), "r"); @@ -1447,10 +1460,11 @@ bool RMSDDriver::MolAlignLib() fmt::print(fg(fmt::color::green) | fmt::emphasis::bold, "\nPlease cite the follow research report!\nJ. Chem. Inf. Model. 2023, 63, 4, 1157–1165 - DOI: 10.1021/acs.jcim.2c01187\n\n"); } FileIterator file("aligned.xyz", true); + file.Next(); m_target_reordered = file.Next(); m_target_aligned = m_target_reordered; m_target = m_target_reordered; - std::filesystem::remove("aligned.xyz"); + // std::filesystem::remove("aligned.xyz"); } else { if (!rndm && !error) { fmt::print(fg(fmt::color::salmon) | fmt::emphasis::bold, "Molalign was not found. Consider getting it from\nhttps://github.com/qcuaeh/molalignlib\nEither adding the location of the binary to the path for executables or append\n-molalignbin /yourpath/molalign\nto your argument list!\n"); diff --git a/src/core/forcefield.cpp b/src/core/forcefield.cpp index 4b55e14..7d37758 100644 --- a/src/core/forcefield.cpp +++ b/src/core/forcefield.cpp @@ -1,6 +1,6 @@ /* * < Generic force field class for curcuma . > - * Copyright (C) 2022 - 2023 Conrad Hübler + * Copyright (C) 2024 Conrad Hübler * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -21,173 +21,167 @@ #include "src/core/qmdff_par.h" #include "src/core/uff_par.h" +#include "forcefieldfunctions.h" + #include "forcefield.h" +#include "forcefieldthread.h" -ForceFieldThread::ForceFieldThread() +ForceField::ForceField(const json& controller) { } -double ForceFieldThread::LJBondStretching() +void ForceField::UpdateGeometry(const Matrix& geometry) { - double factor = m_final_factor * m_bond_scaling; - double energy = 0.0; - Eigen::Vector3d dx = { m_d, 0, 0 }; - Eigen::Vector3d dy = { 0, m_d, 0 }; - Eigen::Vector3d dz = { 0, 0, m_d }; - - for (int index = 0; index < m_bonds.size(); ++index) { - const auto& bond = m_bonds[index]; - - Vector i = Position(bond.i); - Vector j = Position(bond.j); - - // Matrix derivate; - energy += StretchEnergy(i - j, bond.r0_ij, bond.fc, bond.exponent) * factor; - if (m_calculate_gradient) { - /*if (m_calc_gradient == 0) { - double diff = 0; - m_gradient.row(a) += diff * derivate.row(0); - m_gradient.row(b) += diff * derivate.row(1); - - } else if (m_calc_gradient == 1) {*/ - double d_x = (StretchEnergy(i + dx - j, bond.r0_ij, bond.fc, bond.exponent) - StretchEnergy(i - dx - j, bond.r0_ij, bond.fc, bond.exponent)) / (2 * m_d) * factor; - double d_y = (StretchEnergy(i + dy - j, bond.r0_ij, bond.fc, bond.exponent) - StretchEnergy(i - dy - j, bond.r0_ij, bond.fc, bond.exponent)) / (2 * m_d) * factor; - double d_z = (StretchEnergy(i + dz - j, bond.r0_ij, bond.fc, bond.exponent) - StretchEnergy(i - dz - j, bond.r0_ij, bond.fc, bond.exponent)) / (2 * m_d) * factor; - m_gradient(bond.a, 0) += d_x; - m_gradient(bond.a, 1) += d_y; - m_gradient(bond.a, 2) += d_z; - - m_gradient(bond.b, 0) -= d_x; - m_gradient(bond.b, 1) -= d_y; - m_gradient(bond.b, 2) -= d_x; - //} - } - } + m_geometry = geometry; +} - return energy; +void ForceField::UpdateGeometry(const double* coord) +{ +#pragma message("replace with raw data") + for (int i = 0; i < m_natoms; ++i) { + m_geometry(i, 0) = coord[3 * i + 0]; + m_geometry(i, 1) = coord[3 * i + 1]; + m_geometry(i, 2) = coord[3 * i + 2]; + } } -double ForceFieldThread::HarmonicBondStretching() +void ForceField::UpdateGeometry(const std::vector>& geometry) { - double factor = m_final_factor * m_bond_scaling; - double energy = 0.0; - - for (int index = 0; index < m_bonds.size(); ++index) { - const auto& bond = m_bonds[index]; - const int i = bond.i; - const int j = bond.j; - - Vector x = Position(i); - Vector y = Position(j); - Matrix derivate; - double r0 = BondStretching(x, y, derivate, m_calculate_gradient); - - energy += (0.5 * bond.fc * (r0 - bond.r0_ij) * (r0 - bond.r0_ij)) * factor; - if (m_calculate_gradient) { - double diff = (bond.fc) * (r0 - bond.r0_ij) * factor; - m_gradient.row(i) += diff * derivate.row(0); - m_gradient.row(j) += diff * derivate.row(1); - } +#pragma message("replace with raw data") + for (int i = 0; i < m_natoms; ++i) { + m_geometry(i, 0) = geometry[i][0]; + m_geometry(i, 1) = geometry[i][1]; + m_geometry(i, 2) = geometry[i][2]; } +} - return energy; +void ForceField::setParameter(const json& parameter) +{ } -double QMDFFThread::CalculateAngleBending() +void ForceField::setBonds(const json& bonds) { - double threshold = 1e-2; - double energy = 0.0; - Eigen::Vector3d dx = { m_d, 0, 0 }; - Eigen::Vector3d dy = { 0, m_d, 0 }; - Eigen::Vector3d dz = { 0, 0, m_d }; - for (int index = 0; index < m_qmdffangle.size(); ++index) { - const auto& angle = m_qmdffangle[index]; - const int a = angle.a; - const int b = angle.b; - const int c = angle.c; - - auto atom_a = Position(a); - auto atom_b = Position(b); - auto atom_c = Position(c); - Matrix derivate; - double costheta = AngleBending(atom_a, atom_b, atom_c, derivate, m_CalculateGradient); - std::function angle_function; - if (std::abs(costheta - pi) < threshold) - - angle_function = [this](const Eigen::Vector3d& a, const Eigen::Vector3d& b, const Eigen::Vector3d& c, double thetae, double kabc, double reAB, double reAC) -> double { - double val = this->LinearAngleBend(a, b, c, thetae, kabc, reAB, reAC); - if (std::isnan(val)) - return 0; - else - return val; - }; - else - angle_function = [this](const Eigen::Vector3d& a, const Eigen::Vector3d& b, const Eigen::Vector3d& c, double thetae, double kabc, double reAB, double reAC) -> double { - double val = this->AngleBend(a, b, c, thetae, kabc, reAB, reAC); - if (std::isnan(val)) - return 0; - else - return val; - }; - - double e = angle_function(atom_a, atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC); //(angle.kijk * (angle.C0 + angle.C1 * costheta + angle.C2 * (2 * costheta * costheta - 1))) * m_final_factor * m_angle_scaling; - - if (m_CalculateGradient) { - if (m_calc_gradient == 0) { - /* - double sintheta = sin(acos(costheta)); - double dEdtheta = -angle.kijk * sintheta * (angle.C1 + 4 * angle.C2 * costheta) * m_final_factor * m_angle_scaling; - m_gradient.row(i) += dEdtheta * derivate.row(0); - m_gradient.row(j) += dEdtheta * derivate.row(1); - m_gradient.row(k) += dEdtheta * derivate.row(2); - */ - } else if (m_calc_gradient == 1) { - m_gradient(a, 0) += (angle_function(AddVector(atom_a, dx), atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(SubVector(atom_a, dx), atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - m_gradient(a, 1) += (angle_function(AddVector(atom_a, dy), atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(SubVector(atom_a, dy), atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - m_gradient(a, 2) += (angle_function(AddVector(atom_a, dz), atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(SubVector(atom_a, dz), atom_b, atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - - m_gradient(b, 0) += (angle_function(atom_a, AddVector(atom_b, dx), atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(atom_a, SubVector(atom_b, dx), atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - m_gradient(b, 1) += (angle_function(atom_a, AddVector(atom_b, dy), atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(atom_a, SubVector(atom_b, dy), atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - m_gradient(b, 2) += (angle_function(atom_a, AddVector(atom_b, dz), atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(atom_a, SubVector(atom_b, dz), atom_c, angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - - m_gradient(c, 0) += (angle_function(atom_a, atom_b, AddVector(atom_c, dx), angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(atom_a, atom_b, SubVector(atom_c, dx), angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - m_gradient(c, 1) += (angle_function(atom_a, atom_b, AddVector(atom_c, dy), angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(atom_a, atom_b, SubVector(atom_c, dy), angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - m_gradient(c, 2) += (angle_function(atom_a, atom_b, AddVector(atom_c, dz), angle.thetae, angle.kabc, angle.reAB, angle.reAC) - angle_function(atom_a, atom_b, SubVector(atom_c, dz), angle.thetae, angle.kabc, angle.reAB, angle.reAC)) / (2 * m_d); - } - } + m_bonds.clear(); + for (int i = 0; i < bonds.size(); ++i) { + json bond = bonds[i].get(); + Bond b; + b.type = bond["type"]; + b.i = bond["i"]; + b.j = bond["j"]; + b.k = bond["k"]; + b.distance = bond["distance"]; + + b.r0_ij = bond["r0_ij"]; + b.r0_ik = bond["r0_ik"]; + + b.fc = bond["fc"]; + m_bonds.push_back(b); } - return energy; } -ForceField::ForceField(const json& controller) +void ForceField::setAngles(const json& angles) { + m_angles.clear(); + for (int i = 0; i < angles.size(); ++i) { + json angle = angles[i].get(); + Angle a; + + a.type = angle["type"]; + + a.i = angle["i"]; + a.j = angle["j"]; + a.k = angle["k"]; + a.C0 = angle["C0"]; + a.C1 = angle["C1"]; + a.C2 = angle["C2"]; + a.fc = angle["fc"]; + a.r0_ij = angle["r0_ij"]; + a.r0_ik = angle["r0_ik"]; + a.theta0_ijk = angle["theta0_ijk"]; + m_angles.push_back(a); + } } -void ForceField::UpdateGeometry(const Matrix& geometry) +void ForceField::setDihedrals(const json& dihedrals) { - m_geometry = geometry; + m_dihedrals.clear(); + for (int i = 0; i < dihedrals.size(); ++i) { + json dihedral = dihedrals[i].get(); + Dihedral d; + d.type = dihedral["type"]; + + d.i = dihedral["i"]; + d.j = dihedral["j"]; + d.k = dihedral["k"]; + d.l = dihedral["l"]; + d.V = dihedral["V"]; + d.n = dihedral["n"]; + d.phi0 = dihedral["phi0"]; + m_dihedrals.push_back(d); + } } -void ForceField::UpdateGeometry(const double* coord) +void ForceField::setInversions(const json& inversions) { -#pragma message("replace with raw data") - for (int i = 0; i < m_natoms; ++i) { - m_geometry(i, 0) = coord[3 * i + 0]; - m_geometry(i, 1) = coord[3 * i + 1]; - m_geometry(i, 2) = coord[3 * i + 2]; + m_inversions.clear(); + for (int i = 0; i < inversions.size(); ++i) { + json inversion = inversions[i].get(); + Inversion inv; + inv.type = inversion["type"]; + + inv.i = inversion["i"]; + inv.j = inversion["j"]; + inv.k = inversion["k"]; + inv.l = inversion["l"]; + inv.fc = inversion["fc"]; + inv.C0 = inversion["C0"]; + inv.C1 = inversion["C1"]; + inv.C2 = inversion["C2"]; + + m_inversions.push_back(inv); } } -void ForceField::UpdateGeometry(const std::vector>& geometry); +void ForceField::setvdWs(const json& vdws) { -#pragma message("replace with raw data") - for (int i = 0; i < m_natoms; ++i) { - m_geometry(i, 0) = geometry[i][0]; - m_geometry(i, 1) = geometry[i][1]; - m_geometry(i, 2) = geometry[i][2]; + m_vdWs.clear(); + for (int i = 0; i < vdws.size(); ++i) { + json vdw = vdws[i].get(); + vdW v; + v.type = vdw["type"]; + + v.i = vdw["i"]; + v.j = vdw["j"]; + v.C_ij = vdw["C_ij"]; + v.r0_ij = vdw["r0_ij"]; + + m_vdWs.push_back(v); } } -void ForceField::setParameter(const json& parameter) +void ForceField::AutoRanges() { + for (int i = 0; i < m_threads; ++i) { + ForceFieldThread* thread = new ForceFieldThread(i, m_threads); + thread->setGeometry(m_geometry); + m_threadpool->addThread(thread); + m_stored_threads.push_back(thread); + for (int j = int(i * m_bonds.size() / double(m_threads)); j < int((i + 1) * m_bonds.size() / double(m_threads)); ++j) + thread->addBond(m_bonds[j]); + + for (int j = int(i * m_angles.size() / double(m_threads)); j < int((i + 1) * m_angles.size() / double(m_threads)); ++j) + thread->addAngle(m_angles[j]); + + for (int j = int(i * m_dihedrals.size() / double(m_threads)); j < int((i + 1) * m_dihedrals.size() / double(m_threads)); ++j) + thread->addDihedral(m_dihedrals[j]); + + for (int j = int(i * m_inversions.size() / double(m_threads)); j < int((i + 1) * m_inversions.size() / double(m_threads)); ++j) + thread->addInversion(m_inversions[j]); + + for (int j = int(i * m_vdWs.size() / double(m_threads)); j < int((i + 1) * m_vdWs.size() / double(m_threads)); ++j) + thread->addvdW(m_vdWs[j]); + + for (int j = int(i * m_EQs.size() / double(m_threads)); j < int((i + 1) * m_EQs.size() / double(m_threads)); ++j) + thread->addEQ(m_EQs[j]); + } } diff --git a/src/core/forcefield.h b/src/core/forcefield.h index 42d221a..9059643 100644 --- a/src/core/forcefield.h +++ b/src/core/forcefield.h @@ -1,6 +1,6 @@ /* * < Generic force field class for curcuma . > - * Copyright (C) 2022 - 2023 Conrad Hübler + * Copyright (C) 2024 Conrad Hübler * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -21,6 +21,8 @@ #include "src/core/global.h" +#include "forcefieldthread.h" + #include "hbonds.h" #include "external/CxxThreadPool/include/CxxThreadPool.h" @@ -36,6 +38,7 @@ #include "src/core/qmdff_par.h" #include "src/core/uff_par.h" +#include #include #include @@ -44,40 +47,6 @@ #include "json.hpp" using json = nlohmann::json; -struct Bond { - int i = 0, j = 0, k = 0, distance = 0; - double fc = 0, exponent = 0, r0_ij = 0, r0_ik = 0; -}; - -struct Angle { - int i = 0, j = 0, k = 0; - double fc = 0, r0_ij = 0, r0_ik = 0, theta0_ijk = 0; -}; - -class ForceFieldThread : public CxxThread { - -public: - ForceFieldThread(); - -private: - double CalculateBondContribution(); - - double HarmonicBondStretching(); - - double LJBondStretching(); - - double m_final_factor = 1; - double m_bond_scaling = 1; - double m_au = 1; - double m_d = 1e-3; - int m_calc_gradient = 1; - bool m_calculate_gradient = true; - - Matrix m_geometry, m_gradient; - - std::vector m_bonds; -}; - class ForceField { public: @@ -93,8 +62,26 @@ class ForceField { Matrix Gradient() const { return m_gradient; } void setParameter(const json& parameter); + void setParameterFile(const std::string& file); private: + void AutoRanges(); + void setBonds(const json& bonds); + void setAngles(const json& angles); + void setDihedrals(const json& dihedrals); + void setInversions(const json& inversions); + + std::vector m_stored_threads; + CxxThreadPool* m_threadpool; + void setvdWs(const json& vdws); + Matrix m_geometry, m_gradient; int m_natoms = 0; + int m_threads = 1; + std::vector m_bonds; + std::vector m_angles; + std::vector m_dihedrals; + std::vector m_inversions; + std::vector m_vdWs; + std::vector m_EQs; }; diff --git a/src/core/forcefieldfunctions.h b/src/core/forcefieldfunctions.h index 23f1795..aff1146 100644 --- a/src/core/forcefieldfunctions.h +++ b/src/core/forcefieldfunctions.h @@ -1,4 +1,146 @@ -#ifndef FORCEFIELDFUNCTIONS_H -#define FORCEFIELDFUNCTIONS_H +/* + * < Generic force field class for curcuma . > + * Copyright (C) 2022 - 2023 Conrad Hübler + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + */ -#endif // FORCEFIELDFUNCTIONS_H +#pragma once + +#include "src/core/qmdff_par.h" +#include "src/core/uff_par.h" + +#include +#include +#include + +#include + +#include "json.hpp" + +namespace UFF { +inline double BondStretching(const Vector& i, const Vector& j, Matrix& derivate, bool gradient) +{ + Vector ij = i - j; + double distance = (ij).norm(); + if (!gradient) + return distance; + derivate = Matrix::Zero(2, 3); + derivate.row(0) = ij / distance; + derivate.row(1) = -ij / distance; + return distance; +} + +inline double AngleBending(const Vector& i, const Vector& j, const Vector& k, Matrix& derivate, bool gradient) +{ + Eigen::Vector3d rij = i - j; + auto nij = rij / rij.norm(); + Eigen::Vector3d rkj = k - j; + auto nkj = rkj / rkj.norm(); + double costheta = (rij.dot(rkj) / (sqrt(rij.dot(rij) * rkj.dot(rkj)))); + + if (!gradient) + return costheta; + derivate = Matrix::Zero(3, 3); + + double sintheta = sin(acos(costheta)); + double dThetadCosTheta = 1 / sintheta; + derivate = Matrix::Zero(3, 3); + derivate.row(0) = -dThetadCosTheta * (nkj - nij * costheta) / (rij.norm()); + derivate.row(2) = -dThetadCosTheta * (nij - nkj * costheta) / (rkj.norm()); + derivate.row(1) = -derivate.row(0) - derivate.row(2); + + return costheta; +} + +inline Eigen::Vector3d NormalVector(const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k) +{ + return (j - i).cross(j - k); +} +/* +inline double Dihedral(const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, const Eigen::Vector3d& l, double V, double n, double phi0) +{ + Eigen::Vector3d nijk = NormalVector(i, j, k); + Eigen::Vector3d njkl = NormalVector(j, k, l); + double n_ijk = (nijk).norm(); + double n_jkl = (njkl).norm(); + double dotpr = nijk.dot(njkl); + Eigen::Vector3d ji = j - i; + double sign = (-1 * ji).dot(njkl) < 0 ? -1 : 1; + double phi = pi + sign * acos(dotpr / (n_ijk * n_jkl)); + double energy = (1 / 2.0 * V * (1 - cos(n * phi0) * cos(n * phi))) * m_final_factor * m_dihedral_scaling; + if (isnan(energy)) + return 0; + else + return energy; +} */ +} + +namespace QMDFF { + +double LJStretchEnergy(double r_ij, double r0_ij, double fc, double exponent) +{ + const double ratio = r0_ij / r_ij; + double energy = fc * (1 + pow(ratio, exponent) - 2 * pow(ratio, exponent * 0.5)); + if (std::isnan(energy)) + return 0; + else + return energy; +} + +inline double AngleDamping(double r_ij, double r_ik, double r0_ij, double r0_ik) +{ + double kdamp = 1; + double ratioij = r_ij / r0_ij; + double ratioik = r_ik / r0_ik; + + double fijinv = 1 + kdamp * pow(ratioij, 4); + double fikinv = 1 + kdamp * pow(ratioik, 4); + double finv = fijinv * fikinv; + return 1 / finv; +} + +inline double AngleBend(const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, double theta0_ijk, double fc, double r0_ij, double r0_ik) +{ + Eigen::Vector3d ij = i - j; + Eigen::Vector3d ik = i - k; + double damp = AngleDamping(ij.norm(), ik.norm(), r0_ij, r0_ik); + double costheta = (ij.dot(ik) / (sqrt(ij.dot(ij) * ik.dot(ik)))); + double costheta0_ijk = cos(theta0_ijk); + double energy = (fc * damp * (costheta0_ijk - costheta) * (costheta0_ijk - costheta)); + if (std::isnan(energy)) + return 0; + else + return energy; +} + +inline double LinearAngleBend(const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, double theta0_ijk, double fc, double r0_ij, double r0_ik) +{ +#pragma message("check") + /* + if (kabc < 0) + return 0; +*/ + Eigen::Vector3d ij = i - j; + Eigen::Vector3d ik = i - k; + double damp = AngleDamping(ij.norm(), ik.norm(), r0_ij, r0_ik); + double theta = acos((ij.dot(ik) / (sqrt(ij.dot(ij) * ik.dot(ik))))); + double energy = (fc * damp * (theta0_ijk - theta) * (theta0_ijk - theta)); + if (std::isnan(energy)) + return 0; + else + return energy; +} +} diff --git a/src/core/forcefieldthread.cpp b/src/core/forcefieldthread.cpp new file mode 100644 index 0000000..e7ba352 --- /dev/null +++ b/src/core/forcefieldthread.cpp @@ -0,0 +1,429 @@ +/* + * < Generic force field class for curcuma . > + * Copyright (C) 2024 Conrad Hübler + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + */ + +#include "src/core/forcefieldderivaties.h" +#include "src/core/qmdff_par.h" +#include "src/core/uff_par.h" + +#include "forcefieldfunctions.h" + +#include "forcefield.h" + +ForceFieldThread::ForceFieldThread(int thread, int threads) + : m_thread(thread) + , m_threads(threads) +{ + setAutoDelete(false); + m_final_factor = 1 / 2625.15 * 4.19; + // m_d = parameters["differential"].get(); + m_d = 1e-7; +} + +int ForceFieldThread::execute() +{ + // m_CalculateGradient = grd; + m_d4_energy = 0; + m_d3_energy = 0; + + m_bond_energy = CalculateBondContribution(); + /* + m_angle_energy = CalculateAngleBending(); + m_dihedral_energy = CalculateDihedral(); + m_inversion_energy = CalculateInversion(); + m_vdw_energy = CalculateNonBonds(); + */ + /* + CalculateElectrostatic(); */ + m_energy = m_bond_energy + m_angle_energy + m_dihedral_energy + m_inversion_energy + m_vdw_energy; + return 0; +} + +void ForceFieldThread::addBond(const Bond& bonds) +{ + if (bonds.type == 1) + m_uff_bonds.push_back(bonds); + else if (bonds.type == 2) + m_qmdff_bonds.push_back(bonds); +} + +void ForceFieldThread::addAngle(const Angle& angles) +{ + if (angles.type == 1) + m_uff_angles.push_back(angles); + else if (angles.type == 2) + m_qmdff_angles.push_back(angles); +} + +void ForceFieldThread::addDihedral(const Dihedral& dihedrals) +{ + if (dihedrals.type == 1) + m_uff_dihedrals.push_back(dihedrals); + else if (dihedrals.type == 2) + m_qmdff_dihedrals.push_back(dihedrals); +} + +void ForceFieldThread::addInversion(const Inversion& inversions) +{ + if (inversions.type == 1) + m_uff_inversions.push_back(inversions); + else if (inversions.type == 2) + m_qmdff_inversions.push_back(inversions); +} + +void ForceFieldThread::addvdW(const vdW& vdWs) +{ + if (vdWs.type == 1) + m_uff_vdWs.push_back(vdWs); +} + +void ForceFieldThread::addEQ(const EQ& EQs) +{ + if (EQs.type == 2) + m_qmdff_EQs.push_back(EQs); +} + +double ForceFieldThread::CalculateUFFBondContribution() +{ + double factor = m_final_factor * m_bond_scaling; + double energy = 0.0; + + for (int index = 0; index < m_uff_bonds.size(); ++index) { + const auto& bond = m_uff_bonds[index]; + + Vector i = Position(bond.i); + Vector j = Position(bond.j); + Matrix derivate; + double rij = UFF::BondStretching(i, j, derivate, m_calculate_gradient); + + energy += (0.5 * bond.fc * (rij - bond.r0_ij) * (rij - bond.r0_ij)) * factor; + if (m_calculate_gradient) { + double diff = (bond.fc) * (rij - bond.r0_ij) * factor; + m_gradient.row(bond.i) += diff * derivate.row(0); + m_gradient.row(bond.j) += diff * derivate.row(1); + } + } + return energy; +} + +double ForceFieldThread::CalculateUFFAngleContribution() +{ + double energy = 0; + for (int index = 0; index < m_uff_angles.size(); ++index) { + const auto& angle = m_uff_angles[index]; + auto i = Position(angle.i); + auto j = Position(angle.j); + auto k = Position(angle.k); + Matrix derivate; + double costheta = UFF::AngleBending(i, j, k, derivate, m_calculate_gradient); + double e = (angle.fc * (angle.C0 + angle.C1 * costheta + angle.C2 * (2 * costheta * costheta - 1))) * m_final_factor * m_angle_scaling; + + if (m_calculate_gradient) { + if (m_calc_gradient == 0) { + double sintheta = sin(acos(costheta)); + double dEdtheta = -angle.fc * sintheta * (angle.C1 + 4 * angle.C2 * costheta) * m_final_factor * m_angle_scaling; + m_gradient.row(angle.i) += dEdtheta * derivate.row(0); + m_gradient.row(angle.j) += dEdtheta * derivate.row(1); + m_gradient.row(angle.k) += dEdtheta * derivate.row(2); + } + } + } + return energy; +} + +double ForceFieldThread::CalculateUFFDihedralContribution() +{ + double energy = 0; + for (int index = 0; index < m_uff_dihedrals.size(); ++index) { + const auto& dihedral = m_uff_dihedrals[index]; + Eigen::Vector3d i = Position(dihedral.i); + Eigen::Vector3d j = Position(dihedral.j); + Eigen::Vector3d k = Position(dihedral.k); + Eigen::Vector3d l = Position(dihedral.l); + Eigen::Vector3d eins = { 1.0, 1.0, 1.0 }; + Eigen::Vector3d nijk = UFF::NormalVector(i, j, k); + Eigen::Vector3d njkl = UFF::NormalVector(j, k, l); + + double n_ijk = (nijk).norm(); + double n_jkl = (njkl).norm(); + double dotpr = nijk.dot(njkl); + Eigen::Vector3d ji = j - i; + + double sign = (-1 * ji).dot(njkl) < 0 ? -1 : 1; + double phi = pi + sign * acos(dotpr / (n_ijk * n_jkl)); + double sinphi = sin(phi); + double tmp_energy = (1 / 2.0 * dihedral.V * (1 - cos(dihedral.n * dihedral.phi0) * cos(dihedral.n * dihedral.phi0))) * m_final_factor * m_dihedral_scaling; + // Dihedral(i, j, k, l, dihedral.V, dihedral.n, dihedral.phi0) * m_final_factor * m_dihedral_scaling; + if (isnan(tmp_energy)) + continue; + energy += tmp_energy; + if (m_calculate_gradient) { + Eigen::Vector3d kj = k - j; + Eigen::Vector3d kl = k - l; + double tmp = 0; + double dEdphi = (1 / 2.0 * dihedral.V * dihedral.n * (cos(dihedral.n * dihedral.phi0) * sin(dihedral.n * phi))); + if (isnan(dEdphi)) + continue; + + Eigen::Vector3d dEdi = dEdphi * kj.norm() / (nijk.norm() * nijk.norm()) * nijk; + Eigen::Vector3d dEdl = -1 * dEdphi * kj.norm() / (njkl.norm() * njkl.norm()) * njkl; + Eigen::Vector3d dEdj = -1 * dEdi + ((-1 * ji).dot(kj) / (kj.norm() * kj.norm()) * dEdi) - (kl.dot(kj) / (kj.norm() * kj.norm()) * dEdl); + Eigen::Vector3d dEdk = -1 * (dEdi + dEdj + dEdl); + + if (isnan(dEdi.sum()) || isnan(dEdj.sum()) || isnan(dEdk.sum()) || isnan(dEdl.sum())) + continue; + m_gradient(dihedral.i, 0) += dEdi(0); + m_gradient(dihedral.i, 1) += dEdi(1); + m_gradient(dihedral.i, 2) += dEdi(2); + + m_gradient(dihedral.l, 0) += dEdl(0); + m_gradient(dihedral.l, 1) += dEdl(1); + m_gradient(dihedral.l, 2) += dEdl(2); + + m_gradient(dihedral.j, 0) += dEdj(0); + m_gradient(dihedral.j, 1) += dEdj(1); + m_gradient(dihedral.j, 2) += dEdj(2); + + m_gradient(dihedral.k, 0) += dEdk(0); + m_gradient(dihedral.k, 1) += dEdk(1); + m_gradient(dihedral.k, 2) += dEdk(2); + } + } + return energy; +} + +double ForceFieldThread::CalculateUFFInversionContribution() +{ + double energy = 0; + + for (int index = 0; index < m_uff_inversions.size(); ++index) { + const auto& inversion = m_uff_inversions[index]; + + Eigen::Vector3d i = Position(inversion.i); + Eigen::Vector3d j = Position(inversion.j); + Eigen::Vector3d k = Position(inversion.k); + Eigen::Vector3d l = Position(inversion.l); + + Eigen::Vector3d ail = SubVector(i, l); + Eigen::Vector3d nijk = UFF::NormalVector(i, j, k); + + double cosY = (nijk.dot(ail) / ((nijk).norm() * (ail).norm())); + + double sinYSq = 1.0 - cosY * cosY; + double sinY = ((sinYSq > 0.0) ? sqrt(sinYSq) : 0.0); + double cos2Y = sinY * sinY - 1.0; + + double tmp_energy = (inversion.fc * (inversion.C0 + inversion.C1 * sinY + inversion.C2 * cos2Y)) * m_final_factor * m_inversion_scaling; + if (std::isnan(tmp_energy)) + continue; + energy += tmp_energy; + + // energy += Inversion(i, j, k, l, d_forceConstant, C0, C1, C2); + if (m_calculate_gradient) { + Eigen::Vector3d ji = (j - i); + Eigen::Vector3d jk = (k - i); + Eigen::Vector3d jl = (l - i); + + if (ji.norm() < 1e-5 || jk.norm() < 1e-5 || jl.norm() < 1e-5) + continue; + + double dji = ji.norm(); + double djk = jk.norm(); + double djl = jl.norm(); + ji /= dji; + jk /= djk; + jl /= djl; + + Eigen::Vector3d nijk = ji.cross(jk); + nijk /= nijk.norm(); + + double cosY = (nijk.dot(jl)); + double sinYSq = 1.0 - cosY * cosY; + double sinY = ((sinYSq > 0.0) ? sqrt(sinYSq) : 0.0); + double cosTheta = (ji.dot(jk)); + double sinThetaSq = std::max(1.0 - cosTheta * cosTheta, 1.0e-8); + double sinTheta = std::max(((sinThetaSq > 0.0) ? sqrt(sinThetaSq) : 0.0), 1.0e-8); + + double dEdY = -1 * (inversion.fc * (inversion.C1 * cosY - 4 * inversion.C2 * cosY * sinY)) * m_final_factor * m_inversion_scaling; + + Eigen::Vector3d p1 = ji.cross(jk); + Eigen::Vector3d p2 = jk.cross(jl); + Eigen::Vector3d p3 = jl.cross(ji); + + const double sin_dl = p1.dot(jl) / sinTheta; + + // the wilson angle: + const double dll = asin(sin_dl); + const double cos_dl = cos(dll); + + Eigen::Vector3d dYdl = (p1 / sinTheta - (jl * sin_dl)) / djl; + Eigen::Vector3d dYdi = ((p2 + (((-ji + jk * cosTheta) * sin_dl) / sinTheta)) / dji) / sinTheta; + Eigen::Vector3d dYdk = ((p3 + (((-jk + ji * cosTheta) * sin_dl) / sinTheta)) / djk) / sinTheta; + Eigen::Vector3d dYdj = -1 * (dYdi + dYdk + dYdl); + + m_gradient(inversion.i, 0) += dEdY * dYdj(0); + m_gradient(inversion.i, 1) += dEdY * dYdj(1); + m_gradient(inversion.i, 2) += dEdY * dYdj(2); + + m_gradient(inversion.j, 0) += dEdY * dYdi(0); + m_gradient(inversion.j, 1) += dEdY * dYdi(1); + m_gradient(inversion.j, 2) += dEdY * dYdi(2); + + m_gradient(inversion.k, 0) += dEdY * dYdk(0); + m_gradient(inversion.k, 1) += dEdY * dYdk(1); + m_gradient(inversion.k, 2) += dEdY * dYdk(2); + + m_gradient(inversion.l, 0) += dEdY * dYdl(0); + m_gradient(inversion.l, 1) += dEdY * dYdl(1); + m_gradient(inversion.l, 2) += dEdY * dYdl(2); + } + } + return energy; +} + +double ForceFieldThread::CalculateUFFvdWContribution() +{ + double energy = 0; + for (int index = 0; index < m_uff_vdWs.size(); ++index) { + const auto& vdw = m_uff_vdWs[index]; + Eigen::Vector3d i = Position(vdw.i); + Eigen::Vector3d j = Position(vdw.j); + double ij = (i - j).norm() * m_au; + double pow6 = pow((vdw.r0_ij / ij), 6); + + energy += vdw.C_ij * (-2 * pow6 * m_vdw_scaling + pow6 * pow6 * m_rep_scaling) * m_final_factor; + if (m_calculate_gradient) { + double diff = 12 * vdw.C_ij * (pow6 * m_vdw_scaling - pow6 * pow6 * m_rep_scaling) / (ij * ij) * m_final_factor; + m_gradient(vdw.i, 0) += diff * (i(0) - j(0)); + m_gradient(vdw.i, 1) += diff * (i(1) - j(1)); + m_gradient(vdw.i, 2) += diff * (i(2) - j(2)); + + m_gradient(vdw.j, 0) -= diff * (i(0) - j(0)); + m_gradient(vdw.j, 1) -= diff * (i(1) - j(1)); + m_gradient(vdw.j, 2) -= diff * (i(2) - j(2)); + } + } + return energy; +} + +double ForceFieldThread::CalculateQMDFFBondContribution() +{ + double factor = m_final_factor * m_bond_scaling; + double energy = 0.0; + Eigen::Vector3d dx = { m_d, 0, 0 }; + Eigen::Vector3d dy = { 0, m_d, 0 }; + Eigen::Vector3d dz = { 0, 0, m_d }; + + for (int index = 0; index < m_qmdff_bonds.size(); ++index) { + const auto& bond = m_qmdff_bonds[index]; + + Vector i = Position(bond.i); + Vector j = Position(bond.j); + + // Matrix derivate; + energy += QMDFF::LJStretchEnergy((i - j).norm(), bond.r0_ij, bond.fc, bond.exponent) * factor; + if (m_calculate_gradient) { + /*if (m_calc_gradient == 0) { + double diff = 0; + m_gradient.row(a) += diff * derivate.row(0); + m_gradient.row(b) += diff * derivate.row(1); + + } else if (m_calc_gradient == 1) {*/ + double d_x = (QMDFF::LJStretchEnergy((i + dx - j).norm(), bond.r0_ij, bond.fc, bond.exponent) - QMDFF::LJStretchEnergy((i - dx - j).norm(), bond.r0_ij, bond.fc, bond.exponent)) / (2 * m_d) * factor; + double d_y = (QMDFF::LJStretchEnergy((i + dy - j).norm(), bond.r0_ij, bond.fc, bond.exponent) - QMDFF::LJStretchEnergy((i - dy - j).norm(), bond.r0_ij, bond.fc, bond.exponent)) / (2 * m_d) * factor; + double d_z = (QMDFF::LJStretchEnergy((i + dz - j).norm(), bond.r0_ij, bond.fc, bond.exponent) - QMDFF::LJStretchEnergy((i - dz - j).norm(), bond.r0_ij, bond.fc, bond.exponent)) / (2 * m_d) * factor; + m_gradient(bond.i, 0) += d_x; + m_gradient(bond.i, 1) += d_y; + m_gradient(bond.i, 2) += d_z; + + m_gradient(bond.j, 0) -= d_x; + m_gradient(bond.j, 1) -= d_y; + m_gradient(bond.j, 2) -= d_x; + //} + } + } + + return energy; +} +/* +double ForceFieldThread::HarmonicBondStretching() +{ + double energy = 0; + + return energy; +} +*/ +double ForceFieldThread::CalculateQMDFFAngleContribution() +{ + double threshold = 1e-2; + double energy = 0.0; + Eigen::Vector3d dx = { m_d, 0, 0 }; + Eigen::Vector3d dy = { 0, m_d, 0 }; + Eigen::Vector3d dz = { 0, 0, m_d }; + for (int index = 0; index < m_qmdff_angles.size(); ++index) { + const auto& angle = m_qmdff_angles[index]; + /* + const int a = angle.a; + const int b = angle.b; + const int c = angle.c; +*/ + auto i = Position(angle.i); + auto j = Position(angle.j); + auto k = Position(angle.k); + Matrix derivate; + double costheta = AngleBending(i, j, k, derivate, m_calculate_gradient); + std::function angle_function; + if (std::abs(costheta - pi) < threshold) + + angle_function = [this](const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, double thetae, double fc, double r0_ij, double r0_ik) -> double { + double val = QMDFF::LinearAngleBend(i, j, k, thetae, fc, r0_ij, r0_ik); + if (std::isnan(val)) + return 0; + else + return val; + }; + else + angle_function = [this](const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, double thetae, double fc, double r0_ij, double r0_ik) -> double { + double val = QMDFF::AngleBend(i, j, k, thetae, fc, r0_ij, r0_ik); + if (std::isnan(val)) + return 0; + else + return val; + }; + + double e = angle_function(i, j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik); //(angle.kijk * (angle.C0 + angle.C1 * costheta + angle.C2 * (2 * costheta * costheta - 1))) * m_final_factor * m_angle_scaling; + + if (m_calculate_gradient) { + if (m_calc_gradient == 0) { + + } else if (m_calc_gradient == 1) { + m_gradient(angle.i, 0) += (angle_function((i + dx), j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function((i - dx), j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + m_gradient(angle.i, 1) += (angle_function((i + dy), j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function((i - dy), j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + m_gradient(angle.i, 2) += (angle_function((i + dz), j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function((i - dz), j, k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + + m_gradient(angle.j, 0) += (angle_function(i, (j + dx), k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function(i, (j - dx), k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + m_gradient(angle.j, 1) += (angle_function(i, (j + dy), k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function(i, (j - dy), k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + m_gradient(angle.j, 2) += (angle_function(i, (j + dz), k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function(i, (j - dz), k, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + + m_gradient(angle.k, 0) += (angle_function(i, j, k + dx, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function(i, j, k - dx, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + m_gradient(angle.k, 1) += (angle_function(i, j, k + dy, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function(i, j, k - dy, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + m_gradient(angle.k, 2) += (angle_function(i, j, k + dz, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik) - angle_function(i, j, k - dz, angle.theta0_ijk, angle.fc, angle.r0_ij, angle.r0_ik)) / (2 * m_d); + } + } + } + return energy; +} diff --git a/src/core/forcefieldthread.h b/src/core/forcefieldthread.h new file mode 100644 index 0000000..808d7c7 --- /dev/null +++ b/src/core/forcefieldthread.h @@ -0,0 +1,140 @@ +/* + * < Generic force field class for curcuma . > + * Copyright (C) 2024 Conrad Hübler + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + */ + +#pragma once + +#include "src/core/global.h" + +#include "hbonds.h" + +#include "external/CxxThreadPool/include/CxxThreadPool.h" + +#ifdef USE_D3 +#include "src/core/dftd3interface.h" +#endif + +#ifdef USE_D4 +#include "src/core/dftd4interface.h" +#endif + +#include "src/core/qmdff_par.h" +#include "src/core/uff_par.h" + +#include +#include +#include + +#include + +#include "json.hpp" +using json = nlohmann::json; + +struct Bond { + int type = 1; // 1 = UFF, 2 = QMDFF + int i = 0, j = 0, k = 0, distance = 0; + double fc = 0, exponent = 0, r0_ij = 0, r0_ik = 0; +}; + +struct Angle { + int type = 1; // 1 = UFF, 2 = QMDFF + int i = 0, j = 0, k = 0; + double fc = 0, r0_ij = 0, r0_ik = 0, theta0_ijk = 0; + double C0 = 0, C1 = 0, C2 = 0; +}; + +struct Dihedral { + int type = 1; // 1 = UFF, 2 = QMDFF + int i = 0, j = 0, k = 0, l = 0; + double V = 0, n = 0, phi0 = 0; +}; + +struct Inversion { + int type = 1; // 1 = UFF, 2 = QMDFF + int i = 0, j = 0, k = 0, l = 0; + double fc = 0, C0 = 0, C1 = 0, C2 = 0; +}; + +struct vdW { + int type = 1; // 1 = UFF, 2 = QMDFF + int i = 0, j = 0; + double C_ij = 0, r0_ij = 0; +}; + +struct EQ { + int type = 1; // 1 = UFF, 2 = QMDFF + int i = 0, j = 0; + double C_ij = 0, r0_ij = 0; +}; +class ForceFieldThread : public CxxThread { + +public: + ForceFieldThread(int thread, int threads); + virtual int execute() override; + + void addBond(const Bond& bonds); + void addAngle(const Angle& angles); + void addDihedral(const Dihedral& dihedrals); + void addInversion(const Inversion& inversions); + void addvdW(const vdW& vdWs); + void addEQ(const EQ& EQs); + + inline void setGeometry(const Matrix& geometry) { m_geometry = geometry; } + +private: + double CalculateUFFBondContribution(); + double CalculateUFFAngleContribution(); + double CalculateUFFDihedralContribution(); + double CalculateUFFInversionContribution(); + double CalculateUFFvdWContribution(); + + double CalculateQMDFFBondContribution(); + double CalculateQMDFFAngleContribution(); + double CalculateQMDFFDihedralContribution(); + + // double HarmonicBondStretching(); + + // double LJBondStretching(); + + std::function CalculateBondContribution; + std::function CalculateAngleContribution; + std::function CalculateTorsionContribution; + std::function CalculateInversionContribution; + std::function CalculateVdWContribution; + std::function CalculateEQContribution; + std::function CalculateHBondContribution; + + double m_energy = 0, m_bond_energy = 0.0, m_angle_energy = 0.0, m_dihedral_energy = 0.0, m_inversion_energy = 0.0, m_vdw_energy = 0.0, m_d4_energy = 0.0, m_d3_energy = 0.0, m_energy_h4 = 0.0, m_energy_hh = 0.0; + + double m_final_factor = 1; + double m_bond_scaling = 1, m_angle_scaling = 1, m_dihedral_scaling = 1, m_inversion_scaling = 1, m_vdw_scaling = 1, m_rep_scaling = 1; + double m_au = 1; + double m_d = 1e-3; + int m_calc_gradient = 1; + int m_thread = 0, m_threads = 0; + bool m_calculate_gradient = true; + + Matrix m_geometry, m_gradient; + + std::vector m_uff_bonds, m_qmdff_bonds; + std::vector m_uff_angles, m_qmdff_angles; + std::vector m_uff_dihedrals, m_qmdff_dihedrals; + std::vector m_uff_inversions, m_qmdff_inversions; + std::vector m_uff_vdWs; + std::vector m_qmdff_EQs; +}; diff --git a/src/core/qmdff.cpp b/src/core/qmdff.cpp index da74c45..747e52e 100644 --- a/src/core/qmdff.cpp +++ b/src/core/qmdff.cpp @@ -208,6 +208,7 @@ double QMDFFThread::CalculateAngleBending() auto atom_b = Position(b); auto atom_c = Position(c); Matrix derivate; + double costheta = AngleBending(atom_a, atom_b, atom_c, derivate, m_CalculateGradient); std::function angle_function; if (std::abs(costheta - pi) < threshold) diff --git a/src/main.cpp b/src/main.cpp index 002f55b..8c263a0 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -182,7 +182,7 @@ int main(int argc, char **argv) { if(strcmp(argv[1], "-rmsd") == 0) { - if (argc < 4) { + if (argc < 3) { std::cerr << "Please use curcuma for rmsd calcultion as follows\ncurcuma -rmsd A.xyz B.xyz" << std::endl; std::cerr << "Additonal arguments are:" << std::endl; std::cerr << "-reorder **** Force reordering of structure!" << std::endl; @@ -197,10 +197,15 @@ int main(int argc, char **argv) { Molecule mol2(argv[3]); // will only take first structure if (mol1.AtomCount() == 0 || mol2.AtomCount() == 0) { - std::cout << "At least one structure is empty:\n"; - std::cout << argv[2] << " " << mol1.AtomCount() << " atoms" << std::endl; - std::cout << argv[3] << " " << mol2.AtomCount() << " atoms" << std::endl; - exit(0); + FileIterator file(argv[2]); + file.Next(); + mol2 = file.Next(); + if (mol2.AtomCount() == 0) { + std::cout << "At least one structure is empty:\n"; + std::cout << argv[2] << " " << mol1.AtomCount() << " atoms" << std::endl; + std::cout << argv[3] << " " << mol2.AtomCount() << " atoms" << std::endl; + exit(0); + } } std::string reffile = argv[2];