refactoring: Create class for coordination number calculation

include/dftd_ncoord.h

@@ -28,6 +28,79 @@
 
 namespace dftd4 {
 
+
+class NCoordBase
+{
+  public:
+    TVector<double> cn;
+    TMatrix<double> dcndr;
+    static const double rad[];
+    double kcn;  // Steepness of counting function 
+    double norm_exp;
+    double cutoff;
+    // Get the coordination number
+    int get_ncoord( // with ghost atoms
+      const TMolecule&,
+      const TMatrix<double>&,
+      const double,
+      TVector<double>&,
+      TMatrix<double>&,
+      bool);
+    int get_ncoord(  // without ghost atoms
+      const TMolecule&,
+      const TIVector&,
+      const TMatrix<double>&,
+      bool);
+    // Calculate the coordination number using the virtual counting function
+    int ncoord_base(
+      const TMolecule&,
+      const TIVector&,
+      const TMatrix<double>&);
+    // Calculate the derivative of the coordination number
+    int dncoord_base(
+      const TMolecule&,
+      const TIVector&,
+      const TMatrix<double>&);
+    // Get the DFT-D4 coordination number
+    int get_ncoord_d4(
+      const TMolecule&,
+      const TMatrix<double>&,
+      bool);
+    int get_ncoord_d4(
+      const TMolecule&,
+      const TIVector&,
+      const TMatrix<double>&,
+      bool);
+    int ncoord_d4(
+      const TMolecule&,
+      const TIVector&,
+      const TMatrix<double>&);
+    int dncoord_d4(
+      const TMolecule&,
+      const TIVector&,
+      const TMatrix<double>&);
+
+    // Counting function
+    virtual double count_fct(double) const = 0;
+    // Derivative of the counting function
+    virtual double dcount_fct(double) const = 0;
+    // Constructor
+    NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0);
+    // Virtual destructor
+    virtual ~NCoordBase() = default;
+};
+
+class NCoordErf : public NCoordBase {
+  public:
+    // erf() based counting function
+    double count_fct(double) const override;
+    // derivative of the erf() based counting function
+    double dcount_fct(double) const override;
+    // Constructor
+    NCoordErf(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0)
+    : NCoordBase(optional_kcn, optional_norm_exp, optional_cutoff){}
+};
+
 /**
  * Calculate all distance pairs and store in matrix.
  *

src/dftd_dispersion.cpp

@@ -102,17 +102,14 @@ int get_dispersion(
   dist.NewMatrix(nat, nat);
   calc_distances(mol, realIdx, dist);
 
-  TVector<double> cn;       // D4 coordination number
   TVector<double> q;        // partial charges from EEQ model
-  TMatrix<double> dcndr;    // derivative of D4-CN
   TMatrix<double> dqdr;     // derivative of partial charges
   TVector<double> gradient; // derivative of dispersion energy
+  NCoordErf ncoord_erf;
   multicharge::EEQModel chrg_model;  // Charge model
 
-  cn.NewVector(nat);
   q.NewVector(nat);
   if (lgrad) {
-    dcndr.NewMatrix(3 * nat, nat);
     dqdr.NewMatrix(3 * nat, nat);
     gradient.NewVector(3 * nat);
   }
@@ -122,7 +119,7 @@ int get_dispersion(
   if (info != EXIT_SUCCESS) return info;
 
   // get the D4 coordination number
-  info = get_ncoord_d4(mol, realIdx, dist, cutoff.cn, cn, dcndr, lgrad);
+  info = ncoord_erf.get_ncoord_d4(mol, realIdx, dist, lgrad);
   if (info != EXIT_SUCCESS) return info;
 
   // maximum number of reference systems
@@ -145,7 +142,7 @@ int get_dispersion(
     dgwdq.NewMatrix(mref, nat);
   }
   info = d4.weight_references(
-    mol, realIdx, cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad
+    mol, realIdx, ncoord_erf.cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad
   );
   if (info != EXIT_SUCCESS) return info;
 
@@ -212,11 +209,10 @@ int get_dispersion(
       dgwdq.NewMatrix(mref, nat);
     }
     info = d4.weight_references(
-      mol, realIdx, cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad
+      mol, realIdx, ncoord_erf.cn, q, refq, gwvec, dgwdcn, dgwdq, lgrad
     );
     if (info != EXIT_SUCCESS) return info;
 
-    cn.Delete();
     q.Delete();
     refq.Delete();
 
@@ -253,7 +249,6 @@ int get_dispersion(
     );
     if (info != EXIT_SUCCESS) return info;
   } else {
-    cn.Delete();
     q.Delete();
     refq.Delete();
     gwvec.Delete();
@@ -266,9 +261,8 @@ int get_dispersion(
   dc6dcn.DelMat();
   dc6dq.DelMat();
 
-  if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, dcndr, dEdcn, false, 1.0); }
+  if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, ncoord_erf.dcndr, dEdcn, false, 1.0); }
 
-  dcndr.DelMat();
   dEdcn.DelVec();
   dEdq.DelVec();
 

src/dftd_eeq.cpp

@@ -72,23 +72,19 @@ int ChargeModel::get_charges(
   bool lverbose{false};
   int nat = realIdx.Max() + 1;
 
-  TVector<double> cn;    // EEQ cordination number
-  TMatrix<double> dcndr; // Derivative of EEQ-CN
-
-  cn.NewVec(nat);
-  if (lgrad) dcndr.NewMat(nat, 3 * nat);
+  dftd4::NCoordErf ncoord_erf;
 
   // get the EEQ coordination number
-  info = get_ncoord_erf(mol, realIdx, dist, cutoff, cn, dcndr, lgrad);
+  info = ncoord_erf.get_ncoord(mol, realIdx, dist, lgrad);
   if (info != EXIT_SUCCESS) return info;
 
   // corresponds to model%solve in Fortran
   info =
-    eeq_chrgeq(mol, realIdx, dist, charge, cn, q, dcndr, dqdr, lgrad, lverbose);
+    eeq_chrgeq(mol, realIdx, dist, charge, ncoord_erf.cn, q, ncoord_erf.dcndr, dqdr, lgrad, lverbose);
   if (info != EXIT_SUCCESS) return info;
 
-  dcndr.DelMat();
-  cn.DelVec();
+  ncoord_erf.dcndr.DelMat();
+  ncoord_erf.cn.DelVec();
 
   return EXIT_SUCCESS;
 };