Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marvinfriede merged 15 commits into from
Dec 10, 2025
Merged

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

Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
7b77909
The charge derivative dqdr is now implemented for EEQ-BC.
Thomas3R Sep 18, 2025
774a10d
Test charge gradient dqdr against numerical gradient.
Thomas3R Nov 27, 2025
57eddf2
Update src/dftd_eeq.cpp
Thomas3R Dec 1, 2025
7a186a1
Update src/dftd_eeq.cpp
Thomas3R Dec 1, 2025
ef62976
Update src/dftd_eeq.cpp
Thomas3R Dec 1, 2025
c97a660
Update test/test_multi.cpp
Thomas3R Dec 1, 2025
eed31c4
Update test/test_multi.cpp
Thomas3R Dec 1, 2025
58afd5a
Update src/dftd_eeq.cpp
Thomas3R Dec 1, 2025
331c647
Update src/dftd_eeq.cpp
Thomas3R Dec 1, 2025
0125db4
Update test/test_multi.cpp
Thomas3R Dec 1, 2025
ef3968c
Update src/dftd_eeq.cpp
Thomas3R Dec 1, 2025
2563ae6
Update scripts.
Thomas3R Dec 2, 2025
ef31ecb
Update scripts
Thomas3R Dec 2, 2025
2d23006
Address review comments from PR
Thomas3R Dec 5, 2025
df77ae2
Swap dcmatdr dimensions. Use BLAS for Mat x Vec product.
Thomas3R Dec 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 47 additions & 1 deletion include/dftd_cblas.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
C.rows == 0) {
exit(EXIT_FAILURE);
};

// check for transpositions
if (!TransposeA) {
Expand Down Expand Up @@ -236,4 +237,49 @@ inline int BLAS_InvertMatrix(TMatrix<double> &a) {
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
140 changes: 114 additions & 26 deletions include/dftd_eeq.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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;

Expand All @@ -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.
Expand All @@ -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;

Expand All @@ -187,6 +201,12 @@ class ChargeModel
* @param Amat Output coulomb matrix.
* @param dAmat 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.
*/
Expand All @@ -197,7 +217,13 @@ class ChargeModel
const TVector<double> &q,
const TMatrix<double> &Amat,
TMatrix<double> &dAmat,
TMatrix<double> &atrace
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;

/**
Expand Down Expand Up @@ -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;

Expand All @@ -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;

Expand All @@ -270,7 +303,13 @@ class EEQModel : public ChargeModel {
const TVector<double> &q,
const TMatrix<double> &Amat,
TMatrix<double> &dAmat,
TMatrix<double> &atrace
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
Expand Down Expand Up @@ -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;

Expand All @@ -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;

Expand All @@ -346,7 +392,13 @@ class EEQBCModel : public ChargeModel {
const TVector<double> &q,
const TMatrix<double> &Amat,
TMatrix<double> &dAmat,
TMatrix<double> &atrace
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
Expand All @@ -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
Expand All @@ -376,23 +430,57 @@ class EEQBCModel : public ChargeModel {
double &c_ij
) const;

// Get the capacitance for bond between atoms i and j
int get_dcpair(
int iat,
int jat,
TVector<double> &vec,
double &dist_ij,
TVector<double> &dcdr_ij
) const;

// Get the capacitance matrix
int get_cmat(
const TMolecule&,
const TMolecule& mol,
const TIVector &realIdx,
const TMatrix<double>& dist,
TMatrix<double>& cmat
);

// Get the capacitance matrix
int get_dcmatdr(
const TMolecule& mol,
const TIVector &realIdx,
const TMatrix<double>&,
TMatrix<double>&
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
);

// 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
);

// numerical gradient of partial charges w.r.t. atom positions
Expand Down
Loading