more on forcefield

Signed-off-by: Conrad Hübler <Conrad.Huebler@gmx.net>
conradhuebler · Feb 18, 2024 · c3f6aa5 · c3f6aa5
1 parent 9fff262
commit c3f6aa5
Show file tree

Hide file tree

Showing 11 changed files with 396 additions and 23 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -164,6 +164,10 @@ 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/forcefield_terms/qmdff_terms.h
         src/tools/formats.h
         src/tools/geometry.h
         src/tools/general.h
@@ -172,6 +176,7 @@ set(curcuma_core_SRC
 
     add_executable(curcuma
             src/main.cpp
+
             )
 
 if(C17)

diff --git a/src/capabilities/hessian.cpp b/src/capabilities/hessian.cpp
@@ -178,6 +178,7 @@ void Hessian::LoadControlJson()
     m_freq_cutoff = Json2KeyWord<double>(m_defaults, "freq_cutoff");
     m_hess = Json2KeyWord<int>(m_defaults, "hess");
     m_method = Json2KeyWord<std::string>(m_defaults, "method");
+    m_threads = Json2KeyWord<int>(m_defaults, "threads");
 }
 
 void Hessian::setMolecule(const Molecule& molecule)

diff --git a/src/capabilities/hessian.h b/src/capabilities/hessian.h
@@ -41,7 +41,8 @@ static const json HessianJson = {
     { "thermo", 298.15 },
     { "freq_cutoff", 50 },
     { "hess", 1 },
-    { "method", "uff" }
+    { "method", "uff" },
+    { "threads", 1 }
 };
 
 class HessianThread : public CxxThread {

diff --git a/src/capabilities/qmdfffit.cpp b/src/capabilities/qmdfffit.cpp
@@ -46,6 +46,7 @@ void QMDFFFit::start()
     std::cout << "Parametrising QMDFF (see S. Grimmme, J. Chem. Theory Comput. 2014, 10, 10, 4497–4514 [10.1021/ct500573f]) for the original publication!" << std::endl;
 
     std::cout << "Starting with the hessian ..." << std::endl;
+    std::cout << m_defaults << m_controller << std::endl;
     Hessian hessian(m_defaults);
     hessian.setMolecule(m_molecule);
     m_atom_types = m_molecule.Atoms();
@@ -83,7 +84,7 @@ void QMDFFFit::start()
                 const_hessian_matrix(j, i) = 0;
             }
         }
-    std::cout << const_hessian_matrix << std::endl;
+    // std::cout << const_hessian_matrix << std::endl;
     int counter = 1;
     for (int start = 0; start < 10 && counter; ++start) {
         parameter["bonds"] = Bonds();

diff --git a/src/capabilities/simplemd.cpp b/src/capabilities/simplemd.cpp
@@ -1313,21 +1313,8 @@ bool SimpleMD::WriteGeometry()
         geometry(i, 2) = m_current_geometry[3 * i + 2];
     }
     TriggerWriteRestart();
-    // int f1 = m_molecule.GetFragments().size();
     m_molecule.setGeometry(geometry);
-    auto m = m_molecule.DistanceMatrix();
 
-    // int f2 = m_molecule.GetFragments().size();
-    //   std::cout << f1 << " ... " << f2 << std::endl;
-    // m_prev_index = std::abs(f2 - f1);
-    /*
-    int difference = (m.second - m_topo_initial).cwiseAbs().sum();
-
-    if (difference > m_max_top_diff) {
-        std::cout << "*** topology changed " << difference << " ***" << std::endl;
-        result = false;
-    }
-    */
     if (m_writeXYZ) {
         m_molecule.setEnergy(m_Epot);
         m_molecule.setName(std::to_string(m_currentStep));

diff --git a/src/core/eigen_uff.h b/src/core/eigen_uff.h
@@ -118,8 +118,8 @@ class UFFThread : public CxxThread
         auto atom_1 = m_geometry->row(atom2);
         auto atom_2 = m_geometry->row(atom3);
 
-        auto vec_1 = atom_0 - atom_1; //{ atom_0[0] - atom_1[0], atom_0[1] - atom_1[1], atom_0[2] - atom_1[2] };
-        auto vec_2 = atom_0 - atom_2; //{ atom_0[0] - atom_2[0], atom_0[1] - atom_2[1], atom_0[2] - atom_2[2] };
+        auto vec_1 = atom_0 - atom_1;
+        auto vec_2 = atom_0 - atom_2;
 
         return acos(vec_1.dot(vec_2) / sqrt(vec_1.dot(vec_1)) * sqrt(vec_2.dot(vec_2))) * 360 / 2.0 / pi;
     }

diff --git a/src/core/forcefield.cpp b/src/core/forcefield.cpp
@@ -0,0 +1,193 @@
+/*
+ * < Generic force field class for curcuma . >
+ * Copyright (C) 2022 - 2023 Conrad Hübler <Conrad.Huebler@gmx.net>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "src/core/forcefieldderivaties.h"
+#include "src/core/qmdff_par.h"
+#include "src/core/uff_par.h"
+
+#include "forcefield.h"
+
+ForceFieldThread::ForceFieldThread()
+{
+}
+
+double ForceFieldThread::LJBondStretching()
+{
+    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;
+            //}
+        }
+    }
+
+    return energy;
+}
+
+double ForceFieldThread::HarmonicBondStretching()
+{
+    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);
+        }
+    }
+
+    return energy;
+}
+
+double QMDFFThread::CalculateAngleBending()
+{
+    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<double(const Eigen::Vector3d&, const Eigen::Vector3d&, const Eigen::Vector3d&, double, double, double, double)> 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);
+            }
+        }
+    }
+    return energy;
+}
+
+ForceField::ForceField(const json& controller)
+{
+}
+
+void ForceField::UpdateGeometry(const Matrix& geometry)
+{
+    m_geometry = geometry;
+}
+
+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];
+    }
+}
+
+void ForceField::UpdateGeometry(const std::vector<std::array<double, 3>>& geometry);
+{
+#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];
+    }
+}
+
+void ForceField::setParameter(const json& parameter)
+{
+}
diff --git a/src/core/forcefield.h b/src/core/forcefield.h
@@ -0,0 +1,100 @@
+/*
+ * < Generic force field class for curcuma . >
+ * Copyright (C) 2022 - 2023 Conrad Hübler <Conrad.Huebler@gmx.net>
+ *
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#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 <set>
+#include <vector>
+
+#include <Eigen/Dense>
+
+#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<Bond> m_bonds;
+};
+
+class ForceField {
+
+public:
+    ForceField(const json& controller);
+    ~ForceField();
+
+    void UpdateGeometry(const Matrix& geometry);
+    inline void UpdateGeometry(const double* coord);
+    inline void UpdateGeometry(const std::vector<std::array<double, 3>>& geometry);
+
+    double Calculate(bool gradient = true, bool verbose = false);
+
+    Matrix Gradient() const { return m_gradient; }
+
+    void setParameter(const json& parameter);
+
+private:
+    Matrix m_geometry, m_gradient;
+    int m_natoms = 0;
+};