add new radius_factor keyword.

MRChemSoft · Aug 7, 2024 · 4992efd · 4992efd
1 parent 6c75108
commit 4992efd
Show file tree

Hide file tree

Showing 12 changed files with 39 additions and 10 deletions.
diff --git a/doc/users/user_ref.rst b/doc/users/user_ref.rst
@@ -467,6 +467,12 @@ User input reference
     **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`

diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py
@@ -296,6 +296,7 @@ def write_scf_properties(user_dict, origin):
             "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"],
         }
     return prop_dict
 

diff --git a/python/mrchem/input_parser/README.md b/python/mrchem/input_parser/README.md
@@ -1,2 +1,2 @@
-This file was automatically generated by parselglossy on 2024-08-01
+This file was automatically generated by parselglossy on 2024-08-07
 Editing is *STRONGLY DISCOURAGED*
diff --git a/python/mrchem/input_parser/__init__.py b/python/mrchem/input_parser/__init__.py
@@ -1,4 +1,4 @@
 # -*- coding: utf-8 -*-
 
-# This file was automatically generated by parselglossy on 2024-08-01
+# This file was automatically generated by parselglossy on 2024-08-07
 # Editing is *STRONGLY DISCOURAGED*
diff --git a/python/mrchem/input_parser/api.py b/python/mrchem/input_parser/api.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# This file was automatically generated by parselglossy on 2024-08-01
+# This file was automatically generated by parselglossy on 2024-08-07
 # Editing is *STRONGLY DISCOURAGED*
 
 from copy import deepcopy
@@ -376,7 +376,10 @@ def stencil() -> JSONDict:
                                                               "in ['low', "
                                                               "'medium', "
                                                               "'high']"],
-                                            'type': 'str'}],
+                                            'type': 'str'},
+                                        {   'default': 0.6,
+                                            'name': 'radius_factor',
+                                            'type': 'float'}],
                         'name': 'Forces'},
                     {   'keywords': [   {   'default': [],
                                             'name': 'electric_field',

diff --git a/python/mrchem/input_parser/cli.py b/python/mrchem/input_parser/cli.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# This file was automatically generated by parselglossy on 2024-08-01
+# This file was automatically generated by parselglossy on 2024-08-07
 # Editing is *STRONGLY DISCOURAGED*
 
 import argparse

diff --git a/python/mrchem/input_parser/docs/user_ref.rst b/python/mrchem/input_parser/docs/user_ref.rst
@@ -467,6 +467,12 @@ User input reference
     **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`

diff --git a/python/mrchem/input_parser/plumbing/lexer.py b/python/mrchem/input_parser/plumbing/lexer.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# This file was automatically generated by parselglossy on 2024-08-01
+# This file was automatically generated by parselglossy on 2024-08-07
 # Editing is *STRONGLY DISCOURAGED*
 
 import json

diff --git a/python/template.yml b/python/template.yml
@@ -424,6 +424,13 @@ sections:
         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.

diff --git a/src/driver.cpp b/src/driver.cpp
@@ -548,7 +548,8 @@ void driver::scf::calc_properties(const json &json_prop, Molecule &mol, const js
             Timer t_surface;
             t_surface.start();
             std::string leb_prec = json_prop["geometric_derivative"]["geom-1"]["surface_integral_precision"];
-            Eigen::MatrixXd surfaceForces = surface_force::surface_forces(mol, Phi, prec, json_fock, leb_prec);
+            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);
             t_surface.stop();
             GeometricDerivative &G = mol.getGeometricDerivative("geom-1");
             auto &nuc = G.getNuclear();

diff --git a/src/surface_forces/SurfaceForce.cpp b/src/surface_forces/SurfaceForce.cpp
@@ -253,7 +253,12 @@ VectorXd distanceToNearestNeighbour(MatrixXd pos){
  * @return The matrix of forces, shape (nAtoms, 3).
  */
 Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi, double prec, const json &json_fock
-        , std::string leb_prec) {
+        , std::string leb_prec, double radius_factor) {
+
+    if (radius_factor > 0.95 && radius_factor < 0.05){
+        MSG_ABORT("Invalid value of radius_factor")
+    }
+    std::cout << "Radius factor " << radius_factor << std::endl;
 
     // setup density
     mrchem::Density rho(false);
@@ -357,7 +362,7 @@ Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi
     Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(numAtoms, 3);
 
     for (int iAtom = 0; iAtom < numAtoms; iAtom++) {
-        radius = dist(iAtom) * .6;
+        radius = dist(iAtom) * radius_factor;
         coord = mol.getNuclei()[iAtom].getCoord();
         center << coord[0], coord[1], coord[2];
 

diff --git a/src/surface_forces/SurfaceForce.h b/src/surface_forces/SurfaceForce.h
@@ -12,7 +12,7 @@ namespace surface_force {
 
 // Function declaration
 Eigen::MatrixXd surface_forces(mrchem::Molecule &mol, mrchem::OrbitalVector &Phi, double prec
-    , const json &json_fock, std::string leb_prec);
+    , const json &json_fock, std::string leb_prec, double radius_factor);
 
 } // namespace surface_force