Real-space cutoffs

include/dftd_cutoff.h

@@ -0,0 +1,51 @@
+/* This file is part of cpp-d4.
+ *
+ * Copyright (C) 2019 Sebastian Ehlert, Marvin Friede
+ *
+ * cpp-d4 is free software: you can redistribute it and/or modify it under
+ * the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * cpp-d4 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with cpp-d4.  If not, see <https://www.gnu.org/licenses/>.
+ */
+#pragma once
+
+namespace dftd {
+
+// Real space cutoff for CN within D4
+static const double cn_default = 30.0;
+
+// Real space cutoff for CN within EEQ
+static const double cn_eeq_default = 25.0;
+
+// Two-body interaction cutoff
+static const double disp2_default = 60.0;
+
+// Three-body interaction cutoff
+static const double disp3_default = 40.0;
+
+
+// Collection of real space cutoffs.
+class TCutoff {
+  public:
+    double disp2;
+    double disp3;
+    double cn;
+    double cn_eeq;
+
+    explicit TCutoff(
+      double cut_disp2 = disp2_default,
+      double cut_disp3 = disp3_default,
+      double cut_cn = cn_default,
+      double cut_cn_eeq = cn_eeq_default
+    );
+};
+
+} // namespace dftd

include/dftd_dispersion.h

@@ -21,6 +21,7 @@
  */
 #pragma once
 
+#include "dftd_cutoff.h"
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
@@ -45,13 +46,15 @@ int d4(const TMolecule& mol, int ndim, double wf, double g_a, double g_c,
        const TVector<double>& cn, TVector<double>& gw, TMatrix<double>& c6ref);
 
 extern
-int edisp(const TMolecule& mol, const TMatrix<double>& dist, const dparam& par, 
-          int ndim, TVector<double>& q,double g_a, double g_c,
-          TVector<double>& gw, TMatrix<double>& c6ref, bool lmbd,
-          double& energy);
+int edisp(const TMolecule& mol, const TMatrix<double>& dist,
+          const double cutoff_disp2, const double cutoff_disp3,
+          const dparam& par, int ndim, TVector<double>& q,
+          double g_a, double g_c, TVector<double>& gw,
+          TMatrix<double>& c6ref, bool lmbd, double& energy);
 
 extern
 int dispgrad(const TMolecule& mol, const TMatrix<double>& dist,
+             const double cutoff_disp2, const double cutoff_disp3,
              const dparam& par, int ndim, const TVector<double>& q, 
              TMatrix<double>& dqdr, TVector<double>& cn, 
              TMatrix<double>& dcndr, double wf, double g_a, double g_c, 
@@ -60,17 +63,28 @@ int dispgrad(const TMolecule& mol, const TMatrix<double>& dist,
 
 extern
 int apprabc(const TMolecule& mol, const TMatrix<double>& dist,
-            const dparam& par, int ndim, TVector<double>& c6ab,
-            double& energy);
+            const double cutoff, const dparam& par, int ndim,
+            TVector<double>& c6ab, double& energy);
 
 extern
 int dabcappr(const TMolecule& mol, const TMatrix<double>& dist, 
-             const dparam& par, int ndim, TVector<double>& gw, 
-             TVector<double>& dgw, TMatrix<double>& c6ref, 
+             const double cutoff, const dparam& par, int ndim,
+             TVector<double>& gw, TVector<double>& dgw, TMatrix<double>& c6ref, 
              TVector<double>& dc6dr, TVector<double>& dc6dcn, double& energy);
 
+/**
+ * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
+ * 
+ * @param mol Molecular geometry.
+ * @param par DFT-D4 parameters.
+ * @param charge Charge of the molecule.
+ * @param cutoff Real-space cutoffs for CN and dispersion.
+ * @param energy Dispersion energy.
+ * @param GRAD Dispersion gradient.
+ * @return Exit status.
+ */
 extern
-int DFTVDW_D4(const TMolecule &mol, const dparam &par, const int &charge,
-              double &energy, double *GRAD);
+int get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
+                   TCutoff cutoff, double &energy, double *GRAD);
 
 }  // namespace dftd

include/dftd_ncoord.h

@@ -43,35 +43,34 @@ extern int calc_distances(const TMolecule& mol, TMatrix<double>& dist);
  * @param mol Molecule object.
  * @param dist Distance matrix.
  * @param cn Vector of coordination numbers.
+ * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
  * @param dcndr Derivative of coordination number.
  * @param lgrad Flag for gradient computation.
- * @param thr Real-space cutoff (default: 1600.0).
  * @return Exit status.
 */
 extern int get_ncoord_erf(
   const TMolecule& mol,
   const TMatrix<double>& dist,
+  const double cutoff,
   TVector<double>& cn,
   TMatrix<double>& dcndr,
-  bool lgrad = false,
-  double thr = 1600.0
+  bool lgrad = false
 );
 
 /**
  * Calculate error function coordination number.
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix.
+ * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
  * @param cn Vector of coordination numbers.
- * @param dcndr Derivative of coordination number.
- * @param thr Real-space cutoff (default: 1600.0).
  * @return Exit status.
 */
 extern int ncoord_erf(
   const TMolecule& mol,
   const TMatrix<double>& dist,
-  TVector<double>& cn,
-  double thr = 1600.0
+  const double cutoff,
+  TVector<double>& cn
 );
 
 /**
@@ -80,17 +79,17 @@ extern int ncoord_erf(
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix.
+ * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
  * @param cn Vector of coordination numbers.
  * @param dcndr Derivative of coordination number.
- * @param thr Real-space cutoff (default: 1600.0).
  * @return Exit status.
 */
 extern int dncoord_erf(
   const TMolecule& mol,
   const TMatrix<double>& dist,
+  const double cutoff,
   TVector<double>& cn,
-  TMatrix<double>& dcndr,
-  double thr = 1600.0
+  TMatrix<double>& dcndr
 );
 
 
@@ -99,33 +98,33 @@ extern int dncoord_erf(
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix.
+ * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
  * @param cn Vector of coordination numbers.
- * @param thr Real-space cutoff (default: 1600.0).
  * @return Exit status.
 */
 extern int get_ncoord_d4(
   const TMolecule& mol,
   const TMatrix<double>& dist,
+  const double cutoff,
   TVector<double>& cn,
   TMatrix<double>& dcndr,
-  bool lgrad = false,
-  double thr = 1600.0
+  bool lgrad = false
 );
 
 /**
  * Calculate covalent coordination number for DFT-D4.
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix.
+ * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
  * @param cn Vector of coordination numbers.
- * @param thr Real-space cutoff (default: 1600.0).
  * @return Exit status.
 */
 extern int ncoord_d4(
   const TMolecule& mol,
   const TMatrix<double>& dist,
-  TVector<double>& cn,
-  double thr = 1600.0
+  const double cutoff,
+  TVector<double>& cn
 );
 
 
@@ -135,17 +134,17 @@ extern int ncoord_d4(
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix.
+ * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
  * @param cn Vector of coordination numbers.
  * @param dcndr Derivative of coordination number.
- * @param thr Real-space cutoff (default: 1600.0).
  * @return Exit status.
 */
 extern int dncoord_d4(
   const TMolecule& mol,
   const TMatrix<double>& dist,
+  const double cutoff,
   TVector<double>& cn,
-  TMatrix<double>& dcndr,
-  double thr = 1600.0
+  TMatrix<double>& dcndr
 );
 
 /**
@@ -161,14 +160,14 @@ extern double erf_count(double k, double rr);
  * Derivative of the counting function w.r.t. the distance.
  * 
  * @param k Steepness of the counting function.
- * @param rr Cutoff radius.
+ * @param rr TCutoff radius.
  * @return Derivative of the counting function.
  */
 extern double derf_count(double k, double rr);
 
 
 /**
- * Cutoff function for large coordination numbers
+ * TCutoff function for large coordination numbers
  * 
  * @param cn_max Maximum CN (not strictly obeyed).
  * @param cn On input coordination number, on output modified CN.

meson.build

@@ -35,10 +35,11 @@ foreach lib: ['cblas', 'lapacke']
 endforeach
 
 srcs = files(
+  'src/dftd_cutoff.cpp',
+  'src/dftd_damping.cpp',
   'src/dftd_dispersion.cpp',
   'src/dftd_eeq.cpp',
   'src/dftd_ncoord.cpp',
-  'src/dftd_damping.cpp',
 )
 
 cpp_d4_inc = include_directories('include')
@@ -53,13 +54,14 @@ cpp_d4_lib = library(
 )
 
 install_headers(
+  'include/dftd_cutoff.h',
+  'include/dftd_damping.h',
   'include/dftd_dispersion.h',
   'include/dftd_eeq.h',
+  'include/dftd_geometry.h',
+  'include/dftd_matrix.h',
   'include/dftd_ncoord.h',
   'include/dftd_parameters.h',
-  'include/dftd_damping.h',
-  'include/dftd_matrix.h',
-  'include/dftd_geometry.h',
 )
 
 cpp_d4_dep = declare_dependency(

src/dftd_cutoff.cpp

@@ -0,0 +1,34 @@
+/* This file is part of cpp-d4.
+ *
+ * Copyright (C) 2019 Sebastian Ehlert, Marvin Friede
+ *
+ * cpp-d4 is free software: you can redistribute it and/or modify it under
+ * the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * cpp-d4 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with cpp-d4.  If not, see <https://www.gnu.org/licenses/>.
+ */
+#include "dftd_cutoff.h"
+
+namespace dftd {
+
+TCutoff::TCutoff(
+  double cut_disp2/* = disp2_default*/,
+  double cut_disp3/* = disp3_default*/,
+  double cut_cn/* = cn_default*/,
+  double cut_cn_eeq/* = cn_eeq_default*/   
+) {
+  disp2 = cut_disp2;
+  disp3 = cut_disp3;
+  cn = cut_cn;
+  cn_eeq = cut_cn_eeq;
+}
+
+} // namespace dftd

src/dftd_damping.cpp

@@ -15,13 +15,11 @@
  * You should have received a copy of the GNU Lesser General Public License
  * along with cpp-d4.  If not, see <https://www.gnu.org/licenses/>.
  */
-
-#include "dftd_damping.h"
-
 #include <algorithm>
 #include <map>
 #include <string>
 
+#include "dftd_damping.h"
 #include "dftd_dispersion.h"
 
 namespace dftd {
@@ -405,8 +403,8 @@ dfunc get_dfunc(std::string name) {
   return none;
 };
 
-dftd::dparam get_d4eeqbjatm_2019_parameter(dfunc num) {
-  dftd::dparam par;
+dparam get_d4eeqbjatm_2019_parameter(dfunc num) {
+  dparam par;
   double s6{0.0}, s8{0.0}, s10{0.0}, s9{1.0}, a1{0.0}, a2{0.0};
   int alp{16};
   switch (num) {
@@ -1116,8 +1114,8 @@ dftd::dparam get_d4eeqbjatm_2019_parameter(dfunc num) {
   return par;
 }
 
-dftd::dparam get_d4eeqbjmbd_2019_parameter(dfunc num) {
-  dftd::dparam par;
+dparam get_d4eeqbjmbd_2019_parameter(dfunc num) {
+  dparam par;
   double s6{0.0}, s8{0.0}, s10{0.0}, s9{1.0}, a1{0.0}, a2{0.0};
   int alp{16};
   switch (num) {
@@ -1661,7 +1659,7 @@ dftd::dparam get_d4eeqbjmbd_2019_parameter(dfunc num) {
   return par;
 }
 
-void d4par(const std::string func, dftd::dparam& par, const bool lmbd) {
+void d4par(const std::string func, dparam& par, const bool lmbd) {
   auto num = get_dfunc(func);
   if (lmbd) {
     par = get_d4eeqbjatm_2019_parameter(num);