Implementation of the charge gradient dqdr for EEQ-BC with unit test and speed up for solve without gradient

include/dftd_cblas.h

@@ -74,7 +74,7 @@ inline int BLAS_Add_Mat_x_Vec(
   };
 
   return EXIT_FAILURE;
-};
+}
 
 /**
  * @brief General matrix-matrix multiplication (`C = alpha * A * B + C`).
@@ -97,8 +97,9 @@ inline int BLAS_Add_Mat_x_Mat(
 ) {
   // check for size 0 matrices
   if (A.cols == 0 || A.rows == 0 || B.cols == 0 || B.rows == 0 || C.cols == 0 ||
-      C.rows == 0)
-    exit(EXIT_FAILURE);
+      C.rows == 0) {
+      exit(EXIT_FAILURE);
+  };
 
   // check for transpositions
   if (!TransposeA) {
@@ -197,7 +198,7 @@ inline int BLAS_Add_Mat_x_Mat(
     }; // B transposed
   };
   return EXIT_SUCCESS;
-};
+}
 
 /**
  * @brief Compute inverse of a matrix using LU decomposition.
@@ -234,6 +235,51 @@ inline int BLAS_InvertMatrix(TMatrix<double> &a) {
   if (info != 0) { return EXIT_FAILURE; }
 
   return EXIT_SUCCESS;
-};
+}
+
+/**
+ * @brief Solve a symmetric linear system A * X = B for X.
+ *
+ * This routine factorizes a symmetric matrix A using Bunch-Kaufman factorization
+ * and solves for the right-hand side vector B. The matrix A is overwritten
+ * by its factorization. The solution overwrites B.
+ *
+ * @param A Symmetric matrix of size (m x m). Overwritten by the factorization.
+ * @param B Right-hand side vector of size m. Overwritten by the solution.
+ * @return int Returns EXIT_SUCCESS (0) on success, EXIT_FAILURE (1) on error.
+ */
+inline int BLAS_SolveSymmetric(
+  TMatrix<double> &A,   // symmetric matrix
+  TVector<double> &B    // RHS vector (becomes solution)
+) {
+  const lapack_int m    = A.rows;
+  const lapack_int nrhs = 1;
+
+  if (A.cols != m || B.N != m) {
+    fprintf(stderr, "BLAS_SolveSymmetric error: dimension mismatch\n");
+    return EXIT_FAILURE;
+  }
+
+  lapack_int info;
+  lapack_int *ipiv = new lapack_int[A.rows];
+
+  // Factorization
+  info = LAPACKE_dsytrf(LAPACK_ROW_MAJOR, 'L', m, A.p, m, ipiv);
+  if (info != 0) {
+    delete[] ipiv;
+    fprintf(stderr, "dsytrf failed: info=%d\n", (int)info);
+    return EXIT_FAILURE;
+  }
+
+  // Solve for all RHS columns
+  info = LAPACKE_dsytrs(LAPACK_ROW_MAJOR, 'L', m, nrhs, A.p, m, ipiv, B.p, nrhs);
+  delete[] ipiv;
+
+  if (info != 0) {
+    fprintf(stderr, "dsytrs failed: info=%d\n", (int)info);
+    return EXIT_FAILURE;
+  }
+  return EXIT_SUCCESS;
+}
 
 } // namespace dftd4

include/dftd_eeq.h

@@ -134,10 +134,15 @@ class ChargeModel
   * @param charge     Target total charge of the molecule.
   * @param dist       Interatomic distance matrix.
   * @param cn         Coordination numbers for each atom.
-  * @param Xvec       Output vector representing the effective driving term for the
+  * @param dcndr      Derivative of the coordination number w.r.t. atom positions
+  * @param xvec       Output vector representing the effective driving term for the
   *                   EEQ system (size: number of atoms).
-  * @param dXvec      Output vector of derivatives of Xvec w.r.t. atomic positions.
+  * @param dxvecdr    Output vector of derivatives of xvec w.r.t. atom positions.
   *                   Only filled if `lgrad` is true.
+  * @param qloc       Local charge (EEQ-BC specific).
+  * @param dqlocdr    Derivative of the local charge (EEQ-BC specific).
+  * @param cmat       Capacitance matrix.
+  * @param dcmatdr    Derivative of the capacitance matrix.
   * @param lgrad      If true, compute derivatives of Xvec.
   *
   * @returns EXIT_SUCCESS (0) on success, or a nonzero error code if the computation fails.
@@ -148,8 +153,13 @@ class ChargeModel
     const int &charge,
     const TMatrix<double> &dist,
     const TVector<double> &cn,
-    TVector<double> &Xvec,
-    TVector<double> &dXvec,
+    const TMatrix<double> &dcndr,
+    TVector<double> &xvec,
+    TMatrix<double> &dxvecdr,
+    TVector<double> &qloc,
+    TMatrix<double> &dqlocdr,
+    TMatrix<double> &cmat,
+    TMatrix<double> &dcmatdr,
     bool lgrad
   ) = 0;
 
@@ -163,6 +173,8 @@ class ChargeModel
   * @param realIdx Mapping of real atom indices (excludes dummy atoms).
   * @param dist    Interatomic distance matrix.
   * @param cn      Coordination numbers for each atom.
+  * @param qloc    Local charge (EEQ-BC specific)
+  * @param cmat    Capacitance matrix (EEQ-BC specific)
   * @param Amat    Output Coulomb matrix (size: number of atoms + 1 for charge constraint).
   * 
   * @returns EXIT_SUCCESS (0) on success, or a nonzero error code if the computation fails.
@@ -172,6 +184,8 @@ class ChargeModel
     const TIVector &realIdx,
     const TMatrix<double> &dist,
     const TVector<double> &cn,
+    const TVector<double> &qloc,
+    const TMatrix<double> &cmat,
     TMatrix<double> &Amat
   ) const = 0;
 
@@ -185,8 +199,14 @@ class ChargeModel
   * @param dist    Interatomic distance matrix.
   * @param q       Atomic charges computed from the EEQ model.
   * @param Amat    Output coulomb matrix.
-  * @param dAmat   Output derivative matrix with respect to Cartesian coordinates.
+  * @param dAmatdr   Output derivative matrix with respect to Cartesian coordinates.
   * @param atrace  Output trace contributions for the derivative.
+  * @param cn         Coordination numbers for each atom.
+  * @param dcndr      Derivative of the coordination number w.r.t. atom positions
+  * @param qloc       Local charge (EEQ-BC specific).
+  * @param dqlocdr    Derivative of the local charge (EEQ-BC specific).
+  * @param cmat       Capacitance matrix.
+  * @param dcmatdr    Derivative of the capacitance matrix.
   * 
   * @returns EXIT_SUCCESS (0) on success, or a nonzero error code if the computation fails.
   */
@@ -196,8 +216,14 @@ class ChargeModel
     const TMatrix<double> &dist,
     const TVector<double> &q,
     const TMatrix<double> &Amat,
-    TMatrix<double> &dAmat,
-    TMatrix<double> &atrace
+    TMatrix<double> &dAmatdr,
+    TMatrix<double> &atrace,
+    const TVector<double> &cn,
+    const TMatrix<double> &dcndr,
+    const TVector<double> &qloc,
+    const TMatrix<double> &dqlocdr,
+    const TMatrix<double> &cmat,
+    const TMatrix<double> &dcmatdr
   ) const = 0;
 
   /**
@@ -240,16 +266,21 @@ class EEQModel : public ChargeModel {
     // initializes its coordination number function.
     EEQModel();
 
-    // Computes the right-hand side vector (Xvec) for the EEQ charge equations.
+    // Computes the right-hand side vector (xvec) for the EEQ charge equations.
     // Optionally computes its derivatives with respect to atomic positions.
     int get_vrhs(
       const TMolecule &mol,
       const TIVector &realIdx,
       const int &charge,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
-      TVector<double> &Xvec,
-      TVector<double> &dXvec,
+      const TMatrix<double> &dcndr,
+      TVector<double> &xvec,
+      TMatrix<double> &dxvecdr,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      TMatrix<double> &cmat,
+      TMatrix<double> &dcmatdr,
       bool lgrad
     ) override;
 
@@ -259,6 +290,8 @@ class EEQModel : public ChargeModel {
       const TIVector &realIdx,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
+      const TVector<double> &qloc,
+      const TMatrix<double> &cmat,
       TMatrix<double> &Amat
     ) const override;
 
@@ -269,8 +302,14 @@ class EEQModel : public ChargeModel {
       const TMatrix<double> &dist,
       const TVector<double> &q,
       const TMatrix<double> &Amat,
-      TMatrix<double> &dAmat,
-      TMatrix<double> &atrace
+      TMatrix<double> &dAmatdr,
+      TMatrix<double> &atrace,
+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,
+      const TVector<double> &qloc,
+      const TMatrix<double> &dqlocdr,
+      const TMatrix<double> &cmat,
+      const TMatrix<double> &dcmatdr
     ) const override;
 
     // Calculate the coordination number, forwarding to get_ncoord
@@ -324,8 +363,13 @@ class EEQBCModel : public ChargeModel {
       const int &charge,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
-      TVector<double> &Xvec,
-      TVector<double> &dXvec,
+      const TMatrix<double> &dcndr,
+      TVector<double> &xvec,
+      TMatrix<double> &dxvecdr,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      TMatrix<double> &cmat,
+      TMatrix<double> &dcmatdr,
       bool lgrad
     ) override;
 
@@ -335,6 +379,8 @@ class EEQBCModel : public ChargeModel {
       const TIVector &realIdx,
       const TMatrix<double> &dist,
       const TVector<double> &cn,
+      const TVector<double> &qloc,
+      const TMatrix<double> &cmat,
       TMatrix<double> &Amat
     ) const override;
 
@@ -345,8 +391,14 @@ class EEQBCModel : public ChargeModel {
       const TMatrix<double> &dist,
       const TVector<double> &q,
       const TMatrix<double> &Amat,
-      TMatrix<double> &dAmat,
-      TMatrix<double> &atrace
+      TMatrix<double> &dAmatdr,
+      TMatrix<double> &atrace,
+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,
+      const TVector<double> &qloc,
+      const TMatrix<double> &dqlocdr,
+      const TMatrix<double> &cmat,
+      const TMatrix<double> &dcmatdr
     ) const override;
 
     // Calculate the coordination number, forwarding to get_ncoord
@@ -361,11 +413,13 @@ class EEQBCModel : public ChargeModel {
 
     // Get purely geometry-dependent local charges
     int get_qloc(
-      const TMolecule&,
-      const TIVector &,
-      const TMatrix<double>&,
-      const double,
-      TVector<double>&
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double>& dist,
+      const double q_tot,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      bool lgrad
     );
 
     // Get the capacitance for bond between atoms i and j
@@ -376,31 +430,57 @@ class EEQBCModel : public ChargeModel {
       double &c_ij
     ) const;
 
+    // Get the capacitance for bond between atoms i and j
+    int get_dcpair(
+      double dist_ij,  // distance between atom j to atom i
+      double rvdw_ijat,  // pairwise van der Waals radii for ij atom types
+      double cap_ij,  // product of bond capacitances for atom types of i and j
+      TVector<double> &vec, // Vector from j to i
+      TVector<double> &dcdr_ij  // Out: Capacitance for bond ij
+    ) const;
+
     // Get the capacitance matrix
     int get_cmat(
-      const TMolecule&,
+      const TMolecule& mol,
       const TIVector &realIdx,
-      const TMatrix<double>&,
-      TMatrix<double>&
+      const TMatrix<double>& dist,
+      TMatrix<double>& cmat
+    );
+
+    // Get the capacitance matrix
+    int get_dcmatdr(
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double> &dist,
+      TMatrix<double> &cmat
     );
 
     // Get the right-hand side (electronegativity) of the set of linear equations
     int get_xvec(
-      const TMolecule&,
-      const TIVector&,
-      const TMatrix<double>&,
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double> &dist,
       const TVector<double> &cn,
-      TMatrix<double>&,
-      int,
-      TVector<double>&
+      TMatrix<double> &cmat,
+      int charge,
+      TVector<double> &qloc,
+      TVector<double> &xvec
     );
 
-    // numerical gradient of partial charges w.r.t. atom positions
-    int num_grad_dqdr(
-      TMolecule&,  // molecular geometry
-      const TIVector&,
-      int,
-      TMatrix<double>& // numerical gradient
+    // Get the derivative of the right-hand side of the set of linear equations
+    int get_xvec_derivs(
+      const TMolecule &mol,
+      const TIVector &realIdx,
+      const TMatrix<double> &dist,
+      const TVector<double> &cn,
+      const TMatrix<double> &dcndr,
+      TMatrix<double> &cmat,
+      int charge,
+      TVector<double> &xvec,
+      TMatrix<double> &dxvecdr,
+      TVector<double> &qloc,
+      TMatrix<double> &dqlocdr,
+      TMatrix<double> &dcmatdr
     );
 
 };