Merge ORCA branch

Author: marvinfriede

Diff for: .gitignore

@@ -35,3 +35,8 @@
 /build*/
 /install*/
 /_*/
+
+# project-specific
+/scripts/old/*
+/scripts/*.h
+/scripts/*.cpp

Diff for: .vscode/settings.json

@@ -0,0 +1,10 @@
+{
+  "files.associations": {
+    "cassert": "cpp",
+    "cmath": "cpp",
+    "complex": "cpp",
+    "functional": "cpp",
+    "iostream": "cpp",
+    "new": "cpp"
+  }
+}

Diff for: README.md

@@ -1,11 +1,13 @@
-# C++ Port of DFT-D4
+# C++ Port of DFT-D4 (for ORCA)
 
 [![License](https://img.shields.io/github/license/dftd4/cpp-d4)](https://github.com/dftd4/cpp-d4/blob/master/COPYING)
 [![Latest Version](https://img.shields.io/github/v/release/dftd4/cpp-d4)](https://github.com/dftd4/cpp-d4/releases/latest)
 
 This project is a port of the [`dftd4`](https://github.com/dftd4/dftd4) project
 to C++ and provides the D4(EEQ)-ATM method.
 
+**NOTE:** This branch contains some adaptations for compatibility with the ORCA Quantum Chemistry program. This inludes some renaming of variables and functions, and, most notably, the handling of ghost atoms.
+
 ## Building This Project
 
 This project is build with `meson`, to setup and perform a build run:

Diff for: app/main.cpp

@@ -163,6 +163,15 @@ int main(int argc, char **argv) {
   dftd4::TCutoff cutoff;
   dftd4::TD4Model d4;
 
+  // masking (nothing excluded)
+  dftd4::TVector<int> realIdx;
+  realIdx.NewVec(mol.NAtoms);
+  int nat = 0;
+  for (int i = 0; i != mol.NAtoms; i++) {
+    realIdx(i) = nat;
+    nat++;
+  }
+
   // analytical gradient
   double *d4grad;
   if (lgrad) {
@@ -174,7 +183,9 @@ int main(int argc, char **argv) {
     d4grad = nullptr;
   }
 
-  info = dftd4::get_dispersion(mol, charge, d4, par, cutoff, energy, d4grad);
+  info = dftd4::get_dispersion(
+    mol, realIdx, charge, d4, par, cutoff, energy, d4grad
+  );
   if (info != EXIT_SUCCESS) return info;
 
   // Print results
@@ -196,6 +207,7 @@ int main(int argc, char **argv) {
   }
 
   mol.FreeMemory();
+  delete[] d4grad;
 
   return EXIT_SUCCESS;
 }

Diff for: app/readxyz.cpp

@@ -56,11 +56,11 @@ void read_xyzfile(const std::string &name, dftd4::TMolecule &mol) {
       printf("Error: Input file could not be read.");
       exit(EXIT_FAILURE);
     }
-    mol.xyz(i, 0) = x * aatoau;
-    mol.xyz(i, 1) = y * aatoau;
-    mol.xyz(i, 2) = z * aatoau;
+    mol.CC(i, 0) = x * aatoau;
+    mol.CC(i, 1) = y * aatoau;
+    mol.CC(i, 2) = z * aatoau;
     at = element(sym);
-    mol.at(i) = at;
+    mol.ATNO(i) = at;
   }
 
   geo.close();

Diff for: include/damping/dftd_atm.h

@@ -29,6 +29,7 @@ namespace dftd4 {
 
 extern int get_atm_dispersion(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   double s9,
@@ -47,6 +48,7 @@ extern int get_atm_dispersion(
 
 extern int get_atm_dispersion_energy(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   double s9,
@@ -59,6 +61,7 @@ extern int get_atm_dispersion_energy(
 
 extern int get_atm_dispersion_derivs(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   double s9,

Diff for: include/damping/dftd_rational.h

@@ -31,6 +31,7 @@ extern inline double fdmprdr_bj(int n, double r, double c);
 
 extern int get_dispersion2(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   const dparam &par,
@@ -46,6 +47,7 @@ extern int get_dispersion2(
 
 extern int get_dispersion3(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   const dparam &par,
@@ -63,6 +65,7 @@ extern int get_dispersion3(
 
 extern int get_dispersion2_energy(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   const dparam &par,
@@ -72,6 +75,7 @@ extern int get_dispersion2_energy(
 
 extern int get_dispersion2_derivs(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   double cutoff,
   const dparam &par,

Diff for: include/dftd_dispersion.h

@@ -42,6 +42,7 @@ class dparam {
  * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
  *
  * @param mol Molecular geometry.
+ * @param realIdx List for real atoms excluding ghost/non atoms
  * @param charge Molecular charge.
  * @param par DFT-D4 parameters.
  * @param d4 Base D4 dispersion model.
@@ -52,6 +53,7 @@ class dparam {
  */
 extern int get_dispersion(
   const TMolecule &mol,
+  const TIVector &realIdx,
   int charge,
   const TD4Model &d4,
   const dparam &par,

Diff for: include/dftd_eeq.h

@@ -29,6 +29,7 @@ namespace dftd4 {
 
 extern int get_charges(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   int charge,
   double cutoff,
@@ -39,6 +40,7 @@ extern int get_charges(
 
 extern int get_vrhs(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const int &charge,
   const TVector<double> &cn,
   TVector<double> &Xvec,
@@ -48,12 +50,14 @@ extern int get_vrhs(
 
 extern int get_amat_0d(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   TMatrix<double> &Amat
 );
 
 extern int get_damat_0d(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   const TVector<double> &q,
   const TMatrix<double> &Amat,
@@ -63,6 +67,7 @@ extern int get_damat_0d(
 
 extern int eeq_chrgeq(
   const TMolecule &mol,
+  const TIVector &realIdx,
   const TMatrix<double> &dist,
   const int &charge,
   const TVector<double> &cn,

Diff for: include/dftd_geometry.h

@@ -25,22 +25,22 @@ namespace dftd4 {
 class TMolecule {
   public:
     int NAtoms;
-    TMatrix<double> xyz; // Cartesian Coordinates: (NAtoms x 3)-matrix
-    TVector<int> at;     // atomic numbers
+    TMatrix<double> CC; // Cartesian Coordinates: (NAtoms x 3)-matrix
+    TVector<int> ATNO;  // atomic numbers
 
     TMolecule() { NAtoms = 0; }
     ~TMolecule() { FreeMemory(); }
 
     void GetMemory(int NumAt_) {
       FreeMemory();
       NAtoms = NumAt_;
-      xyz.New(NAtoms, 3);
-      at.New(NAtoms);
+      CC.New(NAtoms, 3);
+      ATNO.New(NAtoms);
     }
 
-    void FreeMemory() {
-      xyz.Delete();
-      at.Delete();
+    void FreeMemory(void) {
+      CC.Delete();
+      ATNO.Delete();
     }
 };
 

Diff for: include/dftd_matrix.h

@@ -63,6 +63,8 @@ template <class T> class TVector {
       p = nullptr;
       N = 0;
     }
+    void DelVec(void) { return Delete(); }
+
     void CopyVec(const TVector &v) {
       long int mem;
       if (N != v.N) {
@@ -73,25 +75,41 @@ template <class T> class TVector {
       mem = (long int)N * ElementSize;
       std::memcpy(p, v.p, mem);
     }
+
     void Init() {
       if (p != nullptr) {
         long int mem = (long int)N * ElementSize;
         std::memset(p, 0, mem);
       }
     }
 
-    void Print(char name[]) {
+    void Print(const char name[]) {
       printf("Vector printed: %s (%d)\n", name, N);
       for (int i = 0; i < N; i++) {
         printf("%+23.15e\n", p[i]);
       }
       printf("\n");
     }
+    void PrintInt(const char name[]) {
+      printf("Vector printed: %s (%d)\n", name, N);
+      for (int i = 0; i < N; i++) {
+        printf("%d\n", p[i]);
+      }
+      printf("\n");
+    }
 
     inline T &operator()(int i) { return p[i]; }
     inline const T &operator()(int i) const { return p[i]; }
     inline T &operator[](int i) { return p[i]; }
     inline const T &operator[](int i) const { return p[i]; }
+
+    // Max element
+    T Max() const {
+      T result = p[0];
+      for (int i = 0; i < N; ++i)
+        if (result < p[i]) result = p[i];
+      return result;
+    };
 };
 
 // Define a normal matrix
@@ -142,6 +160,7 @@ template <class T> class TMatrix {
       cols = 0;
       p = nullptr;
     }
+    void DelMat(void) { return Delete(); }
 
     void Init() {
       long int mem;
@@ -211,4 +230,6 @@ template <class T> class TMatrix {
     inline T *operator[](int i) { return p + i * cols; }
 };
 
+typedef TVector<int> TIVector;
+
 } // namespace dftd4

Diff for: include/dftd_model.h

@@ -34,9 +34,9 @@ static const double wf_default = 6.0;
 
 class TD4Model {
   public:
+    double wf;
     double ga;
     double gc;
-    double wf;
 
     explicit TD4Model(
       double ga_scale = ga_default,
@@ -48,6 +48,7 @@ class TD4Model {
 
     int weight_references(
       const TMolecule &mol,
+      const TIVector &realIdx,
       const TVector<double> &cn,
       const TVector<double> &q,
       const TMatrix<double> &refq,
@@ -59,6 +60,7 @@ class TD4Model {
 
     int get_atomic_c6(
       const TMolecule &mol,
+      const TIVector &realIdx,
       const TMatrix<double> &gwvec,
       const TMatrix<double> &dgwdcn,
       const TMatrix<double> &dgwdq,
@@ -68,10 +70,17 @@ class TD4Model {
       bool lgrad = false
     ) const;
 
-    virtual int set_refq_eeq(const TMolecule &mol, TMatrix<double> &refq) const;
+    virtual int set_refq_eeq(
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      TMatrix<double> &refq
+    ) const;
 
-    virtual int
-      set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha) const;
+    virtual int set_refalpha_eeq(
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      TMatrix<double> &alpha
+    ) const;
 };
 
 extern inline double trapzd(const double a[23], const double b[23]);