From 9fff262d5c81c1709ed6f19a09386d2eec5404bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Conrad=20H=C3=BCbler?= Date: Sun, 18 Feb 2024 13:57:44 +0100 Subject: [PATCH] first qmdff protoyp, without repulsion and torsion MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Conrad Hübler --- CMakeLists.txt | 8 +- external/simple-dftd3 | 2 +- src/capabilities/hessian.cpp | 20 +- src/capabilities/hessian.h | 9 +- src/capabilities/optimiser/LevMarQMDFFFit.h | 59 ++- src/capabilities/qmdfffit.cpp | 160 +++++-- src/capabilities/qmdfffit.h | 2 +- src/core/dftd3interface.cpp | 97 ++-- src/core/dftd3interface.h | 3 + src/core/dftd4interface.cpp | 1 + src/core/eigen_uff.cpp | 1 - src/core/energycalculator.cpp | 8 +- src/core/qmdff.cpp | 472 +++++++------------- src/core/qmdff.h | 51 +-- src/core/qmdff_par.h | 46 ++ 15 files changed, 436 insertions(+), 503 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9938554..9f792d4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -272,12 +272,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 ) @@ -288,8 +288,8 @@ endif() if(USE_D4) if(NOT DEFINED D4LIBS_DIR AND NOT DEFINED D4INCLUDE_DIR) - find_package(LAPACKE REQUIRED) - find_package(CBLAS CONFIG REQUIRED) + #find_package(LAPACKE REQUIRED) + #find_package(CBLAS CONFIG REQUIRED) else() include_directories(${INCLUDE_DIR}) endif(NOT DEFINED D4LIBS_DIR AND NOT DEFINED D4INCLUDE_DIR) diff --git a/external/simple-dftd3 b/external/simple-dftd3 index d5b72de..3a50cbd 160000 --- a/external/simple-dftd3 +++ b/external/simple-dftd3 @@ -1 +1 @@ -Subproject commit d5b72de4bb0fe892410180a67c23556db6726e67 +Subproject commit 3a50cbdf8566b3caa8760499ae21746ec9a5e535 diff --git a/src/capabilities/hessian.cpp b/src/capabilities/hessian.cpp index b559bbb..21b967a 100644 --- a/src/capabilities/hessian.cpp +++ b/src/capabilities/hessian.cpp @@ -25,9 +25,8 @@ #include "hessian.h" -HessianThread::HessianThread(const std::string& method, const json& controller, int i, int j, int xi, int xj, bool fullnumerical) - : m_method(method) - , m_controller(controller) +HessianThread::HessianThread(const json& controller, int i, int j, int xi, int xj, bool fullnumerical) + : m_controller(controller) , m_i(i) , m_j(j) , m_xi(xi) @@ -35,6 +34,7 @@ HessianThread::HessianThread(const std::string& method, const json& controller, , m_fullnumerical(fullnumerical) { setAutoDelete(true); + m_method = m_controller["method"]; if (m_fullnumerical) m_schema = [this, i, j, xi, xj]() { @@ -105,6 +105,7 @@ void HessianThread::Seminumerical() m_geom_im_jp = m_molecule.Coords(); m_geom_ip_jp[m_i][m_xi] += m_d; + // std::cout << m_controller << std::endl; EnergyCalculator energy(m_method, m_controller); energy.setMolecule(m_molecule); @@ -122,13 +123,14 @@ void HessianThread::Seminumerical() m_gradient = (gradientp - gradientm) / (2 * m_d) / au / au; } -Hessian::Hessian(const std::string& method, const json& controller, int threads, bool silent) +Hessian::Hessian(const std::string& method, const json& controller, bool silent) : CurcumaMethod(HessianJson, controller, silent) , m_method(method) , m_controller(controller) - , m_threads(threads) + { UpdateController(controller); + m_threads = m_controller["threads"]; /* Yeah, thats not really correct, but it works a bit */ if (m_method.compare("gfnff") == 0) { m_scale_functions = [](double val) -> double { @@ -228,14 +230,14 @@ void Hessian::start() LoadHessian(m_read_file); } - auto eigenvalues = ConvertHessian(m_hessian); + m_frequencies = ConvertHessian(m_hessian); // PrintVibrationsPlain(eigenvalues); auto hessian2 = ProjectHessian(m_hessian); auto projected = ConvertHessian(hessian2); if (!m_silent) - PrintVibrations(eigenvalues, projected); + PrintVibrations(m_frequencies, projected); } Matrix Hessian::ProjectHessian(const Matrix& hessian) @@ -356,7 +358,7 @@ void Hessian::CalculateHessianNumerical() for (int j = 0; j < m_molecule.AtomCount(); ++j) { for (int xi = 0; xi < 3; ++xi) for (int xj = 0; xj < 3; ++xj) { - HessianThread* thread = new HessianThread(m_method, m_controller, i, j, xi, xj, true); + HessianThread* thread = new HessianThread(m_controller, i, j, xi, xj, true); thread->setMolecule(m_molecule); thread->setParameter(m_parameter); pool->addThread(thread); @@ -386,7 +388,7 @@ void Hessian::CalculateHessianSemiNumerical() for (int i = 0; i < m_molecule.AtomCount(); ++i) { for (int xi = 0; xi < 3; ++xi) { - HessianThread* thread = new HessianThread(m_method, m_controller, i, 0, xi, 0, false); + HessianThread* thread = new HessianThread(m_controller, i, 0, xi, 0, false); thread->setMolecule(m_molecule); thread->setParameter(m_parameter); pool->addThread(thread); diff --git a/src/capabilities/hessian.h b/src/capabilities/hessian.h index fb6805b..08b85b7 100644 --- a/src/capabilities/hessian.h +++ b/src/capabilities/hessian.h @@ -46,7 +46,7 @@ static const json HessianJson = { class HessianThread : public CxxThread { public: - HessianThread(const std::string& method, const json& controller, int i, int j, int xi, int xj, bool fullnumerical = true); + HessianThread(const json& controller, int i, int j, int xi, int xj, bool fullnumerical = true); ~HessianThread(); void setMolecule(const Molecule& molecule); @@ -65,7 +65,7 @@ class HessianThread : public CxxThread { void Seminumerical(); std::function m_schema; - EnergyCalculator* m_calculator; + // EnergyCalculator* m_calculator; std::string m_method; json m_controller, m_parameter; Molecule m_molecule; @@ -79,7 +79,7 @@ class HessianThread : public CxxThread { class Hessian : public CurcumaMethod { public: - Hessian(const std::string& method, const json& controller, int threads, bool silent = true); + Hessian(const std::string& method, const json& controller, bool silent = true); Hessian(const json& controller, bool silent = true); void setMolecule(const Molecule& molecule); @@ -89,6 +89,7 @@ class Hessian : public CurcumaMethod { m_hess_calc = type; start(); } + void PrintVibrations() { PrintVibrationsPlain(m_frequencies); } void PrintVibrations(Vector& eigenvalues, const Vector& projected); void PrintVibrationsPlain(const Vector& eigenvalues); @@ -97,6 +98,7 @@ class Hessian : public CurcumaMethod { Matrix getHessian() const { return m_hessian; } void setParameter(const json& parameter) { m_parameter = parameter; } + Vector Frequencies() { return m_frequencies; } private: /* Lets have this for all modules */ @@ -126,6 +128,7 @@ class Hessian : public CurcumaMethod { Vector ConvertHessian(Matrix& hessian); Matrix m_eigen_geometry, m_eigen_gradient, m_hessian; + Vector m_frequencies; Molecule m_molecule; std::string m_method; json m_controller, m_parameter; diff --git a/src/capabilities/optimiser/LevMarQMDFFFit.h b/src/capabilities/optimiser/LevMarQMDFFFit.h index f736a51..e0a4f76 100644 --- a/src/capabilities/optimiser/LevMarQMDFFFit.h +++ b/src/capabilities/optimiser/LevMarQMDFFFit.h @@ -1,6 +1,6 @@ /* * - * Copyright (C) 2023 Conrad Hübler + * Copyright (C) 2023 - 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 @@ -34,6 +34,9 @@ #include "src/tools/geometry.h" +#include "json.hpp" +using json = nlohmann::json; + template struct FCFunctor { @@ -68,11 +71,7 @@ struct MyFCFunctor : FCFunctor { inline ~MyFCFunctor() {} inline int operator()(const Eigen::VectorXd& fc, Eigen::VectorXd& fvec) const { - // std::cout << fc.transpose() << std::endl; - json hc = HessianJson; - hc["method"] = "qmdff"; - hc["threads"] = m_threads; - Hessian he2(hc, true); + Hessian he2(m_controller, true); json bonds = m_parameter["bonds"]; json angles = m_parameter["angles"]; @@ -94,43 +93,59 @@ struct MyFCFunctor : FCFunctor { Matrix hessian = he2.getHessian(); index = 0; + double diff = 0; for (int i = 0; i < m_hessian.rows(); ++i) { - for (int j = 0; j < m_hessian.cols(); ++j) { - fvec(index++) = (m_hessian(i, j) - hessian(i, j)) * (m_hessian(i, j) - hessian(i, j)); // + PseudoFF::DistancePenalty(m_host->Atom(i), guest.Atom(j)); + for (int j = i; j < m_hessian.cols(); ++j) { + index++; + if (index >= fvec.size()) + break; + fvec(index) = (m_hessian(i, j) - (hessian(i, j) + m_const_hessian(i, j))); + diff += (m_hessian(i, j) - (hessian(i, j) + m_const_hessian(i, j))); } } - return 0; } + + void Controller(const json& controller) + { + m_controller = MergeJson(HessianJson, m_controller); + m_controller["method"] = "qmdff"; + } + int no_parameter; int no_points; - int m_threads = 1; int inputs() const { return no_parameter; } int values() const { return no_points; } - Matrix m_hessian; - json m_parameter; + Matrix m_hessian, m_const_hessian; + json m_parameter, m_controller; Molecule m_molecule; }; struct MyFunctorNumericalDiff : Eigen::NumericalDiff { }; -inline Vector OptimiseFC(const Molecule& molecule, const Matrix& hessian, const Vector& fc, const json& parameters, int threads) +inline Vector OptimiseFC(const Molecule& molecule, const Matrix& hessian, const Matrix& const_hessian, const Vector& fc, const json& parameters, const json& controller) { - - // Vector parameter = PositionPair2Vector(anchor, rotation); Vector parameter = fc; - MyFCFunctor functor(fc.size(), hessian.cols() * hessian.rows()); + MyFCFunctor functor(fc.size(), hessian.cols() * hessian.rows() / 2 + hessian.cols() / 2); functor.m_hessian = hessian; + functor.m_const_hessian = const_hessian; functor.m_molecule = molecule; functor.m_parameter = parameters; - functor.m_threads = threads; + functor.m_controller = controller; + + std::cout << controller << std::endl; Eigen::NumericalDiff numDiff(functor); Eigen::LevenbergMarquardt> lm(numDiff); int iter = 0; - + /* + std::cout << "Target hessian " << std::endl + << hessian << std::endl; + std::cout << "Const part of qmdff hessian " << std::endl + << const_hessian << std::endl; + */ /* lm.parameters.factor = config["LevMar_Factor"].toInt(); //step bound for the diagonal shift, is this related to damping parameter, lambda? lm.parameters.maxfev = config["LevMar_MaxFEv"].toDouble(); //max number of function evaluations @@ -141,21 +156,17 @@ inline Vector OptimiseFC(const Molecule& molecule, const Matrix& hessian, const */ Eigen::LevenbergMarquardtSpace::Status status = lm.minimizeInit(parameter); - // std::cout << parameter.transpose() << std::endl; - int MaxIter = 300; + int MaxIter = 30; Vector old_param = parameter; - // std::cout << parameter.transpose() << std::endl; for (; iter < MaxIter /*&& ((qAbs(error_0 - error_2) > ErrorConvergence) || norm > DeltaParameter)*/; ++iter) { std::cout << "Step " << iter << " of maximal " << MaxIter << std::endl; status = lm.minimizeOneStep(parameter); - + std::cout << "Norm " << (old_param - parameter).norm() << std::endl; if ((old_param - parameter).norm() < 1e-8) break; std::cout << parameter.transpose() << std::endl; old_param = parameter; } - // std::cout << parameter.transpose() << std::endl; - // std::cout << "took " << iter << " optimisation steps." << std::endl; return old_param; } diff --git a/src/capabilities/qmdfffit.cpp b/src/capabilities/qmdfffit.cpp index fcd33b1..93dd655 100644 --- a/src/capabilities/qmdfffit.cpp +++ b/src/capabilities/qmdfffit.cpp @@ -1,6 +1,6 @@ /* * - * Copyright (C) 2023 Conrad Hübler + * Copyright (C) 2023 - 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 @@ -42,10 +42,11 @@ void QMDFFFit::LoadControlJson() 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; - Hessian hessian(m_method); + Hessian hessian(m_defaults); hessian.setMolecule(m_molecule); m_atom_types = m_molecule.Atoms(); m_geometry = m_molecule.getGeometry(); @@ -54,23 +55,55 @@ void QMDFFFit::start() // std::cout << m_hessian << std::endl; Initialise(); json parameter; - for (int start = 0; start < 2; ++start) { + parameter["bonds"] = Bonds(); + parameter["angles"] = Angles(); + std::ofstream parameterfile_init("qmdff_init_param.json"); + parameterfile_init << parameter; + json qmdff_init = MergeJson(m_controller, QMDFFFitJson); + qmdff_init["threads"] = m_threads; + qmdff_init["variable"] = false; + qmdff_init["method"] = "qmdff"; + EnergyCalculator calculator("qmdff", qmdff_init); + calculator.setMolecule(m_molecule); + calculator.setParameter(parameter); + calculator.CalculateEnergy(false, false); + Hessian const_hessian("qmdff", qmdff_init); + /* + parameter["bonds"] = bonds; + parameter["angles"] = angles; + */ + const_hessian.setMolecule(m_molecule); + const_hessian.setParameter(parameter); + const_hessian.start(); + Matrix const_hessian_matrix = const_hessian.getHessian(); + for (int i = 0; i < const_hessian_matrix.cols(); ++i) + for (int j = i; j < const_hessian_matrix.cols(); ++j) { + if (std::isnan(const_hessian_matrix(i, j))) { + const_hessian_matrix(i, j) = 0; + const_hessian_matrix(j, i) = 0; + } + } + std::cout << const_hessian_matrix << std::endl; + int counter = 1; + for (int start = 0; start < 10 && counter; ++start) { parameter["bonds"] = Bonds(); parameter["angles"] = Angles(); m_fc_parameter.resize(m_qmdffbonds.size() + m_qmdffangle.size()); int i = 0; for (const auto& b : m_qmdffbonds) { - m_fc_parameter(i) = b.kAB; + m_fc_parameter(i) = std::abs(b.kAB); ++i; } for (const auto& b : m_qmdffangle) { - m_fc_parameter(i) = b.kabc; + m_fc_parameter(i) = std::abs(b.kabc); ++i; } - // std::cout << m_fc_parameter << std::endl; - Vector vec = OptimiseFC(m_molecule, m_hessian, m_fc_parameter, parameter, m_threads); + qmdff_init["variable"] = true; + qmdff_init["const"] = false; + Vector vec = OptimiseFC(m_molecule, m_hessian, const_hessian_matrix, m_fc_parameter, parameter, qmdff_init); + std::cout << vec << std::endl; json bonds = parameter["bonds"]; json angles = parameter["angles"]; int index = 0; @@ -87,27 +120,45 @@ void QMDFFFit::start() json hc = HessianJson; hc["method"] = "qmdff"; hc["threads"] = m_threads; - Hessian he2(hc, false); - parameter["bonds"] = bonds; - parameter["angles"] = angles; - he2.setMolecule(m_molecule); - he2.setParameter(parameter); - he2.start(); auto cache = m_qmdffbonds; m_qmdffbonds.clear(); + counter = 0; for (auto c : cache) { - if (c.kAB > 0) + if ((c.kAB > 0 && c.distance == 1) || c.distance == 0) m_qmdffbonds.push_back(c); + else + counter++; } auto cache2 = m_qmdffangle; m_qmdffangle.clear(); for (auto c : cache2) { if (c.kabc > 0) m_qmdffangle.push_back(c); + else + counter++; } + std::cout << "Skipping " << counter << " force constants " << std::endl; // std::cout << he2.getHessian() << std::endl; } + + std::cout << "Eigenvectors" << std::endl; + std::cout << hessian.Frequencies().transpose() << std::endl; + qmdff_init["variable"] = true; + qmdff_init["const"] = true; + + Hessian he2(qmdff_init, false); + + parameter["bonds"] = Bonds(); + parameter["angles"] = Angles(); + he2.setMolecule(m_molecule); + he2.setParameter(parameter); + he2.start(); + std::cout << he2.Frequencies().transpose() << std::endl; + + std::cout << parameter << std::endl; + std::ofstream parameterfile("qmdff_param.json"); + parameterfile << parameter; } bool QMDFFFit::Initialise() @@ -168,8 +219,6 @@ bool QMDFFFit::Initialise() nonbonds.clean(); setvdWs(ignored_vdw); - m_h4correction.allocate(m_atom_types.size()); - */ return true; } @@ -190,8 +239,8 @@ void QMDFFFit::setBonds(const TContainer& bonds, std::vector>& ign double za = (m_geometry)(b.a, 2) * m_au; double zb = (m_geometry)(b.b, 2) * m_au; - b.reAB = Topology::Distance(xa, xb, ya, yb, za, zb); // BondRestLength(b.i, b.j, bond_order); - b.kAB = 10; // 0.5 * m_bond_force * cZi * cZj / (b.r0 * b.r0 * b.r0); + b.reAB = Topology::Distance(xa, xb, ya, yb, za, zb); + b.kAB = 1; double kaA = ka(m_atom_types[b.a]); double kaB = ka(m_atom_types[b.b]); double dEN = Elements::PaulingEN[m_atom_types[b.a]] - Elements::PaulingEN[m_atom_types[b.b]]; @@ -249,40 +298,61 @@ void QMDFFFit::setAngles(const TContainer& angles, const std::vector - * Copyright (C) 2023 Conrad Hübler + * Copyright (C) 2023 - 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 diff --git a/src/core/dftd3interface.cpp b/src/core/dftd3interface.cpp index 7ed22a5..0e94946 100644 --- a/src/core/dftd3interface.cpp +++ b/src/core/dftd3interface.cpp @@ -46,7 +46,7 @@ DFTD3Interface::DFTD3Interface(const json& controller) m_atm = parameter["d_atm"]; m_damping = parameter["d_damping"]; m_functional = parameter["d_func"]; - + CreateParameter(); m_error = dftd3_new_error(); } @@ -55,6 +55,8 @@ DFTD3Interface::~DFTD3Interface() dftd3_delete_model(&m_disp); dftd3_delete_structure(&m_mol); dftd3_delete_error(&m_error); + dftd3_delete_param(&m_param); + delete m_coord; delete m_attyp; } @@ -63,32 +65,39 @@ void DFTD3Interface::PrintParameter() const { std::cout << m_d3_s6 << " " << m_d3_s8 << " " << m_d3_s9 << " " << m_d3_a1 << " " << m_d3_a2 << " " << m_d3_alp << std::endl; } -/* -void DFTD3Interface::LoadParameter() + +void DFTD3Interface::CreateParameter() { - dftd3_delete_param(&m_param); - char *cstr = new char[m_functional.length() + 1]; - strcpy(cstr, m_functional.c_str()); - if(m_damping.compare("bj") == 0) - { - m_param = dftd3_load_zero_damping(m_error, cstr, m_atm); - }else if(m_damping.compare("zero") == 0) - { - m_param = dftd3_load_rational_damping(m_error, cstr, m_atm); - }else if(m_damping.compare("bjm") == 0) - { - m_param = dftd3_load_mrational_damping(m_error, cstr, m_atm); - }else if(m_damping.compare("zerom") == 0) - { - m_param = dftd3_load_mzero_damping(m_error, cstr, m_atm); - }else if(m_damping.compare("op") == 0) - { - m_param = dftd3_load_optimizedpower_damping(m_error, cstr, m_atm); + if (m_d3_a1 > 1e-8 || m_d3_a2 > 1e-8 || m_d3_s6 > 1e-8 || m_d3_s8 > 1e-8 || m_d3_s9 > 1e-8) { + if (m_damping.compare("bj")) { + m_param = dftd3_new_rational_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp); + } else if (m_damping.compare("zero")) { + m_param = dftd3_new_zero_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp); + } else if (m_damping.compare("bjm")) { + m_param = dftd3_new_mrational_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp); + } else if (m_damping.compare("zerom")) { + m_param = dftd3_new_mzero_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp, m_bet); + } else if (m_damping.compare("op")) { + m_param = dftd3_new_optimizedpower_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp, m_bet); + } + } else { + char* cstr = new char[m_functional.length() + 1]; + strcpy(cstr, m_functional.c_str()); + if (m_damping.compare("bj") == 0) { + m_param = dftd3_load_rational_damping(m_error, cstr, m_atm); + } else if (m_damping.compare("zero") == 0) { + m_param = dftd3_load_zero_damping(m_error, cstr, m_atm); + } else if (m_damping.compare("bjm") == 0) { + m_param = dftd3_load_mrational_damping(m_error, cstr, m_atm); + } else if (m_damping.compare("zerom") == 0) { + m_param = dftd3_load_mzero_damping(m_error, cstr, m_atm); + } else if (m_damping.compare("op") == 0) { + m_param = dftd3_load_optimizedpower_damping(m_error, cstr, m_atm); + } + delete[] cstr; } - delete [] cstr; - PrintParameter(); } -*/ + void DFTD3Interface::UpdateParameters(const json& controller) { json parameter = MergeJson(DFTD3Settings, controller); @@ -100,6 +109,7 @@ void DFTD3Interface::UpdateParameters(const json& controller) m_d3_s8 = parameter["d_s8"]; m_d3_s9 = parameter["d_s9"]; + CreateParameter(); PrintParameter(); } @@ -109,9 +119,9 @@ bool DFTD3Interface::InitialiseMolecule(const std::vector& atomtypes) m_attyp = new int[atomtypes.size()]; m_coord = new double[3 * atomtypes.size()]; for (int i = 0; i < atomtypes.size(); ++i) { - m_coord[3 * i + 0] = 0 / au; - m_coord[3 * i + 1] = 0 / au; - m_coord[3 * i + 2] = 0 / au; + m_coord[3 * i + 0] = (3 * i + 0) / au; + m_coord[3 * i + 1] = (3 * i + 1) / au; + m_coord[3 * i + 2] = (3 * i + 2) / au; m_attyp[i] = atomtypes[i]; } @@ -132,40 +142,9 @@ double DFTD3Interface::DFTD3Calculation(double* grad) { double energy = 0; double sigma[9]; - dftd3_param param; - if (m_d3_a1 > 1e-8 || m_d3_a2 > 1e-8 || m_d3_s6 > 1e-8 || m_d3_s8 > 1e-8 || m_d3_s9 > 1e-8) { - if (m_damping.compare("bj")) { - param = dftd3_new_rational_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp); - } else if (m_damping.compare("zero")) { - param = dftd3_new_zero_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp); - } else if (m_damping.compare("bjm")) { - param = dftd3_new_mrational_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp); - } else if (m_damping.compare("zerom")) { - param = dftd3_new_mzero_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp, m_bet); - } else if (m_damping.compare("op")) { - param = dftd3_new_optimizedpower_damping(m_error, m_d3_s6, m_d3_s8, m_d3_s9, m_d3_a1, m_d3_a2, m_d3_alp, m_bet); - } - } else { - char* cstr = new char[m_functional.length() + 1]; - strcpy(cstr, m_functional.c_str()); - if (m_damping.compare("bj") == 0) { - param = dftd3_load_rational_damping(m_error, cstr, m_atm); - } else if (m_damping.compare("zero") == 0) { - param = dftd3_load_zero_damping(m_error, cstr, m_atm); - } else if (m_damping.compare("bjm") == 0) { - param = dftd3_load_mrational_damping(m_error, cstr, m_atm); - } else if (m_damping.compare("zerom") == 0) { - param = dftd3_load_mzero_damping(m_error, cstr, m_atm); - } else if (m_damping.compare("op") == 0) { - param = dftd3_load_optimizedpower_damping(m_error, cstr, m_atm); - } - delete[] cstr; - } dftd3_update_structure(m_error, m_mol, m_coord, NULL); - dftd3_get_dispersion(m_error, m_mol, m_disp, param, &energy, grad, sigma); - - dftd3_delete_param(¶m); + dftd3_get_dispersion(m_error, m_mol, m_disp, m_param, &energy, grad, sigma); return energy; } diff --git a/src/core/dftd3interface.h b/src/core/dftd3interface.h index c029f7d..01258de 100644 --- a/src/core/dftd3interface.h +++ b/src/core/dftd3interface.h @@ -63,6 +63,8 @@ class DFTD3Interface { inline double ParameterS9() const { return m_d3_s9; } private: + void CreateParameter(); + std::string m_damping; std::string m_functional; @@ -83,4 +85,5 @@ class DFTD3Interface { dftd3_error m_error; dftd3_structure m_mol; dftd3_model m_disp; + dftd3_param m_param; }; diff --git a/src/core/dftd4interface.cpp b/src/core/dftd4interface.cpp index ce87855..b2e89cb 100644 --- a/src/core/dftd4interface.cpp +++ b/src/core/dftd4interface.cpp @@ -101,6 +101,7 @@ double DFTD4Interface::DFTD4Calculation(double* grad) double energy = 0; dftd4::TCutoff cutoff; dftd4::TD4Model d4; + // exit(0); dftd4::get_dispersion(m_mol, m_charge, d4, m_par, cutoff, energy, grad); return energy; } diff --git a/src/core/eigen_uff.cpp b/src/core/eigen_uff.cpp index 11bd66e..920075c 100644 --- a/src/core/eigen_uff.cpp +++ b/src/core/eigen_uff.cpp @@ -560,7 +560,6 @@ double UFFThread::CalculateElectrostatic() eigenUFF::eigenUFF(const json& controller) { json parameter = MergeJson(UFFParameterJson, controller); - std::cout << parameter["threads"] << std::endl; m_threadpool = new CxxThreadPool(); m_threadpool->setProgressBar(CxxThreadPool::ProgressBarType::None); #ifdef USE_D3 diff --git a/src/core/energycalculator.cpp b/src/core/energycalculator.cpp index 1f4b280..bc7ff27 100644 --- a/src/core/energycalculator.cpp +++ b/src/core/energycalculator.cpp @@ -109,7 +109,8 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro #endif } else if (std::find(m_qmdff_method.begin(), m_qmdff_method.end(), m_method) != m_qmdff_method.end()) { // Just D4 energy calculator requested m_qmdff = new QMDFF(controller); - m_qmdff->setParameter(m_parameter); + if (m_parameter.size()) + m_qmdff->setParameter(m_parameter); m_ecengine = [this](bool gradient, bool verbose) { this->CalculateQMDFF(gradient, verbose); }; @@ -118,6 +119,7 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro m_uff = new eigenUFF(controller); } } + EnergyCalculator::~EnergyCalculator() { if (std::find(m_uff_methods.begin(), m_uff_methods.end(), m_method) != m_uff_methods.end()) { // UFF energy calculator requested @@ -139,6 +141,8 @@ EnergyCalculator::~EnergyCalculator() #ifdef USE_D4 delete m_d4; #endif + } else if (std::find(m_qmdff_method.begin(), m_qmdff_method.end(), m_method) != m_qmdff_method.end()) { // Just D4 energy calculator requested + delete m_qmdff; } else { // Fall back to UFF? delete m_uff; } @@ -197,7 +201,7 @@ void EnergyCalculator::setMolecule(const Molecule& molecule) #endif } else if (std::find(m_qmdff_method.begin(), m_qmdff_method.end(), m_method) != m_qmdff_method.end()) { // m_qmdff->setMolecule(atoms, geom); - // m_qmdff->Initialise(); + m_qmdff->Initialise(); } else { // Fall back to UFF? } diff --git a/src/core/qmdff.cpp b/src/core/qmdff.cpp index d1a3589..5ae8f85 100644 --- a/src/core/qmdff.cpp +++ b/src/core/qmdff.cpp @@ -28,9 +28,8 @@ */ #include "hbonds.h" -#ifdef USE_D4 +#include "src/core/dftd3interface.h" #include "src/core/dftd4interface.h" -#endif #include "src/core/qmdff_par.h" #include "src/core/topology.h" @@ -60,7 +59,7 @@ int QMDFFThread::execute() return 0; } -void QMDFFThread::readUFF(const json& parameter) +void QMDFFThread::readFF(const json& parameter) { // json parameter = MergeJson(UFFParameterJson, parameters); /* @@ -90,7 +89,6 @@ double QMDFFThread::StretchEnergy(double rAB, double reAB, double kAB, double a) { double ratio = reAB / rAB; double energy = kAB * (1 + pow(ratio, a) - 2 * pow(ratio, a * 0.5)); - // double energy = (0.5 * kij * (distance - r) * (distance - r)) * m_final_factor * m_bond_scaling; if (isnan(energy)) return 0; else @@ -118,12 +116,10 @@ double QMDFFThread::CalculateStretchEnergy() Vector x = Position(a); Vector y = Position(b); Matrix derivate; - double r0 = // BondStretching(x, y, derivate, m_CalculateGradient); - - energy += StretchEnergy(Distance(xa, xb, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA); //(0.5 * bond.kij * (r0 - bond.r0) * (r0 - bond.r0)) * m_final_factor * m_bond_scaling; + energy += StretchEnergy(Distance(xa, xb, ya, yb, za, zb), bond.reAB, bond.kAB, bond.exponA); if (m_CalculateGradient) { if (m_calc_gradient == 0) { - double diff = 0; //(bond.kij) * (r0 - bond.r0) * m_final_factor * m_bond_scaling; + double diff = 0; m_gradient.row(a) += diff * derivate.row(0); m_gradient.row(b) += diff * derivate.row(1); @@ -213,11 +209,19 @@ double QMDFFThread::CalculateAngleBending() 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 { - return this->LinearAngleBend(a, b, c, thetae, kabc, reAB, reAC); + 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 { - return this->AngleBend(a, b, c, thetae, kabc, reAB, reAC); + 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; @@ -629,43 +633,54 @@ double QMDFFThread::CalculateElectrostatic() QMDFF::QMDFF(const json& controller) { - json parameter = MergeJson(UFFParameterJson, controller); + json parameter = MergeJson(QMDFFParameterJson, controller); + // std::cout << parameter << std::endl; + m_threadpool = new CxxThreadPool(); m_threadpool->setProgressBar(CxxThreadPool::ProgressBarType::None); std::string param_file = parameter["param_file"]; - std::string uff_file = parameter["uff_file"]; + std::string uff_file = parameter["ff_file"]; if (param_file.compare("none") != 0) { readParameterFile(param_file); } if (uff_file.compare("none") != 0) { - readUFFFile(uff_file); + readFFFile(uff_file); } m_final_factor = 1 / 2625.15 * 4.19; m_d = parameter["differential"].get(); + json d3settings = DFTD3Settings; + // the following parameter are simply taken from the HF calculation with def2-SVP basis set (later should be not important) + d3settings["d_s6"] = 1.0; + d3settings["d_s8"] = s8; + d3settings["d_s9"] = 1.0; - readUFF(parameter); + d3settings["d_a1"] = a1; + d3settings["d_a2"] = a2; + d3settings["d_alp"] = 1; + + m_d3 = new DFTD3Interface(d3settings); + m_h4correction.allocate(m_atom_types.size()); + + readFF(parameter); m_writeparam = parameter["writeparam"]; - m_writeuff = parameter["writeuff"]; + m_writeuff = parameter["writeff"]; m_verbose = parameter["verbose"]; m_rings = parameter["rings"]; m_threads = parameter["threads"]; - m_scaling = 1.4; - // m_au = au; -} - -QMDFF::QMDFF() -{ - m_threadpool = new CxxThreadPool(); + m_const = parameter["const"]; + m_variable = parameter["variable"]; + // std::cout << "Calculating variable " << int(m_variable) << " and constants terms " << int(m_const) << std::endl; m_scaling = 1.4; // m_au = au; } QMDFF::~QMDFF() { + delete m_d3; delete m_threadpool; for (int i = 0; i < m_stored_threads.size(); ++i) delete m_stored_threads[i]; @@ -686,7 +701,7 @@ void QMDFF::setParameter(const json& parameter) json QMDFF::writeParameter() const { - json parameters = writeUFF(); + json parameters = writeFF(); parameters["bonds"] = Bonds(); // parameters["angles"] = Angles(); // parameters["dihedrals"] = Dihedrals(); @@ -698,7 +713,7 @@ void QMDFF::Initialise() { if (m_initialised) return; - std::cout << "Initialising QMDFF (see S. Grimmme, J. Chem. Theory Comput. 2014, 10, 10, 4497–4514 [10.1021/ct500573f]) for the original publication!" << std::endl; + // std::cout << "Initialising QMDFF (see S. Grimmme, J. Chem. Theory Comput. 2014, 10, 10, 4497–4514 [10.1021/ct500573f]) for the original publication!" << std::endl; m_uff_atom_types = std::vector(m_atom_types.size(), 0); m_coordination = std::vector(m_atom_types.size(), 0); @@ -738,293 +753,57 @@ void QMDFF::Initialise() if (m_rings) m_identified_rings = Topology::FindRings(m_stored_bonds, m_atom_types.size()); + /* + bonds.clean(); + setBonds(bonds, ignored_vdw, angles, dihedrals, inversions); - bonds.clean(); - setBonds(bonds, ignored_vdw, angles, dihedrals, inversions); - - angles.clean(); - setAngles(angles, ignored_vdw); - - dihedrals.clean(); - setDihedrals(dihedrals); + angles.clean(); + setAngles(angles, ignored_vdw); - inversions.clean(); - setInversions(inversions); + dihedrals.clean(); + setDihedrals(dihedrals); - nonbonds.clean(); - setvdWs(ignored_vdw); + inversions.clean(); + setInversions(inversions); + nonbonds.clean(); + setvdWs(ignored_vdw); + */ if (m_writeparam.compare("none") != 0) writeParameterFile(m_writeparam + ".json"); if (m_writeuff.compare("none") != 0) - writeUFFFile(m_writeuff + ".json"); + writeFFFile(m_writeuff + ".json"); AutoRanges(); m_initialised = true; } -void QMDFF::setBonds(const TContainer& bonds, std::vector>& ignored_vdw, TContainer& angels, TContainer& dihedrals, TContainer& inversions) -{ - for (const auto& bond : bonds.Storage()) { - QMDFFBond b; - - b.a = bond[0]; - b.b = bond[1]; - double xa = (m_geometry)(b.a, 0) * m_au; - double xb = (m_geometry)(b.b, 0) * m_au; - - double ya = (m_geometry)(b.a, 1) * m_au; - double yb = (m_geometry)(b.b, 1) * m_au; - - double za = (m_geometry)(b.a, 2) * m_au; - double zb = (m_geometry)(b.b, 2) * m_au; - - b.reAB = Topology::Distance(xa, xb, ya, yb, za, zb); // BondRestLength(b.i, b.j, bond_order); - b.kAB = 100; // 0.5 * m_bond_force * cZi * cZj / (b.r0 * b.r0 * b.r0); - double kaA = ka(m_atom_types[b.a]); - double kaB = ka(m_atom_types[b.b]); - double dEN = Elements::PaulingEN[m_atom_types[b.a]] - Elements::PaulingEN[m_atom_types[b.b]]; - b.exponA = kaA * kaB + kEN * dEN * dEN; - b.distance = 0; - m_qmdffbonds.push_back(b); - - int i = bond[0]; - int j = bond[1]; - - std::vector k_bodies; - for (auto t : m_stored_bonds[i]) { - k_bodies.push_back(t); - - if (t == j) - continue; - angels.insert({ i, std::min(t, j), std::max(j, t) }); - ignored_vdw[i].insert(t); - } - - std::vector l_bodies; - for (auto t : m_stored_bonds[j]) { - l_bodies.push_back(t); - - if (t == i) - continue; - angels.insert({ std::min(i, t), j, std::max(t, i) }); - ignored_vdw[j].insert(t); - } - - for (int k : k_bodies) { - for (int l : l_bodies) { - if (k == i || k == j || k == l || i == j || i == l || j == l) - continue; - dihedrals.insert({ k, i, j, l }); - - ignored_vdw[i].insert(k); - ignored_vdw[i].insert(l); - ignored_vdw[j].insert(k); - ignored_vdw[j].insert(l); - ignored_vdw[k].insert(l); - ignored_vdw[l].insert(k); - } - } - if (m_stored_bonds[i].size() == 3) { - inversions.insert({ i, m_stored_bonds[i][0], m_stored_bonds[i][1], m_stored_bonds[i][2] }); - } - if (m_stored_bonds[j].size() == 3) { - inversions.insert({ j, m_stored_bonds[j][0], m_stored_bonds[j][1], m_stored_bonds[j][2] }); - } - } -} - -void QMDFF::setAngles(const TContainer& angles, const std::vector>& ignored_vdw) -{ - for (const auto& angle : angles.Storage()) { - - QMDFFAngle a; - - a.a = angle[0]; - a.b = angle[1]; - a.c = angle[2]; - if (a.a == a.b || a.a == a.c || a.b == a.c) - continue; - - QMDFFBond b; - - b.a = angle[1]; - b.b = angle[2]; - b.distance = 1; -#pragma message("which distance should be used") - double xa = (m_geometry)(a.a, 0) * m_au; - double xb = (m_geometry)(a.b, 0) * m_au; - double xc = (m_geometry)(a.c, 0) * m_au; - - double ya = (m_geometry)(a.a, 1) * m_au; - double yb = (m_geometry)(a.b, 1) * m_au; - double yc = (m_geometry)(a.c, 1) * m_au; - - double za = (m_geometry)(a.a, 2) * m_au; - double zb = (m_geometry)(a.b, 2) * m_au; - double zc = (m_geometry)(a.c, 2) * m_au; - - b.reAB = Topology::Distance(xb, xc, yb, yc, zb, zc); // BondRestLength(b.i, b.j, bond_order); - b.kAB = 100; // 0.5 * m_bond_force * cZi * cZj / (b.r0 * b.r0 * b.r0); - double kaA = ka(m_atom_types[b.a]); - double kaB = ka(m_atom_types[b.b]); - double dEN = Elements::PaulingEN[m_atom_types[b.a]] - Elements::PaulingEN[m_atom_types[b.b]]; - b.exponA = ka13 + kb13 * kaA * kaB; - b.distance = 1; - m_qmdffbonds.push_back(b); - - a.reAB = Topology::Distance(xa, xb, ya, yb, za, zb); - a.reAC = Topology::Distance(xa, xc, ya, yc, za, zc); - auto atom_0 = m_geometry.row(a.a); - auto atom_1 = m_geometry.row(a.b); - auto atom_2 = m_geometry.row(a.c); - - 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] }; - - a.thetae = (vec_1.dot(vec_2) / sqrt(vec_1.dot(vec_1)) * sqrt(vec_2.dot(vec_2))) * 360 / 2.0 / pi; - a.kabc = 1; - - m_qmdffangle.push_back(a); - } -} - -void QMDFF::setDihedrals(const TContainer& dihedrals) -{ - for (const auto& dihedral : dihedrals.Storage()) { - UFFDihedral d; - d.i = dihedral[0]; - d.j = dihedral[1]; - d.k = dihedral[2]; - d.l = dihedral[3]; - - d.n = 2; - double f = pi / 180.0; - double bond_order = 1; - d.V = 2; - d.n = 3; - d.phi0 = 180 * f; - - if (std::find(Conjugated.cbegin(), Conjugated.cend(), m_uff_atom_types[d.k]) != Conjugated.cend() && std::find(Conjugated.cbegin(), Conjugated.cend(), m_uff_atom_types[d.j]) != Conjugated.cend()) - bond_order = 2; - else if (std::find(Triples.cbegin(), Triples.cend(), m_uff_atom_types[d.k]) != Triples.cend() || std::find(Triples.cbegin(), Triples.cend(), m_uff_atom_types[d.j]) != Triples.cend()) - bond_order = 3; - else - bond_order = 1; - - if (m_coordination[d.j] == 4 && m_coordination[d.k] == 4) // 2*sp3 - { - d.V = sqrt(UFFParameters[m_uff_atom_types[d.j]][cV] * UFFParameters[m_uff_atom_types[d.k]][cV]); - d.phi0 = 180 * f; - d.n = 3; - } - if (m_coordination[d.j] == 3 && m_coordination[d.k] == 3) // 2*sp2 - { - d.V = 5 * sqrt(UFFParameters[m_uff_atom_types[d.j]][cU] * UFFParameters[m_uff_atom_types[d.k]][cU]) * (1 + 4.18 * log(bond_order)); - d.phi0 = 180 * f; - d.n = 2; - } else if ((m_coordination[d.j] == 4 && m_coordination[d.k] == 3) || (m_coordination[d.j] == 3 && m_coordination[d.k] == 4)) { - d.V = sqrt(UFFParameters[m_uff_atom_types[d.j]][cV] * UFFParameters[m_uff_atom_types[d.k]][cV]); - d.phi0 = 0 * f; - d.n = 6; - - } else { - d.V = 5 * sqrt(UFFParameters[m_uff_atom_types[d.j]][cU] * UFFParameters[m_uff_atom_types[d.k]][cU]) * (1 + 4.18 * log(bond_order)); - d.phi0 = 90 * f; - } - - m_uffdihedral.push_back(d); - } -} - -void QMDFF::setInversions(const TContainer& inversions) +void QMDFF::setMolecule(const std::vector& atom_types, const std::vector>& geometry) { - for (const auto& inversion : inversions.Storage()) { - const int i = inversion[0]; - if (m_coordination[i] != 3) - continue; - - UFFInversion inv; - inv.i = i; - inv.j = inversion[1]; - inv.k = inversion[2]; - inv.l = inversion[3]; - - double C0 = 0.0; - double C1 = 0.0; - double C2 = 0.0; - double f = pi / 180.0; - double kijkl = 0; - if (6 <= m_atom_types[i] && m_atom_types[i] <= 8) { - C0 = 1.0; - C1 = -1.0; - C2 = 0.0; - kijkl = 6; - if (m_atom_types[inv.j] == 8 || m_atom_types[inv.k] == 8 || m_atom_types[inv.l] == 8) - kijkl = 50; - } else { - double w0 = pi / 180.0; - switch (m_atom_types[i]) { - // if the central atom is phosphorous - case 15: - w0 *= 84.4339; - break; - - // if the central atom is arsenic - case 33: - w0 *= 86.9735; - break; - - // if the central atom is antimonium - case 51: - w0 *= 87.7047; - break; - - // if the central atom is bismuth - case 83: - w0 *= 90.0; - break; - } - C2 = 1.0; - C1 = -4.0 * cos(w0 * f); - C0 = -(C1 * cos(w0 * f) + C2 * cos(2.0 * w0 * f)); - kijkl = 22.0 / (C0 + C1 + C2); - } - inv.C0 = C0; - inv.C1 = C1; - inv.C2 = C2; - inv.kijkl = kijkl; - m_uffinversion.push_back(inv); + m_atom_types = atom_types; + m_geometry = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); + for (int i = 0; i < m_atom_types.size(); ++i) { + m_geometry(i, 0) = geometry[i][0]; + m_geometry(i, 1) = geometry[i][1]; + m_geometry(i, 2) = geometry[i][2]; } + m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); + m_d3->InitialiseMolecule(m_atom_types); + m_h4correction.allocate(m_atom_types.size()); } -void QMDFF::setvdWs(const std::vector>& ignored_vdw) +void QMDFF::setMolecule(const std::vector& atom_types, const Matrix& geometry) { - for (int i = 0; i < m_atom_types.size(); ++i) { - for (int j = i + 1; j < m_atom_types.size(); ++j) { - // if (std::find(ignored_vdw[i].begin(), ignored_vdw[i].end(), j) != ignored_vdw[i].end() || std::find(ignored_vdw[j].begin(), ignored_vdw[j].end(), i) != ignored_vdw[j].end()) - // continue; - UFFvdW v; - v.i = i; - v.j = j; - - double cDi = UFFParameters[m_uff_atom_types[v.i]][cD]; - double cDj = UFFParameters[m_uff_atom_types[v.j]][cD]; - double cxi = UFFParameters[m_uff_atom_types[v.i]][cx]; - double cxj = UFFParameters[m_uff_atom_types[v.j]][cx]; - v.Dij = sqrt(cDi * cDj) * 2; - - v.xij = sqrt(cxi * cxj); - - m_uffvdwaals.push_back(v); - } - } + m_atom_types = atom_types; + m_geometry = geometry; + m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); + m_d3->InitialiseMolecule(m_atom_types); } void QMDFF::writeParameterFile(const std::string& file) const { - json parameters = writeUFF(); + json parameters = writeFF(); parameters["bonds"] = Bonds(); parameters["angles"] = Angles(); parameters["dihedrals"] = Dihedrals(); @@ -1034,11 +813,11 @@ void QMDFF::writeParameterFile(const std::string& file) const parameterfile << parameters; } -void QMDFF::writeUFFFile(const std::string& file) const +void QMDFF::writeFFFile(const std::string& file) const { - json parameters = writeUFF(); + json parameters = writeFF(); std::ofstream parameterfile(file); - parameterfile << writeUFF(); + parameterfile << writeFF(); } json QMDFF::Bonds() const @@ -1128,21 +907,22 @@ json QMDFF::vdWs() const return vdws; } -json QMDFF::writeUFF() const +json QMDFF::writeFF() const { json parameters; return parameters; } -void QMDFF::readUFF(const json& parameters) +void QMDFF::readFF(const json& parameters) { json parameter = MergeJson(UFFParameterJson, parameters); + setParameter(parameter); } void QMDFF::readParameter(const json& parameters) { m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); - readUFF(parameters); + readFF(parameters); AutoRanges(); m_initialised = true; @@ -1220,7 +1000,7 @@ void QMDFF::setInversions(const json& inversions) m_uffinversion.push_back(inv); } } - +/* void QMDFF::setvdWs(const json& vdws) { m_uffvdwaals.clear(); @@ -1236,8 +1016,8 @@ void QMDFF::setvdWs(const json& vdws) m_uffvdwaals.push_back(v); } } - -void QMDFF::readUFFFile(const std::string& file) +*/ +void QMDFF::readFFFile(const std::string& file) { nlohmann::json parameters; std::ifstream parameterfile(file); @@ -1246,7 +1026,7 @@ void QMDFF::readUFFFile(const std::string& file) } catch (nlohmann::json::type_error& e) { } catch (nlohmann::json::parse_error& e) { } - readUFF(parameters); + readFF(parameters); } void QMDFF::readParameterFile(const std::string& file) @@ -1265,31 +1045,35 @@ void QMDFF::AutoRanges() { for (int i = 0; i < m_threads; ++i) { QMDFFThread* thread = new QMDFFThread(i, m_threads); - thread->readUFF(writeUFF()); + thread->readFF(writeFF()); thread->setMolecule(m_atom_types, &m_geometry); m_threadpool->addThread(thread); m_stored_threads.push_back(thread); - for (int j = int(i * m_qmdffbonds.size() / double(m_threads)); j < int((i + 1) * m_qmdffbonds.size() / double(m_threads)); ++j) - thread->AddBond(m_qmdffbonds[j]); - - for (int j = int(i * m_qmdffangle.size() / double(m_threads)); j < int((i + 1) * m_qmdffangle.size() / double(m_threads)); ++j) - thread->AddAngle(m_qmdffangle[j]); - - for (int j = int(i * m_uffdihedral.size() / double(m_threads)); j < int((i + 1) * m_uffdihedral.size() / double(m_threads)); ++j) - thread->AddDihedral(m_uffdihedral[j]); - - for (int j = int(i * m_uffinversion.size() / double(m_threads)); j < int((i + 1) * m_uffinversion.size() / double(m_threads)); ++j) - thread->AddInversion(m_uffinversion[j]); - - for (int j = int(i * m_uffvdwaals.size() / double(m_threads)); j < int((i + 1) * m_uffvdwaals.size() / double(m_threads)); ++j) - thread->AddvdW(m_uffvdwaals[j]); + if (m_variable) + for (int j = int(i * m_qmdffbonds.size() / double(m_threads)); j < int((i + 1) * m_qmdffbonds.size() / double(m_threads)); ++j) + thread->AddBond(m_qmdffbonds[j]); + if (m_variable) + for (int j = int(i * m_qmdffangle.size() / double(m_threads)); j < int((i + 1) * m_qmdffangle.size() / double(m_threads)); ++j) + thread->AddAngle(m_qmdffangle[j]); + + if (m_const) + for (int j = int(i * m_uffdihedral.size() / double(m_threads)); j < int((i + 1) * m_uffdihedral.size() / double(m_threads)); ++j) + thread->AddDihedral(m_uffdihedral[j]); + + if (m_const) + for (int j = int(i * m_uffinversion.size() / double(m_threads)); j < int((i + 1) * m_uffinversion.size() / double(m_threads)); ++j) + thread->AddInversion(m_uffinversion[j]); + /* + for (int j = int(i * m_uffvdwaals.size() / double(m_threads)); j < int((i + 1) * m_uffvdwaals.size() / double(m_threads)); ++j) + thread->AddvdW(m_uffvdwaals[j]); + */ } m_uff_bond_end = m_qmdffbonds.size(); m_uff_angle_end = m_qmdffangle.size(); m_uff_dihedral_end = m_uffdihedral.size(); m_uff_inv_end = m_uffinversion.size(); - m_uff_vdw_end = m_uffvdwaals.size(); + // m_uff_vdw_end = m_uffvdwaals.size(); } void QMDFF::UpdateGeometry(const double* coord) @@ -1378,6 +1162,21 @@ Eigen::MatrixXd QMDFF::NumGrad() double QMDFF::Calculate(bool grd, bool verbose) { m_CalculateGradient = grd; + hbonds4::atom_t geometry[m_atom_types.size()]; + for (int i = 0; i < m_atom_types.size(); ++i) { + geometry[i].x = m_geometry(i, 0) * m_au; + geometry[i].y = m_geometry(i, 1) * m_au; + geometry[i].z = m_geometry(i, 2) * m_au; + geometry[i].e = m_atom_types[i]; + m_h4correction.GradientH4()[i].x = 0; + m_h4correction.GradientH4()[i].y = 0; + m_h4correction.GradientH4()[i].z = 0; + + m_h4correction.GradientHH()[i].x = 0; + m_h4correction.GradientHH()[i].y = 0; + m_h4correction.GradientHH()[i].z = 0; + m_d3->UpdateAtom(i, m_geometry(i, 0), m_geometry(i, 1), m_geometry(i, 2)); + } double energy = 0.0; @@ -1385,7 +1184,9 @@ double QMDFF::Calculate(bool grd, bool verbose) double angle_energy = 0.0; double dihedral_energy = 0.0; double inversion_energy = 0.0; - double vdw_energy = 0.0; + // double vdw_energy = 0.0; + double d3_energy = 0; + m_threadpool->setActiveThreadCount(m_threads); for (int i = 0; i < m_stored_threads.size(); ++i) { @@ -1400,19 +1201,52 @@ double QMDFF::Calculate(bool grd, bool verbose) angle_energy += m_stored_threads[i]->AngleEnergy(); dihedral_energy += m_stored_threads[i]->DihedralEnergy(); inversion_energy += m_stored_threads[i]->InversionEnergy(); - vdw_energy += m_stored_threads[i]->VdWEnergy(); + // vdw_energy += m_stored_threads[i]->VdWEnergy(); m_gradient += m_stored_threads[i]->Gradient(); } /* + CalculateElectrostatic(); */ - energy = bond_energy + angle_energy + dihedral_energy + inversion_energy + vdw_energy; + double energy_h4 = 0; + double energy_hh = 0; + + if (m_const) { + + if (m_h4_scaling > 1e-8) + energy_h4 = m_h4correction.energy_corr_h4(m_atom_types.size(), geometry); + if (m_hh_scaling > 1e-8) + m_h4correction.energy_corr_hh_rep(m_atom_types.size(), geometry); + for (int i = 0; i < m_atom_types.size(); ++i) { + m_gradient(i, 0) += m_final_factor * m_h4_scaling * m_h4correction.GradientH4()[i].x + m_final_factor * m_hh_scaling * m_h4correction.GradientHH()[i].x; + m_gradient(i, 1) += m_final_factor * m_h4_scaling * m_h4correction.GradientH4()[i].y + m_final_factor * m_hh_scaling * m_h4correction.GradientHH()[i].y; + m_gradient(i, 2) += m_final_factor * m_h4_scaling * m_h4correction.GradientH4()[i].z + m_final_factor * m_hh_scaling * m_h4correction.GradientHH()[i].z; + } + + if (grd) { + double grad[3 * m_atom_types.size()]; + d3_energy = m_d3->DFTD3Calculation(grad); + for (int i = 0; i < m_atom_types.size(); ++i) { + double val = grad[3 * i + 0] * au; + if (!std::isnan(val) && std::abs(val) < 1e10) + m_gradient(i, 0) += val; + val = grad[3 * i + 1] * au; + if (!std::isnan(val) && std::abs(val) < 1e10) + m_gradient(i, 1) += val; + val = grad[3 * i + 2] * au; + if (!std::isnan(val) && std::abs(val) < 1e10) + m_gradient(i, 2) += val; + } + } else + d3_energy = m_d3->DFTD3Calculation(0); + energy += m_final_factor * m_h4_scaling * energy_h4 + m_final_factor * m_hh_scaling * energy_hh + d3_energy; + } + energy = bond_energy + angle_energy + dihedral_energy + inversion_energy; // + vdw_energy; if (verbose) { std::cout << "Total energy " << energy << " Eh. Sum of " << std::endl << "Bond Energy " << bond_energy << " Eh" << std::endl << "Angle Energy " << angle_energy << " Eh" << std::endl << "Dihedral Energy " << dihedral_energy << " Eh" << std::endl << "Inversion Energy " << inversion_energy << " Eh" << std::endl - << "Nonbonded Energy " << vdw_energy << " Eh" << std::endl + << "Nonbonded Energy " << d3_energy + m_final_factor * m_h4_scaling * energy_h4 + m_final_factor * m_hh_scaling * energy_hh << " Eh" << std::endl << std::endl; for (int i = 0; i < m_atom_types.size(); ++i) { diff --git a/src/core/qmdff.h b/src/core/qmdff.h index 828ef8d..2a87edb 100644 --- a/src/core/qmdff.h +++ b/src/core/qmdff.h @@ -33,6 +33,9 @@ #include "external/CxxThreadPool/include/CxxThreadPool.h" +#include "src/core/dftd3interface.h" +#include "src/core/dftd4interface.h" + #include "src/core/uff_par.h" #include #include @@ -62,7 +65,7 @@ class QMDFFThread : public CxxThread { void UpdateGeometry(Matrix* geometry) { m_geometry = geometry; - m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); + m_gradient = Eigen::MatrixXd::Zero(geometry->rows(), 3); } Matrix Gradient() const { return m_gradient; } @@ -73,7 +76,8 @@ class QMDFFThread : public CxxThread { m_geometry = geometry; m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); } - void readUFF(const json& parameters); + + void readFF(const json& parameters); inline double Energy() const { return m_energy; } inline double BondEnergy() const { return m_bond_energy; } @@ -211,31 +215,14 @@ class QMDFFThread : public CxxThread { class QMDFF { public: QMDFF(const json& controller); - QMDFF(); ~QMDFF(); void UpdateGeometry(const double* coord); void UpdateGeometry(const std::vector>& geometry); - void setMolecule(const std::vector& atom_types, const std::vector>& geometry) - { - m_atom_types = atom_types; - m_geometry = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); - for (int i = 0; i < m_atom_types.size(); ++i) { - m_geometry(i, 0) = geometry[i][0]; - m_geometry(i, 1) = geometry[i][1]; - m_geometry(i, 2) = geometry[i][2]; - } - m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); - } - - void setMolecule(const std::vector& atom_types, const Matrix& geometry) - { - m_atom_types = atom_types; - m_geometry = geometry; - m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); - } + void setMolecule(const std::vector& atom_types, const std::vector>& geometry); + void setMolecule(const std::vector& atom_types, const Matrix& geometry); void setParameter(const json& parameter); json writeParameter() const; @@ -253,10 +240,10 @@ class QMDFF { Eigen::MatrixXd NumGrad(); void writeParameterFile(const std::string& file) const; - void writeUFFFile(const std::string& file) const; + void writeFFFile(const std::string& file) const; // json writeParameter() const; - json writeUFF() const; + json writeFF() const; json Bonds() const; json Angles() const; @@ -265,25 +252,16 @@ class QMDFF { json vdWs() const; void setBonds(const json& bonds); - void setBonds(const TContainer& bonds, std::vector>& ignored_vdw, TContainer& angels, TContainer& dihedrals, TContainer& inversions); - void setAngles(const json& angles); - void setAngles(const TContainer& angles, const std::vector>& ignored_vdw); - void setDihedrals(const json& dihedrals); - void setDihedrals(const TContainer& dihedrals); - void setInversions(const json& inversions); - void setInversions(const TContainer& inversions); - - void setvdWs(const json& vdws); - void setvdWs(const std::vector>& ignored_vdw); + // void setvdWs(const json& vdws); void readParameterFile(const std::string& file); - void readUFFFile(const std::string& file); + void readFFFile(const std::string& file); void readParameter(const json& parameters); - void readUFF(const json& parameters); + void readFF(const json& parameters); void setInitialisation(bool initialised) { m_initialised = initialised; } @@ -363,6 +341,9 @@ class QMDFF { int m_calc_gradient = 0; int m_threads = 1; + bool m_variable = true, m_const = true; std::vector m_stored_threads; CxxThreadPool* m_threadpool; + hbonds4::H4Correction m_h4correction; + DFTD3Interface* m_d3; }; diff --git a/src/core/qmdff_par.h b/src/core/qmdff_par.h index 166f4ec..dbe8794 100644 --- a/src/core/qmdff_par.h +++ b/src/core/qmdff_par.h @@ -107,3 +107,49 @@ inline double kZ(int element) else // row 5 return 0.6; } + +const json QMDFFParameterJson{ + { "bond_scaling", 1 }, + { "angle_scaling", 1 }, + { "dihedral_scaling", 1 }, + { "inversion_scaling", 1 }, + { "vdw_scaling", 1 }, + { "rep_scaling", 1 }, + { "coulomb_scaling", 1 }, + { "bond_force", 664.12 }, + { "angle_force", 664.12 }, + { "differential", 1e-7 }, + { "h4_scaling", 0 }, + { "hh_scaling", 0 }, + { "h4_oh_o", 2.32 }, + { "h4_oh_n", 3.10 }, + { "h4_nh_o", 1.07 }, + { "h4_nh_n", 2.01 }, + { "h4_wh_o", 0.42 }, + { "h4_nh4", 3.61 }, + { "h4_coo", 1.41 }, + { "hh_rep_k", 0.42 }, + { "hh_rep_e", 12.7 }, + { "hh_rep_r0", 2.3 }, + { "d4", 0 }, + { "d3", 0 }, + { "d_s6", 1.00 }, + { "d_s8", 1.20065498 }, + { "d_s10", 0 }, + { "d_s9", 1 }, + { "d_a1", 0.40085597 }, + { "d_a2", 5.02928789 }, + { "d_alp", 16 }, + { "d_func", "pbe0" }, + { "d_atm", true }, + { "param_file", "none" }, + { "ff_file", "none" }, + { "writeparam", "none" }, + { "writeff", "none" }, + { "verbose", false }, + { "rings", false }, + { "threads", 1 }, + { "gradient", 0 }, + { "variable", true }, + { "const", true } +};