refactoring: Create class for coordination number calculation

include/dftd_cblas.h

@@ -41,8 +41,8 @@ inline int BLAS_Add_Mat_x_Vec(
   bool Transpose,
   double alpha
 ) {
-  if (A.rows == C.N && A.cols == V.N) {
     if (Transpose) {
+      if (A.cols == C.N && A.rows == V.N) {
       cblas_dgemv(
         CblasRowMajor,
         CblasTrans,
@@ -58,7 +58,9 @@ inline int BLAS_Add_Mat_x_Vec(
         1
       );
       return EXIT_SUCCESS;
-    } else {
+    };
+  } else {
+      if (A.rows == C.N && A.cols == V.N) {
       cblas_dgemv(
         CblasRowMajor,
         CblasNoTrans,

include/dftd_ncoord.h

@@ -28,6 +28,107 @@
 
 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>&);
+
+    // Electronegativity factor
+    virtual double get_en_factor(int, int) const;
+    // Counting function
+    virtual double count_fct(double) const = 0;
+    // Derivative of the counting function
+    virtual double dcount_fct(double) const = 0;
+    // Soft maximum/cutoff for coordination number
+    virtual int cut_coordination_number(const double, TVector<double>&, TMatrix<double>&, bool) = 0;
+    // Constructor
+    NCoordBase(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0);
+    // Virtual destructor
+    virtual ~NCoordBase() {
+      cn.DelVec();
+      dcndr.DelMat();
+    }
+};
+
+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;
+    // Soft maximum/cutoff for coordination number
+    int cut_coordination_number(const double, TVector<double>&, TMatrix<double>&, bool) 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){}
+    // Use default destructor; base class handles cleanup
+    ~NCoordErf() override = default;
+};
+
+class NCoordErfD4 : 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;
+    // coordination number scaling factor based on electronegativity difference
+    double get_en_factor(int, int) const override;
+    // Soft maximum/cutoff for coordination number
+    int cut_coordination_number(const double, TVector<double>&, TMatrix<double>&, bool) override;
+    // Constructor
+    NCoordErfD4(double optional_kcn = 7.5, double optional_norm_exp = 1.0, double optional_cutoff = 25.0)
+    : NCoordBase(optional_kcn, optional_norm_exp, optional_cutoff){}
+    // Use default destructor; base class handles cleanup
+    ~NCoordErfD4() override = default;
+};
+
 /**
  * Calculate all distance pairs and store in matrix.
  *
@@ -138,86 +239,6 @@ extern int dncoord_erf(
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 
-/**
- * Wrapper for error function coordination number for DFT-D4.
- *
- * @param mol Molecule object.
- * @param realIdx List for real atoms excluding ghost/non atoms.
- * @param dist Distance matrix.
- * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
- * @param cn Vector of coordination numbers.
- * @return Exit status.
- */
-extern int get_ncoord_d4(
-  const TMolecule &mol,
-  const TIVector &realIdx,
-  const TMatrix<double> &dist,
-  double cutoff,
-  TVector<double> &cn,
-  TMatrix<double> &dcndr,
-  bool lgrad = false
-);
-
-/**
- * Wrapper for error function 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.
- * @return Exit status.
- */
-extern int get_ncoord_d4(
-  const TMolecule &mol,
-  const TMatrix<double> &dist,
-  double cutoff,
-  TVector<double> &cn,
-  TMatrix<double> &dcndr,
-  bool lgrad = false
-);
-
-/**
- * Calculate covalent coordination number for DFT-D4.
- *
- * @param mol Molecule object.
- * @param realIdx List for real atoms excluding ghost/non atoms.
- * @param dist Distance matrix.
- * @param cutoff Real-space cutoff (default: @see {dftd_cutoff.h}).
- * @param cn Vector of coordination numbers.
- * @return Exit status.
- */
-extern int ncoord_d4(
-  const TMolecule &mol,
-  const TIVector &realIdx,
-  const TMatrix<double> &dist,
-  double cutoff,
-  TVector<double> &cn
-);
-
-/**
- * Calculate covalent coordination number and derivative
- * w.r.t. nuclear coordinates
- *
- * @param mol Molecule object.
- * @param realIdx List for real atoms excluding ghost/non atoms.
- * @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.
- * @return Exit status.
- */
-extern int dncoord_d4(
-  const TMolecule &mol,
-  const TIVector &realIdx,
-  const TMatrix<double> &dist,
-  double cutoff,
-  TVector<double> &cn,
-  TMatrix<double> &dcndr
-);
-
-///////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////////////////////
-
 /**
  * Error function counting function for coordination number contributions.
  *

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
+  NCoordErfD4 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(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, true, 1.0);
 
-  dcndr.DelMat();
   dEdcn.DelVec();
   dEdq.DelVec();
 

src/dftd_eeq.cpp

@@ -72,24 +72,17 @@ 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();
-
   return EXIT_SUCCESS;
 };
 
@@ -194,9 +187,9 @@ int ChargeModel::eeq_chrgeq(
         jj = realIdx(j);
         if (jj < 0) continue;
 
-        dAmat(3 * jj, ii) -= dcndr(jj, 3 * ii) * dxdcn(ii);
-        dAmat(3 * jj + 1, ii) -= dcndr(jj, 3 * ii + 1) * dxdcn(ii);
-        dAmat(3 * jj + 2, ii) -= dcndr(jj, 3 * ii + 2) * dxdcn(ii);
+        dAmat(3 * jj,     ii) -= dcndr(ii, 3 * jj    ) * dxdcn(ii);
+        dAmat(3 * jj + 1, ii) -= dcndr(ii, 3 * jj + 1) * dxdcn(ii);
+        dAmat(3 * jj + 2, ii) -= dcndr(ii, 3 * jj + 2) * dxdcn(ii);
       }
     }
 
@@ -364,10 +357,10 @@ int EEQModel::get_damat_0d(
       atrace(jj, 1) -= dgy * q(ii);
       atrace(jj, 2) -= dgz * q(ii);
 
-      dAmat(3 * ii, jj) = dgx * q(ii);
-      dAmat(3 * ii + 1, jj) = dgy * q(ii);
-      dAmat(3 * ii + 2, jj) = dgz * q(ii);
-      dAmat(3 * jj, ii) = -dgx * q(jj);
+      dAmat(3 * ii,     jj) =  dgx * q(ii);
+      dAmat(3 * ii + 1, jj) =  dgy * q(ii);
+      dAmat(3 * ii + 2, jj) =  dgz * q(ii);
+      dAmat(3 * jj,     ii) = -dgx * q(jj);
       dAmat(3 * jj + 1, ii) = -dgy * q(jj);
       dAmat(3 * jj + 2, ii) = -dgz * q(jj);
     }