first qmdff protoyp, without repulsion and torsion

Signed-off-by: Conrad Hübler <Conrad.Huebler@gmx.net>

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)

src/capabilities/hessian.cpp

@@ -25,16 +25,16 @@
 
 #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)
     , m_xj(xj)
     , 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);

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<void(void)> 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;

src/capabilities/optimiser/LevMarQMDFFFit.h

@@ -1,6 +1,6 @@
 /*
  * <LevenbergMarquardt QMDFF FC Fit. >
- * Copyright (C) 2023 Conrad Hübler <Conrad.Huebler@gmx.net>
+ * Copyright (C) 2023 - 2024 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
@@ -34,6 +34,9 @@
 
 #include "src/tools/geometry.h"
 
+#include "json.hpp"
+using json = nlohmann::json;
+
 template <typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic>
 
 struct FCFunctor {
@@ -68,11 +71,7 @@ struct MyFCFunctor : FCFunctor<double> {
     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<double> {
         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<MyFCFunctor> {
 };
 
-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<MyFCFunctor> numDiff(functor);
     Eigen::LevenbergMarquardt<Eigen::NumericalDiff<MyFCFunctor>> 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;
 }