Refactor and adapt some ORCA conventions

include/dftd_cblas.h

@@ -24,7 +24,7 @@ extern "C" {
 #include "lapacke.h"
 }
 
-namespace dftd {
+namespace dftd4 {
 
 inline int BLAS_Add_Mat_x_Vec(TVector<double>& C, TMatrix<double>& A,
                               TVector<double>& V, bool Transpose,
@@ -167,4 +167,4 @@ inline int BLAS_InvertMatrix(TMatrix<double>& a) {
   return EXIT_SUCCESS;
 };
 
-}  // namespace dftd
+}  // namespace dftd4

include/dftd_cutoff.h

@@ -17,7 +17,7 @@
  */
 #pragma once
 
-namespace dftd {
+namespace dftd4 {
 
 // Collection of real space cutoffs.
 class TCutoff {
@@ -36,4 +36,4 @@ class TCutoff {
 
 };
 
-}; // namespace dftd
+}; // namespace dftd4

include/dftd_damping.h

@@ -26,7 +26,7 @@
 
 #include "dftd_dispersion.h"
 
-namespace dftd {
+namespace dftd4 {
 
 /**
  * @brief Collect the D4 parameters.
@@ -38,7 +38,7 @@ namespace dftd {
  */
 extern int d4par(
   const std::string func,
-  dftd::dparam &par,
+  dftd4::dparam &par,
   const bool latm = true
 );
 

include/dftd_dispersion.h

@@ -25,11 +25,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
-
-extern inline double fdmpr_bj(const int n, const double r, const double c);
-extern inline double fdmprdr_bj(const int n, const double r, const double c);
-
+namespace dftd4 {
 
 class dparam {
  public:
@@ -42,6 +38,26 @@ class dparam {
   int alp;
 };
 
+/**
+ * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
+ * 
+ * @param Dat Molecular information (only for consistency with ORCA).
+ * @param mol Molecular geometry.
+ * @param par DFT-D4 parameters.
+ * @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 TMolInfo &Dat,
+  const TMolecule &mol,
+  const dparam &par,
+  TCutoff cutoff,
+  double &energy,
+  double *GRAD
+);
+
 // Generic wrappers for two- and three-body dispersion
 
 extern int get_dispersion2(
@@ -150,19 +166,7 @@ extern int get_atm_dispersion_derivs(
 
 extern double triple_scale(int ii, int jj, int kk);
 
-/**
- * @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 get_dispersion(const TMolecule &mol, const dparam &par, const int &charge,
-                   TCutoff cutoff, double &energy, double *GRAD);
+extern inline double fdmpr_bj(const int n, const double r, const double c);
+extern inline double fdmprdr_bj(const int n, const double r, const double c);
 
-}  // namespace dftd
+}  // namespace dftd4

include/dftd_eeq.h

@@ -25,7 +25,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 extern int get_charges(
   const TMolecule& mol,

include/dftd_geometry.h

@@ -20,7 +20,7 @@
 
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 // Input of the molecular geometry
 class TMolecule {
  public:
@@ -44,4 +44,13 @@ class TMolecule {
   }
 };
 
-}  // namespace dftd
+class TMolInfo {
+  public:
+    int Charge;
+
+    TMolInfo(int c = 0) {
+      Charge = c;
+    }
+};
+
+}  // namespace dftd4

include/dftd_matrix.h

@@ -22,7 +22,7 @@
 #include <cstring>
 #include <iostream>
 
-namespace dftd {
+namespace dftd4 {
 
 // Define a vector
 template <class T>
@@ -40,7 +40,7 @@ class TVector {
   ~TVector() {
     if (p != 0) Delete();
   }
-  void New(int VectorLength) {
+  void NewVector(int VectorLength) {
     if (VectorLength < 0) {
       std::exit(EXIT_FAILURE);
     }
@@ -60,6 +60,11 @@ class TVector {
       Init();
     }
   }
+
+  // alias for NewVector
+  void New(int VectorLength) { return NewVector(VectorLength); }
+  void NewVec(int VectorLength) { return NewVector(VectorLength); }
+
   void Delete(void) {
     if (p != 0 && N != 0) {
       delete[] p;
@@ -116,7 +121,8 @@ class TMatrix {
     if (p != 0) Delete();
   }
 
-  void New(int r, int c) {
+
+  void NewMatrix(int r, int c) {
     if (r < 0 || c < 0) std::exit(EXIT_FAILURE);
     if (p != 0 && r == rows && c == cols) {
       Init();
@@ -135,7 +141,11 @@ class TMatrix {
     return;
   }
 
-  void New(const TMatrix& v) { New(v.rows, v.cols); }
+  void NewMatrix(const TMatrix& v) { NewMatrix(v.rows, v.cols); }
+
+  // alias for NewMatrix
+  void New(int r, int c) { return NewMatrix(r, c); }
+  void NewMat(int r, int c) { return NewMatrix(r, c); }
 
   void Delete(void) {
     if (p != 0 && rows * cols != 0) {
@@ -172,7 +182,7 @@ class TMatrix {
         // for non-square matrix, we need an additional copy
         TMatrix<T> temp;
         temp.CopyMat(*this);
-        New(cols, rows);
+        NewMatrix(cols, rows);
         for (i = 0; i < rows; i++) {
           for (j = 0; j < cols; j++) {
             p[i * cols + j] = temp.p[j * cols + i];
@@ -214,4 +224,4 @@ class TMatrix {
   inline T* operator[](int i) { return p + i * cols; }
 };
 
-}  // namespace dftd
+}  // namespace dftd4

include/dftd_model.h

@@ -21,7 +21,7 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 extern inline double trapzd(const double a[23], const double b[23]);
 
@@ -32,6 +32,8 @@ extern inline double zeta(const double a, const double c, const double qref,
 extern inline double dzeta(const double a, const double c, const double qref,
                     const double qmod);
 
+extern int get_max_ref(const TMolecule &mol, int &mref);
+
 extern int weight_references(
   const TMolecule& mol,
   const TVector<double>& cn,
@@ -57,4 +59,4 @@ extern int set_refalpha_eeq(const TMolecule& mol, TMatrix<double>& alpha);
 
 extern bool is_exceptional(double val);
 
-} // namespace dftd
+} // namespace dftd4

include/dftd_ncoord.h

@@ -26,10 +26,10 @@
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 
-namespace dftd {
+namespace dftd4 {
 
 /**
- * Calculate all distance pairs and store in matrix
+ * Calculate all distance pairs and store in matrix.
  * 
  * @param mol Molecule object.
  * @param dist Distance matrix (inout).
@@ -186,4 +186,4 @@ extern inline double log_cn_cut(const double cn_max, const double cn);
 
 extern inline double dlog_cn_cut(const double cn_max, const double cn);
 
-};  // namespace dftd
+};  // namespace dftd4

include/dftd_parameters.h

@@ -26,7 +26,7 @@
  */
 #pragma once
 
-namespace dftd
+namespace dftd4
 {
 
 static const double zeff[119] { 0,

include/dftd_readxyz.h

@@ -21,6 +21,6 @@
 
 #include "dftd_geometry.h"
 
-extern void read_xyzfile(const std::string&, dftd::TMolecule&);
+extern void read_xyzfile(const std::string&, dftd4::TMolecule&);
 
 extern int element(const std::string&);

src/dftd_cutoff.cpp

@@ -17,7 +17,7 @@
  */
 #include "dftd_cutoff.h"
 
-namespace dftd {
+namespace dftd4 {
 
 // Real space cutoff for CN within D4
 static const double cn_default = 30.0;
@@ -47,4 +47,4 @@ void TCutoff::disable() {
   cn_eeq = new_cutoff;
 };
 
-}; // namespace dftd
+}; // namespace dftd4

src/dftd_damping.cpp

@@ -22,7 +22,7 @@
 #include "dftd_damping.h"
 #include "dftd_dispersion.h"
 
-namespace dftd {
+namespace dftd4 {
 
 enum dfunc {
   none = 0,
@@ -1536,7 +1536,7 @@ dparam get_d4eeqbjmbd_2019_parameter(dfunc num) {
 
 int d4par(
   const std::string func,
-  dftd::dparam &par,
+  dftd4::dparam &par,
   const bool latm/* = true*/
 ) {
   auto num = get_dfunc(func);
@@ -1548,4 +1548,4 @@ int d4par(
   return EXIT_SUCCESS;
 }
 
-}  // namespace dftd
+}  // namespace dftd4