update tblite/xtb and extract dipole from xtb

Signed-off-by: Conrad Huebler <[email protected]>

CMakeLists.txt

@@ -206,6 +206,7 @@ if(USE_XTB)
       )
       FetchContent_MakeAvailable("mctc-lib")
     endif()
+    include_directories(SYSTEM ${PROJECT_BINARY_DIR}/_deps/tblite-src/include/)
 
     if(HELPERS)
 	add_executable(xtb_helper

src/capabilities/curcumaopt.cpp

@@ -110,6 +110,7 @@ void CurcumaOpt::start()
 void CurcumaOpt::ProcessMoleculesSerial(const std::vector<Molecule>& molecules)
 {
     EnergyCalculator interface(Json2KeyWord<std::string>(m_defaults, "method"), m_controller["sp"]);
+    std::string method = Json2KeyWord<std::string>(m_defaults, "method");
 
     auto iter = molecules.begin();
     interface.setMolecule(*iter);
@@ -122,6 +123,15 @@ void CurcumaOpt::ProcessMoleculesSerial(const std::vector<Molecule>& molecules)
         auto start = std::chrono::system_clock::now();
         interface.updateGeometry(iter->Coords());
         double energy = interface.CalculateEnergy(true, true);
+        std::cout << "dipole?" << std::endl;
+#ifdef USE_TBLITE
+        if (method.compare("gfn2") == 0) {
+            std::vector<double> dipole = interface.Dipole();
+            std::cout << std::endl
+                      << std::endl
+                      << "Dipole momement " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) << std::endl;
+        }
+#endif
         if (m_hessian) {
             Hessian hess(m_method, m_defaults, m_threads);
             hess.setMolecule(*iter);
@@ -215,7 +225,18 @@ double CurcumaOpt::SinglePoint(const Molecule* initial, const json& controller,
     EnergyCalculator interface(method, controller);
     interface.setMolecule(*initial);
 
-    return interface.CalculateEnergy(true, true);
+    double energy = interface.CalculateEnergy(true, true);
+    std::cout << "dipole1" << std::endl;
+
+#ifdef USE_TBLITE
+    if (method.compare("gfn2") == 0) {
+        std::vector<double> dipole = interface.Dipole();
+        std::cout << std::endl
+                  << std::endl
+                  << "Dipole momement " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
+    }
+#endif
+    return energy;
 }
 
 Molecule CurcumaOpt::LBFGSOptimise(const Molecule* initial, const json& controller, std::string& output, std::vector<Molecule>* intermediate)

src/capabilities/simplemd.cpp

@@ -85,6 +85,7 @@ void SimpleMD::LoadControlJson()
 
     m_writerestart = Json2KeyWord<int>(m_defaults, "writerestart");
     m_respa = Json2KeyWord<int>(m_defaults, "respa");
+    m_dipole = Json2KeyWord<bool>(m_defaults, "dipole");
 
     m_writeXYZ = Json2KeyWord<bool>(m_defaults, "writeXYZ");
     m_writeinit = Json2KeyWord<bool>(m_defaults, "writeinit");
@@ -530,7 +531,16 @@ void SimpleMD::start()
     }
     if (m_thermostat.compare("csvr") == 0)
         std::cout << "Exchange with heat bath " << m_Ekin_exchange << "Eh" << std::endl;
-
+    if (m_dipole) {
+        /*
+        double dipole = 0.0;
+        for( auto d : m_collected_dipole)
+            dipole += d;
+        dipole /= m_collected_dipole.size();
+        std::cout << dipole*2.5418 << " average dipole in Debye and " << dipole*2.5418*3.3356e-30 << " Cm" << std::endl;
+        */
+        std::cout << "Calculated averaged dipole moment " << m_aver_dipol * 2.5418 << " Debye and " << m_aver_dipol * 2.5418 * 3.3356 << " Cm [e-30]" << std::endl;
+    }
     std::ofstream restart_file("curcuma_final.json");
     restart_file << WriteRestartInformation() << std::endl;
 
@@ -779,7 +789,10 @@ void SimpleMD::PrintStatus() const
 #endif
     } else {
 #ifdef GCC
-        std::cout << fmt::format("{1: ^{0}f} {2: ^{0}f} {3: ^{0}f} {4: ^{0}f} {5: ^{0}f} {6: ^{0}f} {7: ^{0}f} {8: ^{0}f} {9: ^{0}f} {10: ^{0}f}\n", 15, m_currentStep / 1000, m_Epot, m_aver_Epot, m_Ekin, m_aver_Ekin, m_Etot, m_aver_Etot, m_T, m_aver_Temp, remaining);
+        if (m_dipole)
+            std::cout << fmt::format("{1: ^{0}f} {2: ^{0}f} {3: ^{0}f} {4: ^{0}f} {5: ^{0}f} {6: ^{0}f} {7: ^{0}f} {8: ^{0}f} {9: ^{0}f} {10: ^{0}f}  {11: ^{0}f}\n", 15, m_currentStep / 1000, m_Epot, m_aver_Epot, m_Ekin, m_aver_Ekin, m_Etot, m_aver_Etot, m_T, m_aver_Temp, m_aver_dipol * 2.5418 * 3.3356, remaining);
+        else
+            std::cout << fmt::format("{1: ^{0}f} {2: ^{0}f} {3: ^{0}f} {4: ^{0}f} {5: ^{0}f} {6: ^{0}f} {7: ^{0}f} {8: ^{0}f} {9: ^{0}f} {10: ^{0}f} \n", 15, m_currentStep / 1000, m_Epot, m_aver_Epot, m_Ekin, m_aver_Ekin, m_Etot, m_aver_Etot, m_T, m_aver_Temp, remaining);
 #else
         std::cout << m_currentStep * m_timestep / fs2amu / 1000 << " " << m_Epot << " " << m_Ekin << " " << m_Epot + m_Ekin << m_T << std::endl;
 
@@ -802,6 +815,11 @@ double SimpleMD::Gradient(const double* coord, double* grad)
 
     double Energy = m_interface->CalculateEnergy(true);
     m_interface->getGradient(grad);
+    if (m_dipole) {
+        std::vector<double> dipole = m_interface->Dipole();
+        m_curr_dipole = sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]);
+        m_collected_dipole.push_back(sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]));
+    }
     return Energy;
 }
 
@@ -819,7 +837,9 @@ double SimpleMD::EKin()
     m_aver_Epot = (m_Epot + (m_currentStep)*m_aver_Epot) / (m_currentStep + 1);
     m_aver_Ekin = (m_Ekin + (m_currentStep)*m_aver_Ekin) / (m_currentStep + 1);
     m_aver_Etot = (m_Etot + (m_currentStep)*m_aver_Etot) / (m_currentStep + 1);
-
+    if (m_dipole) {
+        m_aver_dipol = (m_curr_dipole + (m_currentStep)*m_aver_dipol) / (m_currentStep + 1);
+    }
     return ekin;
 }
 

src/capabilities/simplemd.h

@@ -62,7 +62,8 @@ static json CurcumaMDJson{
     { "rattle", false },
     { "rattle_tolerance", 1e-5 },
     { "thermostat", "csvr" },
-    { "respa", 1 }
+    { "respa", 1 },
+    { "dipole", false }
 };
 
 class SimpleMD : public CurcumaMethod {
@@ -135,7 +136,7 @@ class SimpleMD : public CurcumaMethod {
     std::string m_basename;
     int m_natoms = 0;
     int m_dumb = 1;
-    double m_T = 0, m_Epot = 0, m_aver_Epot = 0, m_Ekin = 0, m_aver_Ekin = 0, m_Etot = 0, m_aver_Etot = 0;
+    double m_T = 0, m_Epot = 0, m_aver_Epot = 0, m_Ekin = 0, m_aver_Ekin = 0, m_Etot = 0, m_aver_Etot = 0, m_aver_dipol = 0, m_curr_dipole = 0;
     int m_hmass = 4;
     double m_single_step = 1;
     double m_timestep = 0.5, m_currentStep = 0, m_maxtime = 1000;
@@ -157,11 +158,12 @@ class SimpleMD : public CurcumaMethod {
     double m_pos_conv = 0, m_scale_velo = 1.0, m_coupling = 10;
     double m_impuls = 0, m_impuls_scaling = 0.75, m_dt2 = 0;
     double m_rattle_tolerance = 1e-4;
-
+    std::vector<double> m_collected_dipole;
     Matrix m_topo_initial;
     std::vector<Molecule*> m_unique_structures;
     std::string m_method = "UFF", m_initfile = "none", m_thermostat = "csvr";
     bool m_unstable = false;
+    bool m_dipole = false;
 };
 
 class MDThread : public CxxThread {

src/core/energycalculator.cpp

@@ -32,6 +32,16 @@
 EnergyCalculator::EnergyCalculator(const std::string& method, const json& controller)
     : m_method(method)
 {
+    m_charges = []() {
+        return std::vector<double>{};
+    };
+    m_dipole = []() {
+        return std::vector<double>{};
+    };
+    m_bonds = []() {
+        return std::vector<std::vector<double>>{ {} };
+    };
+
     if (std::find(m_uff_methods.begin(), m_uff_methods.end(), m_method) != m_uff_methods.end()) { // UFF energy calculator requested
         m_uff = new eigenUFF(controller);
         m_ecengine = [this](bool gradient, bool verbose) {
@@ -43,6 +53,15 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro
         m_ecengine = [this](bool gradient, bool verbose) {
             this->CalculateTBlite(gradient, verbose);
         };
+        m_charges = [this]() {
+            return this->m_tblite->Charges();
+        };
+        m_dipole = [this]() {
+            return this->m_tblite->Dipole();
+        };
+        m_bonds = [this]() {
+            return this->m_tblite->BondOrders();
+        };
 #else
         std::cout << "TBlite was not included ..." << std::endl;
         exit(1);
@@ -54,6 +73,15 @@ EnergyCalculator::EnergyCalculator(const std::string& method, const json& contro
         m_ecengine = [this](bool gradient, bool verbose) {
             this->CalculateXTB(gradient, verbose);
         };
+        m_charges = [this]() {
+            return this->m_xtb->Charges();
+        };
+        m_dipole = [this]() {
+            return this->m_xtb->Dipole();
+        };
+        m_bonds = [this]() {
+            return this->m_xtb->BondOrders();
+        };
 #else
         std::cout << "XTB was not included ..." << std::endl;
         exit(1);
@@ -307,3 +335,18 @@ Matrix EnergyCalculator::Gradient() const
 {
     return m_eigen_gradient;
 }
+
+std::vector<double> EnergyCalculator::Charges() const
+{
+    return m_charges();
+}
+
+std::vector<double> EnergyCalculator::Dipole() const
+{
+    return m_dipole();
+}
+
+std::vector<std::vector<double>> EnergyCalculator::BondOrders() const
+{
+    return m_bonds();
+}

src/core/energycalculator.h

@@ -92,6 +92,11 @@ class EnergyCalculator {
     }
 #endif
 
+    std::vector<double> Charges() const;
+    std::vector<double> Dipole() const;
+
+    std::vector<std::vector<double>> BondOrders() const;
+
 private:
     void InitialiseUFF();
     void CalculateUFF(bool gradient, bool verbose = false);
@@ -131,6 +136,9 @@ class EnergyCalculator {
     StringList m_d3_methods = { "d3" };
     StringList m_d4_methods = { "d4" };
     std::function<void(bool, bool)> m_ecengine;
+    std::function<std::vector<double>()> m_charges, m_dipole;
+    std::function<std::vector<std::vector<double>>()> m_bonds;
+
     std::string m_method;
     std::vector<std::array<double, 3>> m_geometry, m_gradient;
     Matrix m_eigen_geometry, m_eigen_gradient;

src/core/tbliteinterface.cpp

@@ -17,7 +17,9 @@
  *
  */
 
-#include "external/tblite/include/tblite.h"
+#ifndef tblite_delete
+#include "tblite.h"
+#endif
 
 #include "src/core/global.h"
 #include "src/tools/general.h"
@@ -181,6 +183,17 @@ std::vector<double> TBLiteInterface::Charges() const
     return charges;
 }
 
+std::vector<double> TBLiteInterface::Dipole() const
+{
+    std::vector<double> dipole(3);
+    double* c = new double[3];
+    tblite_get_result_dipole(m_error, m_tblite_res, c);
+    for (int i = 0; i < 3; ++i)
+        dipole[i] = c[i];
+    delete[] c;
+    return dipole;
+}
+
 std::vector<std::vector<double>> TBLiteInterface::BondOrders() const
 {
     std::vector<std::vector<double>> bond_orders(m_atomcount);

src/core/tbliteinterface.h

@@ -21,7 +21,9 @@
 
 #include "src/tools/general.h"
 
-#include "external/tblite/include/tblite.h"
+#ifndef tblite_delete
+#include "tblite.h"
+#endif
 
 #include "src/core/molecule.h"
 
@@ -59,6 +61,8 @@ class TBLiteInterface {
     void clear();
 
     std::vector<double> Charges() const;
+    std::vector<double> Dipole() const;
+
     std::vector<std::vector<double>> BondOrders() const;
 
 private:

src/core/xtbinterface.cpp

@@ -142,6 +142,43 @@ double XTBInterface::GFNCalculation(int parameter, double* grad)
     return energy;
 }
 
+std::vector<double> XTBInterface::Charges() const
+{
+    std::vector<double> charges(m_atomcount);
+    double* c = new double[m_atomcount];
+    xtb_getCharges(m_env, m_xtb_res, c);
+    for (int i = 0; i < m_atomcount; ++i)
+        charges[i] = c[i];
+    delete[] c;
+    return charges;
+}
+
+std::vector<double> XTBInterface::Dipole() const
+{
+    std::vector<double> dipole(3);
+    double* c = new double[3];
+    xtb_getDipole(m_env, m_xtb_res, c);
+    for (int i = 0; i < 3; ++i)
+        dipole[i] = c[i];
+    delete[] c;
+    return dipole;
+}
+
+std::vector<std::vector<double>> XTBInterface::BondOrders() const
+{
+    std::vector<std::vector<double>> bond_orders(m_atomcount);
+    double* bonds = new double[m_atomcount * m_atomcount];
+    xtb_getBondOrders(m_env, m_xtb_res, bonds);
+    for (int i = 0; i < m_atomcount; ++i) {
+        std::vector<double> b(m_atomcount);
+        for (int j = 0; j < m_atomcount; ++j)
+            b[j] = bonds[i * j];
+        bond_orders[i] = b;
+    }
+    delete[] bonds;
+    return bond_orders;
+}
+
 void XTBInterface::clear()
 {
     xtb_delResults(&m_xtb_res);

src/core/xtbinterface.h

@@ -53,6 +53,10 @@ class XTBInterface {
     double GFNCalculation(int parameter = 2, double* grad = 0);
 
     void clear();
+    std::vector<double> Charges() const;
+    std::vector<double> Dipole() const;
+
+    std::vector<std::vector<double>> BondOrders() const;
 
 private:
     int m_atomcount = 0;