Surface forces (after rebasing #488)

doc/users/user_ref.rst

@@ -458,6 +458,33 @@ User input reference
 
     **Default** ``user['GeometryOptimizer']['run']``
 
+ :Forces: Define parameters for the computation of forces. 
+
+  :red:`Keywords`
+   :method: Method for computing forces. ``surface_integrals`` (more accurate) uses surface integrals over the quantum mechanical stress tensor, while ``hellmann_feynman`` uses the Hellmann-Feynman theorem. 
+
+    **Type** ``str``
+
+    **Default** ``surface_integrals``
+
+    **Predicates**
+      - ``value.lower() in ['surface_integrals', 'hellmann_feynman']``
+
+   :surface_integral_precision: Precision of the surface integrals used in the computation of forces. Determines the number of Lebedev grid points used in the surface integration. 
+
+    **Type** ``str``
+
+    **Default** ``medium``
+
+    **Predicates**
+      - ``value.lower() in ['low', 'medium', 'high']``
+
+   :radius_factor: Sets the radius of the surface used in the computation of forces. The radius is given by this factor times the distance to the neariest neighbour. Must be between 0.1 and 0.9. This should rarely need to be changed. Different values can change the accuracy of the forces. 
+
+    **Type** ``float``
+
+    **Default** ``0.6``
+
  :ExternalFields: Define external electromagnetic fields. 
 
   :red:`Keywords`

python/mrchem/helpers.py

@@ -287,7 +287,7 @@ def write_scf_solver(user_dict, wf_dict):
 def write_scf_properties(user_dict, origin):
     prop_dict = {}
     if user_dict["Properties"]["dipole_moment"]:
-        prop_dict["dipole_moment"] = {}
+        prop_dict["dipole_moment"] = {} 
         prop_dict["dipole_moment"]["dip-1"] = {
             "operator": "h_e_dip",
             "precision": user_dict["world_prec"],
@@ -306,6 +306,9 @@ def write_scf_properties(user_dict, origin):
             "operator": "h_nuc_grad",
             "precision": user_dict["world_prec"],
             "smoothing": user_dict["Precisions"]["nuclear_prec"],
+            "method": user_dict["Forces"]["method"],
+            "surface_integral_precision": user_dict["Forces"]["surface_integral_precision"],
+            "radius_factor": user_dict["Forces"]["radius_factor"],
         }
     if user_dict["Properties"]["hirshfeld_charges"]:
         prop_dict["hirshfeld_charges"] = {}

python/mrchem/input_parser/api.py

@@ -370,6 +370,24 @@ def stencil() -> JSONDict:
                                             'name': 'geometric_derivative',
                                             'type': 'bool'}],
                         'name': 'Properties'},
+                    {   'keywords': [   {   'default': 'surface_integrals',
+                                            'name': 'method',
+                                            'predicates': [   'value.lower() '
+                                                              'in '
+                                                              "['surface_integrals', "
+                                                              "'hellmann_feynman']"],
+                                            'type': 'str'},
+                                        {   'default': 'medium',
+                                            'name': 'surface_integral_precision',
+                                            'predicates': [   'value.lower() '
+                                                              "in ['low', "
+                                                              "'medium', "
+                                                              "'high']"],
+                                            'type': 'str'},
+                                        {   'default': 0.5,
+                                            'name': 'radius_factor',
+                                            'type': 'float'}],
+                        'name': 'Forces'},
                     {   'keywords': [   {   'default': [],
                                             'name': 'electric_field',
                                             'predicates': [   'len(value) == 0 '

python/mrchem/input_parser/docs/user_ref.rst

@@ -458,6 +458,33 @@ User input reference
 
     **Default** ``user['GeometryOptimizer']['run']``
 
+ :Forces: Define parameters for the computation of forces. 
+
+  :red:`Keywords`
+   :method: Method for computing forces. ``surface_integrals`` (more accurate) uses surface integrals over the quantum mechanical stress tensor, while ``hellmann_feynman`` uses the Hellmann-Feynman theorem. 
+
+    **Type** ``str``
+
+    **Default** ``surface_integrals``
+
+    **Predicates**
+      - ``value.lower() in ['surface_integrals', 'hellmann_feynman']``
+
+   :surface_integral_precision: Precision of the surface integrals used in the computation of forces. Determines the number of Lebedev grid points used in the surface integration. 
+
+    **Type** ``str``
+
+    **Default** ``medium``
+
+    **Predicates**
+      - ``value.lower() in ['low', 'medium', 'high']``
+
+   :radius_factor: Sets the radius of the surface used in the computation of forces. The radius is given by this factor times the distance to the neariest neighbour. Must be between 0.1 and 0.9. This should rarely need to be changed. Different values can change the accuracy of the forces. 
+
+    **Type** ``float``
+
+    **Default** ``0.6``
+
  :ExternalFields: Define external electromagnetic fields. 
 
   :red:`Keywords`

python/template.yml

@@ -416,6 +416,34 @@ sections:
         default: false
         docstring: |
           Compute Hirshfeld charges.
+  - name: Forces
+    docstring: |
+      Define parameters for the computation of forces.
+    keywords:
+      - name: method
+        type: str
+        default: surface_integrals
+        predicates:
+          - value.lower() in ['surface_integrals', 'hellmann_feynman']
+        docstring: |
+          Method for computing forces. ``surface_integrals`` (more accurate) uses surface
+          integrals over the quantum mechanical stress tensor, while ``hellmann_feynman`` uses the Hellmann-Feynman
+          theorem.
+      - name: surface_integral_precision
+        type: str
+        default: medium
+        predicates:
+          - value.lower() in ['low', 'medium', 'high']
+        docstring: |
+          Precision of the surface integrals used in the computation of forces. Determines the number of Lebedev grid points used in the
+          surface integration.
+      - name: radius_factor
+        type: float
+        default: 0.6
+        docstring: |
+          Sets the radius of the surface used in the computation of forces.
+          The radius is given by this factor times the distance to the neariest neighbour. Must be between 0.1 and 0.9.
+          This should rarely need to be changed. Different values can change the accuracy of the forces.
   - name: ExternalFields
     docstring: |
       Define external electromagnetic fields.

src/CMakeLists.txt

@@ -36,17 +36,17 @@ add_subdirectory(properties)
 add_subdirectory(qmfunctions)
 add_subdirectory(qmoperators)
 add_subdirectory(scf_solver)
+add_subdirectory(surface_forces)
 add_subdirectory(tensor)
 add_subdirectory(utils)
 
 target_compile_definitions(mrchem PRIVATE 
     HIRSHFELD_SOURCE_DIR="${HIRSHFELD_SOURCE_DIR}"
     HIRSHFELD_INSTALL_DIR="${HIRSHFELD_INSTALL_DIR}"
-)
-
-target_compile_definitions(mrchem PRIVATE 
     AZORA_POTENTIALS_SOURCE_DIR="${AZORA_POTENTIALS_SOURCE_DIR}"
     AZORA_POTENTIALS_INSTALL_DIR="${AZORA_POTENTIALS_INSTALL_DIR}"
+    LEBVEDEV_SOURCE_DIR="${LEBVEDEV_SOURCE_DIR}"
+    LEBVEDEV_INSTALL_DIR="${LEBVEDEV_INSTALL_DIR}"
 )
 
 target_link_libraries(mrchem

src/driver.cpp

@@ -57,6 +57,7 @@
 #include "qmoperators/one_electron/NuclearOperator.h"
 #include "qmoperators/one_electron/ZoraOperator.h"
 #include "qmoperators/one_electron/AZoraPotential.h"
+#include "qmoperators/one_electron/NablaOperator.h"
 
 #include "qmoperators/one_electron/H_BB_dia.h"
 #include "qmoperators/one_electron/H_BM_dia.h"
@@ -83,8 +84,10 @@
 #include "environment/LPBESolver.h"
 #include "environment/PBESolver.h"
 #include "environment/Permittivity.h"
-
+#include "surface_forces/SurfaceForce.h"
 #include "properties/hirshfeld/HirshfeldPartition.h"
+#include <fstream>
+#include <filesystem>
 
 #include "mrdft/Factory.h"
 
@@ -117,7 +120,7 @@ namespace scf {
 bool guess_orbitals(const json &input, Molecule &mol);
 bool guess_energy(const json &input, Molecule &mol, FockBuilder &F);
 void write_orbitals(const json &input, Molecule &mol);
-void calc_properties(const json &input, Molecule &mol);
+void calc_properties(const json &input, Molecule &mol, const json &json_fock);
 void plot_quantities(const json &input, Molecule &mol);
 } // namespace scf
 
@@ -318,7 +321,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
 
     if (json_out["success"]) {
         if (json_scf.contains("write_orbitals")) scf::write_orbitals(json_scf["write_orbitals"], mol);
-        if (json_scf.contains("properties")) scf::calc_properties(json_scf["properties"], mol);
+        if (json_scf.contains("properties")) scf::calc_properties(json_scf["properties"], mol, json_fock);
         if (json_scf.contains("plots")) scf::plot_quantities(json_scf["plots"], mol);
     }
 
@@ -471,7 +474,7 @@ void driver::scf::write_orbitals(const json &json_orbs, Molecule &mol) {
  * input section, and will compute all properties which are present in this input.
  * This includes the diamagnetic contributions to the magnetic response properties.
  */
-void driver::scf::calc_properties(const json &json_prop, Molecule &mol) {
+void driver::scf::calc_properties(const json &json_prop, Molecule &mol, const json &json_fock) {
     Timer t_tot, t_lap;
     auto plevel = Printer::getPrintLevel();
     if (plevel == 1) mrcpp::print::header(1, "Computing ground state properties");
@@ -519,12 +522,29 @@ void driver::scf::calc_properties(const json &json_prop, Molecule &mol) {
     if (json_prop.contains("geometric_derivative")) {
         t_lap.start();
         mrcpp::print::header(2, "Computing geometric derivative");
-        for (const auto &item : json_prop["geometric_derivative"].items()) {
-            const auto &id = item.key();
-            const double &prec = item.value()["precision"];
-            const double &smoothing = item.value()["smoothing"];
-            GeometricDerivative &G = mol.getGeometricDerivative(id);
-            auto &nuc = G.getNuclear();
+        GeometricDerivative &G = mol.getGeometricDerivative("geom-1");
+        auto &nuc = G.getNuclear();
+
+        // calculate nuclear gradient
+        for (auto k = 0; k < mol.getNNuclei(); ++k) {
+            const auto nuc_k = nuclei[k];
+            auto Z_k = nuc_k.getCharge();
+            auto R_k = nuc_k.getCoord();
+            nuc.row(k) = Eigen::RowVector3d::Zero();
+            for (auto l = 0; l < mol.getNNuclei(); ++l) {
+                if (l == k) continue;
+                const auto nuc_l = nuclei[l];
+                auto Z_l = nuc_l.getCharge();
+                auto R_l = nuc_l.getCoord();
+                std::array<double, 3> R_kl = {R_k[0] - R_l[0], R_k[1] - R_l[1], R_k[2] - R_l[2]};
+                auto R_kl_3_2 = std::pow(math_utils::calc_distance(R_k, R_l), 3.0);
+                nuc.row(k) -= Eigen::Map<Eigen::RowVector3d>(R_kl.data()) * (Z_k * Z_l / R_kl_3_2);
+            }
+        }
+        // calculate electronic gradient using the Hellmann-Feynman theorem
+        if (json_prop["geometric_derivative"]["geom-1"]["method"] == "hellmann_feynman") {
+            const double &prec = json_prop["geometric_derivative"]["geom-1"]["precision"];
+            const double &smoothing = json_prop["geometric_derivative"]["geom-1"]["smoothing"];
             auto &el = G.getElectronic();
 
             for (auto k = 0; k < mol.getNNuclei(); ++k) {
@@ -534,19 +554,24 @@ void driver::scf::calc_properties(const json &json_prop, Molecule &mol) {
                 double c = detail::nuclear_gradient_smoothing(smoothing, Z_k, mol.getNNuclei());
                 NuclearGradientOperator h(Z_k, R_k, prec, c);
                 h.setup(prec);
-                nuc.row(k) = Eigen::RowVector3d::Zero();
-                for (auto l = 0; l < mol.getNNuclei(); ++l) {
-                    if (l == k) continue;
-                    const auto nuc_l = nuclei[l];
-                    auto Z_l = nuc_l.getCharge();
-                    auto R_l = nuc_l.getCoord();
-                    std::array<double, 3> R_kl = {R_k[0] - R_l[0], R_k[1] - R_l[1], R_k[2] - R_l[2]};
-                    auto R_kl_3_2 = std::pow(math_utils::calc_distance(R_k, R_l), 3.0);
-                    nuc.row(k) -= Eigen::Map<Eigen::RowVector3d>(R_kl.data()) * (Z_k * Z_l / R_kl_3_2);
-                }
                 el.row(k) = h.trace(Phi).real();
                 h.clear();
             }
+        // calculate electronic gradient using the surface integrals method
+        } else if (json_prop["geometric_derivative"]["geom-1"]["method"] == "surface_integrals") {
+            double prec = json_prop["geometric_derivative"]["geom-1"]["precision"];
+            std::string leb_prec = json_prop["geometric_derivative"]["geom-1"]["surface_integral_precision"];
+            double radius_factor = json_prop["geometric_derivative"]["geom-1"]["radius_factor"];
+            Eigen::MatrixXd surfaceForces = surface_force::surface_forces(mol, Phi, prec, json_fock, leb_prec, radius_factor);
+            GeometricDerivative &G = mol.getGeometricDerivative("geom-1");
+            auto &nuc = G.getNuclear();
+            auto &el = G.getElectronic();
+            // set electronic gradient
+            for (int k = 0; k < mol.getNNuclei(); k++) {
+                el.row(k) = surfaceForces.row(k) - nuc.row(k);
+            }
+        } else {
+            MSG_ABORT("Invalid method for geometric derivative");
         }
         mrcpp::print::footer(2, t_lap, 2);
         if (plevel == 1) mrcpp::print::time(1, "Geometric derivative", t_lap);
@@ -783,7 +808,7 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) {
     F_0.setup(unpert_prec);
     if (plevel == 1) mrcpp::print::footer(1, t_unpert, 2);
 
-    if (json_rsp.contains("properties")) scf::calc_properties(json_rsp["properties"], mol);
+    if (json_rsp.contains("properties")) scf::calc_properties(json_rsp["properties"], mol, unpert_fock);
 
     ///////////////////////////////////////////////////////////
     //////////////   Preparing Perturbed System   /////////////

src/mrdft/Functional.h

@@ -62,6 +62,8 @@ class Functional {
     }
     double XCenergy = 0.0;
 
+    Eigen::MatrixXd evaluate(Eigen::MatrixXd &inp) const;
+    Eigen::MatrixXd evaluate_transposed(Eigen::MatrixXd &inp) const;
     friend class MRDFT;
 
 protected:
@@ -78,8 +80,6 @@ class Functional {
     virtual int getCtrInputLength() const = 0;
     virtual int getCtrOutputLength() const = 0;
 
-    Eigen::MatrixXd evaluate(Eigen::MatrixXd &inp) const;
-    Eigen::MatrixXd evaluate_transposed(Eigen::MatrixXd &inp) const;
     Eigen::MatrixXd contract(Eigen::MatrixXd &xc_data, Eigen::MatrixXd &d_data) const;
     Eigen::MatrixXd contract_transposed(Eigen::MatrixXd &xc_data, Eigen::MatrixXd &d_data) const;
 

src/properties/GeometricDerivative.h

@@ -56,6 +56,8 @@ class GeometricDerivative final {
         mrcpp::print::separator(0, '-');
         print_utils::matrix(0, "Electronic", getElectronic());
         print_utils::scalar(0, "Norm", getElectronic().norm(), "(au)");
+        mrcpp::print::separator(0, '-');
+        print_utils::scalar(0, "Est. var. of force error", variance(), "(au)", 4, true);
         mrcpp::print::separator(0, '=', 2);
     }
 
@@ -68,6 +70,18 @@ class GeometricDerivative final {
                 {"electronic_norm", getElectronic().norm()}};
     }
 
+    /**
+     * @brief Compute the variance of the total force using the formula
+     * sigma = sqrt(1/(3N) * sum_i sum_j (F_{ij})^2) presented in
+     * Gubler et al.: Journal of Computational Physics: X Volume 17, November 2023, 100131
+    */
+    double variance() const {
+        DoubleMatrix forces = getTensor();
+        Eigen::Vector3d total_force = forces.colwise().sum();
+        double force_variance = total_force.norm() * std::sqrt( 1.0 / (3.0 * forces.rows()) );
+        return force_variance;
+    }
+
 protected:
     DoubleMatrix nuclear;
     DoubleMatrix electronic;

src/qmoperators/one_electron/HessianOperator.h

@@ -0,0 +1,44 @@
+#pragma once
+
+#include "tensor/RankOneOperator.h"
+
+#include "qmoperators/QMDerivative.h"
+#include "qmoperators/one_electron/NablaOperator.h"
+
+namespace mrchem {
+
+/** @class HessianOperator
+ *
+ * @brief Hessian operator
+ *
+ * This Hessian operator stores the second derivatives of Orbitals in an OrbitalVector. The ordering of the indices is
+ * [xx, yy, zz, yz, xz, xy].
+ */
+class HessianOperator final : public RankOneOperator<6> {
+public:
+    HessianOperator(std::shared_ptr<mrcpp::DerivativeOperator<3>> D1, std::shared_ptr<mrcpp::DerivativeOperator<3>> D2, double prec) {
+        bool imag = false;
+        NablaOperator N1 = NablaOperator(D1, imag);
+        NablaOperator N2 = NablaOperator(D2, imag);
+        N1.setup(prec);
+        N2.setup(prec);
+
+        // Invoke operator= to assign *this operator
+        RankOneOperator<6> &d = (*this);
+        d[0] = N2[0];
+        d[1] = N2[1];
+        d[2] = N2[2];
+        d[3] = N1[1] * N1[2];
+        d[4] = N1[0] * N1[2];
+        d[5] = N1[0] * N1[1];
+
+        d[0].name() = "del[x]del[x]";
+        d[1].name() = "del[y]del[y]";
+        d[2].name() = "del[z]del[z]";
+        d[3].name() = "del[y]del[z]";
+        d[4].name() = "del[x]del[z]";
+        d[5].name() = "del[x]del[y]";
+    }
+};
+
+}

src/qmoperators/two_electron/XCOperator.h

@@ -70,8 +70,10 @@ class XCOperator final : public RankZeroOperator {
     }
     void clearSpin() { this->potential->setReal(nullptr); }
 
+    std::shared_ptr<XCPotential> getPotential() { return potential; }
+
 private:
     std::shared_ptr<XCPotential> potential{nullptr};
 };
 
-} // namespace mrchem
+} // namespace mrchem