diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f792d4..453bbc0 100644 --- 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 index 21b967a..1d2b49b 100644 --- a/src/capabilities/hessian.cpp +++ b/src/capabilities/hessian.cpp @@ -178,6 +178,7 @@ void Hessian::LoadControlJson() m_freq_cutoff = Json2KeyWord(m_defaults, "freq_cutoff"); m_hess = Json2KeyWord(m_defaults, "hess"); m_method = Json2KeyWord(m_defaults, "method"); + m_threads = Json2KeyWord(m_defaults, "threads"); } void Hessian::setMolecule(const Molecule& molecule) diff --git a/src/capabilities/hessian.h b/src/capabilities/hessian.h index 08b85b7..8023da3 100644 --- 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 index 93dd655..b2aa320 100644 --- 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 index 43a1c34..8c67584 100644 --- 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 index 9206dac..6616d22 100644 --- 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 new file mode 100644 index 0000000..4b55e14 --- /dev/null +++ b/src/core/forcefield.cpp @@ -0,0 +1,193 @@ +/* + * < 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 . + * + */ + +#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 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>& 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 new file mode 100644 index 0000000..42d221a --- /dev/null +++ b/src/core/forcefield.h @@ -0,0 +1,100 @@ +/* + * < 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 . + * + */ + +#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 "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: + ForceField(const json& controller); + ~ForceField(); + + void UpdateGeometry(const Matrix& geometry); + inline void UpdateGeometry(const double* coord); + inline void UpdateGeometry(const std::vector>& 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; +}; diff --git a/src/core/forcefield_terms/qmdff_terms.h b/src/core/forcefield_terms/qmdff_terms.h new file mode 100644 index 0000000..1f0d265 --- /dev/null +++ b/src/core/forcefield_terms/qmdff_terms.h @@ -0,0 +1,77 @@ +/* + * + * 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 . + * + */ + +#pragma once + +#include + +namespace QMDFFTerms { + +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 (isnan(energy)) + return 0; + else + return energy; +} + +double AngleDamping(double r_ij, double r_ik, double r0_ij, double r0_ik) +{ + double kdamp = 1; + double ratio_ij = r_ij / r0_ij; + double ratio_ik = r_ik / r0_ik; + + double f_ij_inv = 1 + kdamp * pow(ratio_ij, 4); + double f_ik_inv = 1 + kdamp * pow(ratio_ik, 4); + double finv = f_ij_inv * f_ik_inv; + return 1 / finv; +} + +double AngleBend(const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, double thetae, 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 costhetae = cos(thetae); + double energy = (fc * damp * (costhetae - costheta) * (costhetae - costheta)); + if (isnan(energy)) + return 0; + else + return energy; +} + +double LinearAngleBend(const Eigen::Vector3d& i, const Eigen::Vector3d& j, const Eigen::Vector3d& k, double thetae, double fc, double r0_ij, double r0_ik) +{ + if (fc < 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 * (thetae - theta) * (thetae - theta)); + if (isnan(energy)) + return 0; + else + return energy; +} + +} diff --git a/src/core/forcefieldfunctions.h b/src/core/forcefieldfunctions.h new file mode 100644 index 0000000..23f1795 --- /dev/null +++ b/src/core/forcefieldfunctions.h @@ -0,0 +1,4 @@ +#ifndef FORCEFIELDFUNCTIONS_H +#define FORCEFIELDFUNCTIONS_H + +#endif // FORCEFIELDFUNCTIONS_H diff --git a/src/core/qmdff.cpp b/src/core/qmdff.cpp index 5ae8f85..da74c45 100644 --- a/src/core/qmdff.cpp +++ b/src/core/qmdff.cpp @@ -132,13 +132,17 @@ double QMDFFThread::CalculateStretchEnergy() double za = (*m_geometry)(a, 2) * m_au; double zb = (*m_geometry)(b, 2) * m_au; - m_gradient(a, 0) += (StretchEnergy(Distance(xa + m_d, xb, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa - m_d, xb, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); - m_gradient(a, 1) += (StretchEnergy(Distance(xa, xb, ya + m_d, yb, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya - m_d, yb, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); - m_gradient(a, 2) += (StretchEnergy(Distance(xa, xb, ya, yb, za + m_d, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya, yb, za - m_d, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); - m_gradient(b, 0) += (StretchEnergy(Distance(xa, xb + m_d, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb - m_d, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); - m_gradient(b, 1) += (StretchEnergy(Distance(xa, xb, ya, yb + m_d, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya, yb - m_d, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); - m_gradient(b, 2) += (StretchEnergy(Distance(xa, xb, ya, yb, za, zb + m_d), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya, yb, za, zb - m_d), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); + double d_x = (StretchEnergy(Distance(xa + m_d, xb, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa - m_d, xb, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); + double d_y = (StretchEnergy(Distance(xa, xb, ya + m_d, yb, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya - m_d, yb, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); + double d_z = (StretchEnergy(Distance(xa, xb, ya, yb, za + m_d, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya, yb, za - m_d, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); + m_gradient(a, 0) += d_x; + m_gradient(a, 1) += d_y; + m_gradient(a, 2) += d_z; + + m_gradient(b, 0) -= d_x; //(StretchEnergy(Distance(xa, xb + m_d, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb - m_d, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); + m_gradient(b, 1) -= d_y; //(StretchEnergy(Distance(xa, xb, ya, yb + m_d, za, zb), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya, yb - m_d, za, zb), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); + m_gradient(b, 2) -= d_z; //(StretchEnergy(Distance(xa, xb, ya, yb, za, zb + m_d), bond.reAB, bond.kAB, bond.exponA) - StretchEnergy(Distance(xa, xb, ya, yb, za, zb - m_d), bond.reAB, bond.kAB, bond.exponA)) / (2 * m_d); } } }