use partial charges to calculate dipole moment + atom-wise scaling

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

src/capabilities/curcumaopt.cpp

@@ -124,7 +124,6 @@ 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();
@@ -133,18 +132,6 @@ void CurcumaOpt::ProcessMoleculesSerial(const std::vector<Molecule>& molecules)
                       << "Dipole momement (GFN2)" << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) << std::endl;
         }
 #endif
-        Molecule mol(*iter);
-        mol.setPartialCharges(interface.Charges());
-        auto dipoles = mol.CalculateDipoleMoments();
-        for (const auto& dipole : dipoles) {
-            std::cout << std::endl
-                      << std::endl
-                      << "Dipole momement for single molecule " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
-        }
-        auto dipole = mol.CalculateDipoleMoment();
-        std::cout << std::endl
-                  << std::endl
-                  << "Dipole momement for whole structure " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
 
         if (m_hessian) {
             Hessian hess(m_method, m_defaults, m_threads);
@@ -243,29 +230,16 @@ double CurcumaOpt::SinglePoint(const Molecule* initial, const json& controller,
     interface.setMolecule(*initial);
 
     double energy = interface.CalculateEnergy(true, true);
-    std::cout << "dipole1" << std::endl;
-
+    double store = 0;
 #ifdef USE_TBLITE
     if (method.compare("gfn2") == 0) {
         std::vector<double> dipole = interface.Dipole();
         std::cout << std::endl
                   << std::endl
                   << "Dipole momement (GNF2)" << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
+        store = sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]);
     }
 #endif
-    Molecule mol(*initial);
-    mol.setPartialCharges(interface.Charges());
-    auto dipoles = mol.CalculateDipoleMoments();
-    for (const auto& dipole : dipoles) {
-        std::cout << std::endl
-                  << std::endl
-                  << "Dipole momement for single molecule " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
-    }
-    auto dipole = mol.CalculateDipoleMoment();
-    std::cout << std::endl
-              << std::endl
-              << "Dipole momement for whole structure " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
-
     return energy;
 }
 

src/capabilities/optimiser/OptimiseDipoleScaling.h

@@ -0,0 +1,128 @@
+/*
+ * <LevenbergMarquardt Optimsation for Dipole Calculation from Partial Charges. >
+ * Copyright (C) 2023 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
+ * 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <unsupported/Eigen/NonLinearOptimization>
+
+#include <iostream>
+
+#include "src/core/elements.h"
+#include "src/core/global.h"
+#include "src/core/molecule.h"
+#include "src/core/pseudoff.h"
+
+#include "src/tools/geometry.h"
+
+#include "json.hpp"
+#include <LBFGS.h>
+#include <LBFGSB.h>
+
+using json = nlohmann::json;
+
+using Eigen::VectorXd;
+using namespace LBFGSpp;
+
+class LBFGSDipoleInterface {
+public:
+    LBFGSDipoleInterface(int n_)
+    {
+    }
+    double operator()(const VectorXd& x, VectorXd& grad)
+    {
+        double fx = 0.0;
+        double dx = 5e-1;
+        std::vector<double> scaling(x.size());
+        for (int i = 0; i < scaling.size(); ++i) {
+            scaling[i] = x(i);
+        }
+
+        auto dipole_vector = m_molecule->CalculateDipoleMoment(scaling);
+        double dipole = sqrt(dipole_vector[0] * dipole_vector[0] + dipole_vector[1] * dipole_vector[1] + dipole_vector[2] * dipole_vector[2]) * 2.5418;
+        fx = std::abs(m_dipole - dipole);
+        // std::cout << dipole << " " << m_dipole << " "<< fx<< std::endl;
+        for (int i = 0; i < scaling.size(); ++i) {
+            // if(std::abs(fx) < std::abs(m_smaller))
+            //     std::cout << scaling[i] << " ";
+            scaling[i] += dx;
+            dipole_vector = m_molecule->CalculateDipoleMoment(scaling);
+            double dipolep = sqrt(dipole_vector[0] * dipole_vector[0] + dipole_vector[1] * dipole_vector[1] + dipole_vector[2] * dipole_vector[2]) * 2.5418;
+            scaling[i] -= 2 * dx;
+            dipole_vector = m_molecule->CalculateDipoleMoment(scaling);
+            double dipolem = sqrt(dipole_vector[0] * dipole_vector[0] + dipole_vector[1] * dipole_vector[1] + dipole_vector[2] * dipole_vector[2]) * 2.5418;
+
+            grad[i] = (std::abs(dipolep - dipole) - std::abs(dipolem - dipole)) / (2.0 * dx);
+            scaling[i] += dx;
+        }
+        // if(std::abs(fx) < std::abs(m_smaller))
+        //     std::cout << std::endl;
+        m_smaller = std::abs(fx);
+        m_parameter = x;
+        return fx;
+    }
+
+    Vector Parameter() const { return m_parameter; }
+    void setMolecule(const Molecule* molecule)
+    {
+        m_molecule = molecule;
+    }
+
+    double m_dipole;
+    double m_smaller = 1;
+
+private:
+    int m_atoms = 0;
+    Vector m_parameter;
+    const Molecule* m_molecule;
+};
+
+inline std::pair<double, std::vector<double>> OptimiseScaling(const Molecule* molecule, double dipole, double initial, double threshold, int maxiter)
+{
+    Vector parameter(molecule->AtomCount());
+    for (int i = 0; i < molecule->AtomCount(); ++i)
+        parameter(i) = initial;
+
+    Vector old_param = parameter;
+    // std::cout << parameter.transpose() << std::endl;
+
+    LBFGSParam<double> param;
+    param.epsilon = 1e-5;
+    LBFGSSolver<double> solver(param);
+    LBFGSDipoleInterface fun(parameter.size());
+    fun.m_dipole = dipole * 2.5418;
+    fun.setMolecule(molecule);
+    int iteration = 0;
+    // std::cout << parameter.transpose() << std::endl;
+    double fx;
+    int converged = solver.InitializeSingleSteps(fun, parameter, fx);
+
+    for (iteration = 1; iteration <= maxiter && fx > threshold; ++iteration) {
+        solver.SingleStep(fun, parameter, fx);
+        std::cout << "Iteration: " << iteration << " Difference: " << fx << std::endl;
+    }
+    std::vector<double> scaling(parameter.size());
+    for (int i = 0; i < scaling.size(); ++i)
+        scaling[i] = parameter(i);
+    auto dipole_vector = molecule->CalculateDipoleMoment(scaling);
+    dipole = sqrt(dipole_vector[0] * dipole_vector[0] + dipole_vector[1] * dipole_vector[1] + dipole_vector[2] * dipole_vector[2]);
+
+    return std::pair<double, std::vector<double>>(dipole, scaling);
+}

src/core/molecule.cpp

@@ -823,10 +823,9 @@ Eigen::Vector3d Molecule::COM(bool protons, int fragment)
     return com;
 }
 
-std::vector<Position> Molecule::CalculateDipoleMoments()
+std::vector<Position> Molecule::CalculateDipoleMoments(const std::vector<double>& scaling) const
 {
     std::vector<Position> dipole_moments;
-    double quickndirty_scaling = 3;
 
     if (m_charges.size() != m_geometry.size()) {
         std::cout << "No partial charges available" << std::endl;
@@ -847,18 +846,22 @@ std::vector<Position> Molecule::CalculateDipoleMoments()
         pos(1) /= mass;
         pos(2) /= mass;
         for (int i : m_fragments[f]) {
-            dipole(0) += m_charges[i] * (m_geometry[i][0] - pos(0)) * quickndirty_scaling;
-            dipole(1) += m_charges[i] * (m_geometry[i][1] - pos(1)) * quickndirty_scaling;
-            dipole(2) += m_charges[i] * (m_geometry[i][2] - pos(2)) * quickndirty_scaling;
+            double scale = 3;
+            if (scaling.size() > i)
+                scale = scaling[i];
+            dipole(0) += m_charges[i] * (m_geometry[i][0] - pos(0)) * scale;
+            dipole(1) += m_charges[i] * (m_geometry[i][1] - pos(1)) * scale;
+            dipole(2) += m_charges[i] * (m_geometry[i][2] - pos(2)) * scale;
+            //  std::cout << scale << " ";
         }
+        // std::cout << std::endl;
         dipole_moments.push_back(dipole);
     }
     return dipole_moments;
 }
 
-Position Molecule::CalculateDipoleMoment()
+Position Molecule::CalculateDipoleMoment(const std::vector<double>& scaling) const
 {
-    double quickndirty_scaling = 3;
     double mass = 0;
 
     Position pos = { 0, 0, 0 }, dipole = { 0, 0, 0 };
@@ -877,10 +880,15 @@ Position Molecule::CalculateDipoleMoment()
     pos(1) /= mass;
     pos(2) /= mass;
     for (int i = 0; i < m_geometry.size(); ++i) {
-        dipole(0) += m_charges[i] * (m_geometry[i][0] - pos(0)) * quickndirty_scaling;
-        dipole(1) += m_charges[i] * (m_geometry[i][1] - pos(1)) * quickndirty_scaling;
-        dipole(2) += m_charges[i] * (m_geometry[i][2] - pos(2)) * quickndirty_scaling;
+        double scale = 3;
+        if (scaling.size() > i)
+            scale = scaling[i];
+        dipole(0) += m_charges[i] * (m_geometry[i][0] - pos(0)) * scale;
+        dipole(1) += m_charges[i] * (m_geometry[i][1] - pos(1)) * scale;
+        dipole(2) += m_charges[i] * (m_geometry[i][2] - pos(2)) * scale;
+        // std::cout << scale << " ";
     }
+    // std::cout << std::endl;
     return dipole;
 }
 

src/core/molecule.h

@@ -128,8 +128,8 @@ class Molecule
     std::string XYZString() const;
     std::string XYZString(const std::vector<int> &order) const;
 
-    std::vector<Position> CalculateDipoleMoments();
-    Position CalculateDipoleMoment();
+    std::vector<Position> CalculateDipoleMoments(const std::vector<double>& scaling = std::vector<double>()) const;
+    Position CalculateDipoleMoment(const std::vector<double>& scaling = std::vector<double>()) const;
 
     std::vector<int> BoundHydrogens(int atom, double scaling = 1.5) const;
     std::map<int, std::vector<int>> getConnectivtiy(double scaling = 1.5, int latest = -1) const;

src/main.cpp

@@ -37,6 +37,8 @@
 #include "src/tools/general.h"
 #include "src/tools/info.h"
 
+#include "src/capabilities/optimiser/OptimiseDipoleScaling.h"
+
 #include <iostream>
 #include <string>
 #include <cstring>
@@ -833,6 +835,56 @@ int main(int argc, char **argv) {
                 std::cout << ":: " << gyr << " " << sum / double(count) << std::endl;
                 count++;
             }
+        } else if (strcmp(argv[1], "-dipole") == 0) {
+            FileIterator file(argv[2]);
+
+            double target = 0, threshold = 1e-1, initial = 3;
+            bool target_set = false;
+            int maxiter = 10;
+            json blob = controller["dipole"];
+            if (blob.contains("target")) {
+                target = blob["target"];
+                target_set = true;
+            }
+            if (blob.contains("threshold")) {
+                threshold = blob["threshold"];
+            }
+            if (blob.contains("initial")) {
+                initial = blob["initial"];
+            }
+            if (blob.contains("maxiter")) {
+                maxiter = blob["maxiter"];
+            }
+
+            while (!file.AtEnd()) {
+                Molecule mol = file.Next();
+                EnergyCalculator interface("gfn2", blob);
+                interface.setMolecule(mol);
+                interface.CalculateEnergy(false, true);
+                mol.setPartialCharges(interface.Charges());
+                std::vector<double> dipole_moment = interface.Dipole();
+                if (!target_set)
+                    target = sqrt(dipole_moment[0] * dipole_moment[0] + dipole_moment[1] * dipole_moment[1] + dipole_moment[2] * dipole_moment[2]);
+                std::cout << "Target dipole moment: " << target << std::endl;
+                auto dipoles = mol.CalculateDipoleMoments();
+                for (const auto& dipole : dipoles) {
+                    std::cout << std::endl
+                              << std::endl
+                              << "Dipole momement for single molecule " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
+                }
+                auto dipole = mol.CalculateDipoleMoment();
+                std::cout << std::endl
+                          << std::endl
+                          << "Dipole momement for whole structure " << dipole[0] << " " << dipole[1] << " " << dipole[2] << " : " << sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]) * 2.5418 << std::endl;
+
+                auto result = OptimiseScaling(&mol, target, initial, threshold, maxiter);
+                std::cout << "Final dipole moment " << result.first * 2.5418 << std::endl;
+                std::cout << "Final dipole moment " << result.first << std::endl;
+
+                for (auto i : result.second)
+                    std::cout << i << " ";
+                std::cout << std::endl;
+            }
         } else {
             bool centered = false;
             for (std::size_t i = 2; i < argc; ++i) {
@@ -854,31 +906,6 @@ int main(int argc, char **argv) {
                           << std::endl;
                 std::cout << mol.COM().transpose() << std::endl;
                 std::cout << mol.GyrationRadius() << std::endl;
-                /*
-                Hessian hess("gfn2", UFFParameterJson);
-                hess.setMolecule(mol);
-                hess.CalculateHessian();
-                */
-                /*
-                UFF forcefield(UFFParameterJson);
-                forcefield.setMolecule(mol.Atoms(), mol.Coords());
-                // forcefield.readParameterFile("parameter.json");
-
-                forcefield.Initialise();
-                std::cout << forcefield.Calculate(true) << std::endl;
-                forcefield.writeParameterFile("parameter.json");
-                forcefield.readParameterFile("parameter.json");
-                std::cout << forcefield.Calculate(true) << std::endl;
-                */
-                /*        double grad[3 * mol.AtomCount()];
-                          forcefield.Gradient(grad);
-                          forcefield.NumGrad(grad);
-                          */
-                // CompactTopo(mol.HydrogenBondMatrix(-1,-1));
-                // std::cout << CompareTopoMatrix(mol.HydrogenBondMatrix(-1,-1),mol.HydrogenBondMatrix(-1,-1) ) << std::endl;
-                // auto m = mol.DistanceMatrix();
-                // std::cout << m.first << std::endl;
-                // std::cout << m.second << std::endl;
             }
         }
     }