Overload get_dispersion for atom-wise energies (#41)

-Also added a typedef for TRVector

---------

Co-authored-by: Hagen Neugebauer <[email protected]>

Author: haneug

include/dftd_dispersion.h

@@ -71,6 +71,34 @@ extern int get_dispersion(
   double *GRAD
 );
 
+/**
+ * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
+ *
+ * This function calculates the atom-wise dispersion energy and gradients
+ * for the given molecular geometry, considering only the atoms specified
+ * in `realIdx`.
+ *
+ * @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.
+ * @param cutoff Real-space cutoffs for CN and dispersion.
+ * @param energies atom-wise dispersion energies (inout).
+ * @param GRAD Dispersion gradient (inout).
+ * @return Exit status.
+ */
+extern int get_dispersion(
+  const TMolecule &mol,
+  const TIVector &realIdx,
+  int charge,
+  const TD4Model &d4,
+  const dparam &par,
+  TCutoff cutoff,
+  TRVector &energies,
+  double *GRAD
+);
+
 /**
  * @brief Wrapper to handle the evaluation of dispersion energy and derivatives.
  *

include/dftd_matrix.h

@@ -231,5 +231,6 @@ template <class T> class TMatrix {
 };
 
 typedef TVector<int> TIVector;
+typedef TVector<double> TRVector;
 
 } // namespace dftd4

src/dftd_dispersion.cpp

@@ -59,6 +59,37 @@ int get_dispersion(
   double &energy,
   double *GRAD
 ) {
+
+  int info{0};
+
+  int nat = realIdx.Max() + 1;
+
+  TRVector energies; // atom-wise energies
+  energies.NewVector(nat);
+
+  info = get_dispersion(mol, realIdx, charge, d4, par, cutoff, energies, GRAD);
+  if (info != EXIT_SUCCESS) return info;
+
+  // sum up atom-wise energies
+  for (int i = 0; i != nat; i++) {
+    energy += energies(i);
+  }
+
+  energies.DelVec();
+
+  return EXIT_SUCCESS;
+}
+
+int get_dispersion(
+  const TMolecule &mol,
+  const TIVector &realIdx,
+  const int charge,
+  const TD4Model &d4,
+  const dparam &par,
+  const TCutoff cutoff,
+  TRVector &energies,
+  double *GRAD
+) {
   // setup variables
   int info{0};
   bool lmbd = (par.s9 != 0.0);
@@ -137,8 +168,6 @@ int get_dispersion(
 
   TVector<double> dEdcn;
   TVector<double> dEdq;
-  TVector<double> energies;
-  energies.NewVector(nat);
   if (lgrad) {
     dEdcn.NewVector(nat);
     dEdq.NewVector(nat);
@@ -242,12 +271,6 @@ int get_dispersion(
   dEdcn.DelVec();
   dEdq.DelVec();
 
-  // sum up atom-wise energies
-  for (int i = 0; i != nat; i++) {
-    energy += energies(i);
-  }
-  energies.DelVec();
-
   // write to input gradient
   if (lgrad) {
     for (int i = 0, ii = 0; i != mol.NAtoms; i++) {