Minor cleanup (#18)

- Remove unused LAPACK/BLAS functions
- add some more documentation
- format

.clang-format

@@ -4,6 +4,7 @@ BasedOnStyle: LLVM
 AccessModifierOffset: 2
 ConstructorInitializerIndentWidth: 2
 ContinuationIndentWidth: 2
+IndentAccessModifiers: true
 IndentWidth: 2
 IndentWrappedFunctionNames: true
 TabWidth: 2
@@ -23,3 +24,4 @@ AllowShortIfStatementsOnASingleLine: AllIfsAndElse
 
 AlignOperands: true
 AllowAllArgumentsOnNextLine: true
+InsertTrailingCommas: Wrapped

include/dftd_cblas.h

@@ -24,6 +24,16 @@
 
 namespace dftd4 {
 
+/**
+ * @brief General matrix vector multiplication (`C = alpha * A * V + C`).
+ * 
+ * @param C Result vector C. Modified in-place.
+ * @param A Matrix A.
+ * @param V Vector V.
+ * @param Transpose Specifies whether to transpose matrix A.
+ * @param alpha Scaling factor for the product of matrix A and vector X.
+ * @return Exit code
+ */
 inline int BLAS_Add_Mat_x_Vec(
   TVector<double> &C,
   TMatrix<double> &A,
@@ -70,34 +80,18 @@ inline int BLAS_Add_Mat_x_Vec(
   return EXIT_FAILURE;
 };
 
-inline void BLAS_Add_Mat_x_Vec(
-  double *C,
-  const double *A,
-  const double *B,
-  const int rows,
-  const int cols,
-  const bool Transpose,
-  const double alpha,
-  const double beta
-) {
-  if (Transpose) {
-    cblas_dgemv(
-      CblasRowMajor, CblasTrans, rows, cols, alpha, A, cols, B, 1, beta, C, 1
-    );
-  } else {
-    cblas_dgemv(
-      CblasRowMajor, CblasNoTrans, rows, cols, alpha, A, cols, B, 1, beta, C, 1
-    );
-  };
-};
-
-inline double
-  BLAS_Vec_x_Vec(const TVector<double> &A, const TVector<double> &B) {
-  if (A.N != B.N) exit(EXIT_FAILURE);
-  return cblas_ddot(A.N, A.p, 1, B.p, 1);
-};
 
-/// C = alpha * A * B + C
+/**
+ * @brief General matrix-matrix multiplication (`C = alpha * A * B + C`).
+ * 
+ * @param C Result matrix C. Modified in-place.
+ * @param A Matrix A.
+ * @param B Matrix B.
+ * @param TransposeA Specifies whether to transpose matrix A.
+ * @param TransposeB Specifies whether to transpose matrix B.
+ * @param alpha Scaling factor for the product of matrix A and matrix B.
+ * @return Exit code.
+ */
 inline int BLAS_Add_Mat_x_Mat(
   TMatrix<double> &C,
   const TMatrix<double> &A,
@@ -209,37 +203,19 @@ inline int BLAS_Add_Mat_x_Mat(
   return EXIT_SUCCESS;
 };
 
-// Linear equation solver
-inline int BLAS_LINEQ(TMatrix<double> &a, TMatrix<double> &b, int m) {
-  lapack_int info, n, nrhs;
-  lapack_int *ipiv;
-
-  if (a.rows != a.cols) return EXIT_FAILURE;
-
-  n = a.rows;
-  nrhs = m;
-
-  a.Transpose();
-  b.Transpose();
-
-  ipiv = new lapack_int[n];
-
-  info = LAPACKE_dgesv(LAPACK_COL_MAJOR, n, nrhs, a.p, n, ipiv, b.p, n);
-
-  delete[] ipiv;
-
-  a.Transpose();
-  b.Transpose();
-
-  return info;
-}
-
+/**
+ * @brief Compute inverse of a matrix using LU decomposition.
+ * 
+ * @param a Matrix a.
+ * @return Exit code.
+ */
 inline int BLAS_InvertMatrix(TMatrix<double> &a) {
   if (a.rows != a.cols) { return EXIT_FAILURE; }
 
   lapack_int info;
   lapack_int *ipiv = new lapack_int[a.rows];
 
+  // LU factorization of a general m-by-n matrix
   info = LAPACKE_dgetrf(
     LAPACK_ROW_MAJOR,
     (lapack_int)a.rows,
@@ -248,13 +224,12 @@ inline int BLAS_InvertMatrix(TMatrix<double> &a) {
     (lapack_int)a.cols,
     ipiv
   );
-
   if (info != 0) { return EXIT_FAILURE; }
 
+  // Inverse of an LU-factored general matrix
   info = LAPACKE_dgetri(
     LAPACK_ROW_MAJOR, (lapack_int)a.rows, a.p, (lapack_int)a.cols, ipiv
   );
-
   if (info != 0) { return EXIT_FAILURE; }
 
   delete[] ipiv;

include/dftd_dispersion.h

@@ -28,14 +28,14 @@
 namespace dftd4 {
 
 class dparam {
-    public:
-  double s6;
-  double s8;
-  double s10;
-  double s9;
-  double a1;
-  double a2;
-  int alp;
+  public:
+    double s6;
+    double s8;
+    double s10;
+    double s9;
+    double a1;
+    double a2;
+    int alp;
 };
 
 /**
@@ -46,16 +46,16 @@ class dparam {
  * @param par DFT-D4 parameters.
  * @param d4 Base D4 dispersion model.
  * @param cutoff Real-space cutoffs for CN and dispersion.
- * @param energy Dispersion energy.
- * @param GRAD Dispersion gradient.
+ * @param energy Dispersion energy (inout).
+ * @param GRAD Dispersion gradient (inout).
  * @return Exit status.
  */
 extern int get_dispersion(
   const TMolecule &mol,
   const int charge,
-  TD4Model &d4,
+  const TD4Model &d4,
   const dparam &par,
-  TCutoff cutoff,
+  const TCutoff cutoff,
   double &energy,
   double *GRAD
 );

include/dftd_model.h

@@ -23,49 +23,49 @@
 
 namespace dftd4 {
 
-// Default weighting factor for coordination number interpolation
-static const double wf_default = 6.0;
-
 // Default maximum charge scaling height for partial charge extrapolation
 static const double ga_default = 3.0;
 
 // Default charge scaling steepness for partial charge extrapolation
 static const double gc_default = 2.0;
 
+// Default weighting factor for coordination number interpolation
+static const double wf_default = 6.0;
+
 class TD4Model {
-    public:
-  double wf;
-  double ga;
-  double gc;
-
-  explicit TD4Model(
-    double wf_scale = wf_default,
-    double ga_scale = ga_default,
-    double gc_scale = gc_default
-  );
-
-  int weight_references(
-    const TMolecule &mol,
-    const TVector<double> &cn,
-    const TVector<double> &q,
-    TMatrix<double> &gwvec,
-    TMatrix<double> &dgwdcn,
-    TMatrix<double> &dgwdq,
-    bool lgrad = false
-  );
-
-  int get_atomic_c6(
-    const TMolecule &mol,
-    const TMatrix<double> &gwvec,
-    const TMatrix<double> &dgwdcn,
-    const TMatrix<double> &dgwdq,
-    TMatrix<double> &c6,
-    TMatrix<double> &dc6dcn,
-    TMatrix<double> &dc6dq,
-    bool lgrad = false
-  );
-
-  int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha);
+  public:
+    double ga;
+    double gc;
+    double wf;
+
+    explicit TD4Model(
+      double ga_scale = ga_default,
+      double gc_scale = gc_default,
+      double wf_scale = wf_default
+    );
+
+    int weight_references(
+      const TMolecule &mol,
+      const TVector<double> &cn,
+      const TVector<double> &q,
+      TMatrix<double> &gwvec,
+      TMatrix<double> &dgwdcn,
+      TMatrix<double> &dgwdq,
+      bool lgrad = false
+    ) const;
+
+    int get_atomic_c6(
+      const TMolecule &mol,
+      const TMatrix<double> &gwvec,
+      const TMatrix<double> &dgwdcn,
+      const TMatrix<double> &dgwdq,
+      TMatrix<double> &c6,
+      TMatrix<double> &dc6dcn,
+      TMatrix<double> &dc6dq,
+      bool lgrad = false
+    ) const;
+
+    int set_refalpha_eeq(const TMolecule &mol, TMatrix<double> &alpha) const;
 };
 
 extern inline double trapzd(const double a[23], const double b[23]);