diff --git a/src/core/dftd3interface.cpp b/src/core/dftd3interface.cpp index bbe3b6e..a5d4d9f 100644 --- a/src/core/dftd3interface.cpp +++ b/src/core/dftd3interface.cpp @@ -47,16 +47,21 @@ DFTD3Interface::DFTD3Interface(const json& controller) m_damping = parameter["d_damping"]; m_functional = parameter["d_func"]; CreateParameter(); +#ifdef USE_D3 m_error = dftd3_new_error(); +#endif } DFTD3Interface::DFTD3Interface() { +#ifdef USE_D3 m_error = dftd3_new_error(); +#endif } DFTD3Interface::~DFTD3Interface() { +#ifdef USE_D3 dftd3_delete_model(&m_disp); dftd3_delete_structure(&m_mol); dftd3_delete_error(&m_error); @@ -64,6 +69,7 @@ DFTD3Interface::~DFTD3Interface() delete m_coord; delete m_attyp; +#endif } void DFTD3Interface::PrintParameter() const @@ -73,6 +79,7 @@ void DFTD3Interface::PrintParameter() const void DFTD3Interface::CreateParameter() { +#ifdef USE_D3 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); @@ -101,6 +108,7 @@ void DFTD3Interface::CreateParameter() } delete[] cstr; } +#endif } void DFTD3Interface::UpdateParameters(const json& controller) @@ -148,9 +156,10 @@ bool DFTD3Interface::InitialiseMolecule(const std::vector& atomtypes) m_coord[3 * i + 2] = (3 * i + 2) / au; m_attyp[i] = atomtypes[i]; } - +#ifdef USE_D3 m_mol = dftd3_new_structure(m_error, atomtypes.size(), m_attyp, m_coord, NULL, NULL); m_disp = dftd3_new_d3_model(m_error, m_mol); +#endif return true; } @@ -167,15 +176,19 @@ double DFTD3Interface::DFTD3Calculation(double* grad) double energy = 0; double sigma[9]; +#ifdef USE_D3 dftd3_update_structure(m_error, m_mol, m_coord, NULL); dftd3_get_dispersion(m_error, m_mol, m_disp, m_param, &energy, grad, sigma); +#endif return energy; } void DFTD3Interface::clear() { +#ifdef USE_D3 dftd3_delete_model(&m_disp); dftd3_delete_structure(&m_mol); dftd3_delete_error(&m_error); +#endif } diff --git a/src/core/dftd3interface.h b/src/core/dftd3interface.h index 379f325..c60f62e 100644 --- a/src/core/dftd3interface.h +++ b/src/core/dftd3interface.h @@ -20,8 +20,9 @@ #pragma once #include "src/tools/general.h" - +#ifdef USE_D3 #include "s-dftd3.h" +#endif #include "src/core/molecule.h" @@ -85,8 +86,11 @@ class DFTD3Interface { double m_bet = 8; bool m_atm = false; + +#ifdef USE_D3 dftd3_error m_error; dftd3_structure m_mol; dftd3_model m_disp; dftd3_param m_param; +#endif }; diff --git a/src/core/forcefieldthread.cpp b/src/core/forcefieldthread.cpp index 8720591..9f4ec39 100644 --- a/src/core/forcefieldthread.cpp +++ b/src/core/forcefieldthread.cpp @@ -423,16 +423,21 @@ D3Thread::D3Thread(int thread, int threads) setAutoDelete(false); m_final_factor = 1 / 2625.15 * 4.19; m_d = 1e-7; +#ifdef USE_D3 m_d3 = new DFTD3Interface(); +#endif } D3Thread::~D3Thread() { +#ifdef USE_D3 delete m_d3; +#endif } int D3Thread::execute() { +#ifdef USE_D3 for (int i = 0; i < m_atom_types.size(); ++i) { m_d3->UpdateAtom(i, m_geometry(i, 0), m_geometry(i, 1), m_geometry(i, 2)); } @@ -447,7 +452,10 @@ int D3Thread::execute() } } else m_vdw_energy = m_d3->DFTD3Calculation(0); - +#else + std::cerr << "D3 is not included, sorry for that" << std::endl; + exit(1); +#endif return 0; } diff --git a/src/core/forcefieldthread.h b/src/core/forcefieldthread.h index a0847be..9ca94f7 100644 --- a/src/core/forcefieldthread.h +++ b/src/core/forcefieldthread.h @@ -168,18 +168,24 @@ class D3Thread : public ForceFieldThread { ~D3Thread(); void setParamater(const json& parameter) { +#ifdef USE_D3 m_d3->UpdateParametersD3(parameter); +#endif } void Initialise(const std::vector& atom_types) { m_atom_types = atom_types; +#ifdef USE_D3 m_d3->InitialiseMolecule(m_atom_types); +#endif } virtual int execute() override; private: +#ifdef USE_D3 DFTD3Interface* m_d3; +#endif std::vector m_atom_types; }; diff --git a/src/core/qmdff.cpp b/src/core/qmdff.cpp index 13bf070..3d0aa8d 100644 --- a/src/core/qmdff.cpp +++ b/src/core/qmdff.cpp @@ -665,8 +665,10 @@ QMDFF::QMDFF(const json& controller) d3settings["d_a1"] = a1; d3settings["d_a2"] = a2; d3settings["d_alp"] = 1; - +#ifdef USE_D3 m_d3 = new DFTD3Interface(d3settings); +#endif + m_h4correction.allocate(m_atom_types.size()); readFF(parameter); @@ -685,7 +687,10 @@ QMDFF::QMDFF(const json& controller) QMDFF::~QMDFF() { +#ifdef USE_D3 delete m_d3; +#endif + delete m_threadpool; for (int i = 0; i < m_stored_threads.size(); ++i) delete m_stored_threads[i]; @@ -805,7 +810,9 @@ void QMDFF::setMolecule(const std::vector& atom_types, const Matrix& geomet m_atom_types = atom_types; m_geometry = geometry; m_gradient = Eigen::MatrixXd::Zero(m_atom_types.size(), 3); +#ifdef USE_D3 m_d3->InitialiseMolecule(m_atom_types); +#endif } void QMDFF::writeParameterFile(const std::string& file) const @@ -1182,7 +1189,9 @@ double QMDFF::Calculate(bool grd, bool verbose) m_h4correction.GradientHH()[i].x = 0; m_h4correction.GradientHH()[i].y = 0; m_h4correction.GradientHH()[i].z = 0; +#ifdef USE_D3 m_d3->UpdateAtom(i, m_geometry(i, 0), m_geometry(i, 1), m_geometry(i, 2)); +#endif } double energy = 0.0; @@ -1230,7 +1239,10 @@ double QMDFF::Calculate(bool grd, bool verbose) if (grd) { double grad[3 * m_atom_types.size()]; +#ifdef USE_D3 d3_energy = m_d3->DFTD3Calculation(grad); +#endif + 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) @@ -1243,7 +1255,10 @@ double QMDFF::Calculate(bool grd, bool verbose) m_gradient(i, 2) += val; } } else +#ifdef USE_D3 d3_energy = m_d3->DFTD3Calculation(0); +#endif + 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; diff --git a/src/core/qmdff.h b/src/core/qmdff.h index cbecf04..528f1a8 100644 --- a/src/core/qmdff.h +++ b/src/core/qmdff.h @@ -344,5 +344,8 @@ class QMDFF { std::vector m_stored_threads; CxxThreadPool* m_threadpool; hbonds4::H4Correction m_h4correction; + +#ifdef USE_D3 DFTD3Interface* m_d3; +#endif };