From 75d089384e7ae1d3728c688dc5edc7c6aa26ff33 Mon Sep 17 00:00:00 2001 From: Leopold Seidler Date: Mon, 9 Mar 2026 16:56:41 +0100 Subject: [PATCH 1/3] fixed dadr for 0d and 3d --- src/multicharge/model/eeqbc.f90 | 44 +++++++++++++++++++++------------ test/unit/test_model.f90 | 11 +-------- test/unit/test_pbc.f90 | 12 +-------- 3 files changed, 30 insertions(+), 37 deletions(-) diff --git a/src/multicharge/model/eeqbc.f90 b/src/multicharge/model/eeqbc.f90 index bd1db2da..251ec11d 100644 --- a/src/multicharge/model/eeqbc.f90 +++ b/src/multicharge/model/eeqbc.f90 @@ -671,7 +671,7 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & real(wp), intent(out) :: dadL(:, :, :) real(wp), intent(out) :: atrace(:, :) - integer :: iat, jat, izp, jzp + integer :: iat, jat, izp, jzp, kat real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3), dgamdL(3, 3) real(wp), allocatable :: dgamdr(:, :) @@ -689,7 +689,7 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & !$omp parallel default(none) & !$omp shared(atrace, dadr, dadL, mol, self, cn, qloc, qvec) & !$omp shared(cmat, dcdr, dcdL, dcndr, dcndL, dqlocdr, dqlocdL) & - !$omp private(iat, izp, jat, jzp, gam, vec, r2, dtmp, norm_cn, arg) & + !$omp private(iat, izp, jat, jzp, kat, gam, vec, r2, dtmp, norm_cn, arg) & !$omp private(radi, radj, dradi, dradj, dgamdr, dgamdL, dG, dS) & !$omp private(atrace_local, dadr_local, dadL_local) allocate(atrace_local, source=atrace) @@ -733,10 +733,16 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & ! Effective charge width derivative dtmp = 2.0_wp * exp(-arg) / (sqrtpi) - atrace_local(:, iat) = -dtmp * qvec(jat) * dgamdr(:, jat) * cmat(jat, iat) + atrace_local(:, iat) - atrace_local(:, jat) = -dtmp * qvec(iat) * dgamdr(:, iat) * cmat(iat, jat) + atrace_local(:, jat) - dadr_local(:, iat, jat) = +dtmp * qvec(iat) * dgamdr(:, iat) * cmat(iat, jat) + dadr_local(:, iat, jat) - dadr_local(:, jat, iat) = +dtmp * qvec(jat) * dgamdr(:, jat) * cmat(jat, iat) + dadr_local(:, jat, iat) + atrace_local(:, iat) = +dtmp * qvec(jat) * dgamdr(:, iat) * cmat(jat, iat) + atrace_local(:, iat) + atrace_local(:, jat) = +dtmp * qvec(iat) * dgamdr(:, jat) * cmat(iat, jat) + atrace_local(:, jat) + do kat = 1, mol%nat + if (kat /= iat) then + dadr_local(:, kat, iat) = +dtmp * qvec(jat) * dgamdr(:, kat) * cmat(jat, iat) + dadr_local(:, kat, iat) + end if + if (kat /= jat) then + dadr_local(:, kat, jat) = +dtmp * qvec(iat) * dgamdr(:, kat) * cmat(iat, jat) + dadr_local(:, kat, jat) + end if + end do dadL_local(:, :, iat) = +dtmp * qvec(jat) * dgamdL(:, :) * cmat(jat, iat) + dadL_local(:, :, iat) dadL_local(:, :, jat) = +dtmp * qvec(iat) * dgamdL(:, :) * cmat(iat, jat) + dadL_local(:, :, jat) @@ -805,11 +811,11 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & real(wp), intent(out) :: dadL(:, :, :) real(wp), intent(out) :: atrace(:, :) - integer :: iat, jat, izp, jzp, img - real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn, rvdw, wsw, dgam - real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3) - real(wp) :: dgamdL(3, 3), capi, capj - real(wp), allocatable :: dgamdr(:, :), dtrans(:, :) + integer :: iat, jat, izp, jzp, img, kat + real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn, rvdw, wsw, dgam + real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3) + real(wp) :: dgamdL(3, 3), capi, capj + real(wp), allocatable :: dgamdr(:, :), dtrans(:, :) ! Thread-private arrays for reduction real(wp), allocatable :: atrace_local(:, :) @@ -826,7 +832,7 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & !$omp parallel default(none) & !$omp shared(self, mol, cn, qloc, qvec, wsc, dadr, dadL, atrace) & !$omp shared (cmat, dcdr, dcdL, dcndr, dcndL, dqlocdr, dqlocdL, dtrans) & - !$omp private(iat, izp, jat, jzp, img, gam, vec, r2, dtmp, norm_cn, arg, rvdw) & + !$omp private(iat, izp, jat, jzp, kat, img, gam, vec, r2, dtmp, norm_cn, arg, rvdw) & !$omp private(radi, radj, dradi, dradj, capi, capj, dgamdr, dgamdL, dG, dS, wsw) & !$omp private(dgam, dadr_local, dadL_local, atrace_local) allocate(atrace_local, source=atrace) @@ -875,10 +881,16 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & dadL_local(:, :, iat) = +dS * qvec(jat) + dadL_local(:, :, iat) ! Effective charge width derivative - atrace_local(:, iat) = +dgam * qvec(jat) * dgamdr(:, jat) + atrace_local(:, iat) - atrace_local(:, jat) = +dgam * qvec(iat) * dgamdr(:, iat) + atrace_local(:, jat) - dadr_local(:, iat, jat) = -dgam * qvec(iat) * dgamdr(:, iat) + dadr_local(:, iat, jat) - dadr_local(:, jat, iat) = -dgam * qvec(jat) * dgamdr(:, jat) + dadr_local(:, jat, iat) + atrace_local(:, iat) = -dgam * qvec(jat) * dgamdr(:, iat) + atrace_local(:, iat) + atrace_local(:, jat) = -dgam * qvec(iat) * dgamdr(:, jat) + atrace_local(:, jat) + do kat = 1, mol%nat + if (kat /= iat) then + dadr_local(:, kat, iat) = -dgam * qvec(jat) * dgamdr(:, kat) + dadr_local(:, kat, iat) + end if + if (kat /= jat) then + dadr_local(:, kat, jat) = -dgam * qvec(iat) * dgamdr(:, kat) + dadr_local(:, kat, jat) + end if + end do dadL_local(:, :, iat) = -dgam * qvec(jat) * dgamdL(:, :) + dadL_local(:, :, iat) dadL_local(:, :, jat) = -dgam * qvec(iat) * dgamdL(:, :) + dadL_local(:, :, jat) diff --git a/test/unit/test_model.f90 b/test/unit/test_model.f90 index 1e109191..529f3aa5 100644 --- a/test/unit/test_model.f90 +++ b/test/unit/test_model.f90 @@ -98,7 +98,6 @@ subroutine test_dadr(error, mol, model) type(error_type), allocatable, intent(out) :: error integer :: iat, ic, jat, kat - real(wp) :: thr2_local real(wp), parameter :: trans(3, 1) = 0.0_wp real(wp), parameter :: step = 1.0e-6_wp real(wp), allocatable :: cn(:) @@ -115,14 +114,6 @@ subroutine test_dadr(error, mol, model) & dqlocdL(3, 3, mol%nat), dadr(3, mol%nat, mol%nat + 1), dadL(3, 3, mol%nat + 1), & & atrace(3, mol%nat), numtrace(3, mol%nat), numgrad(3, mol%nat, mol%nat + 1), qvec(mol%nat)) - ! Set tolerance higher if testing eeqbc model - select type (model) - type is (eeqbc_model) - thr2_local = 3.0_wp*thr2 - class default - thr2_local = thr2 - end select - ! Obtain the vector of charges call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) @@ -195,7 +186,7 @@ subroutine test_dadr(error, mol, model) dadr(:, iat, iat) = atrace(:, iat) + dadr(:, iat, iat) end do - if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2_local)) then + if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2)) then call test_failed(error, "Derivative of the A matrix does not match") print'(a)', "dadr:" print'(3es21.12)', dadr diff --git a/test/unit/test_pbc.f90 b/test/unit/test_pbc.f90 index 07118ae7..1d4c2821 100644 --- a/test/unit/test_pbc.f90 +++ b/test/unit/test_pbc.f90 @@ -467,7 +467,6 @@ subroutine test_dadr(error, mol, model) type(error_type), allocatable, intent(out) :: error integer :: iat, ic, jat, kat - real(wp) :: thr2_local real(wp), parameter :: cutoff = 25.0_wp real(wp), parameter :: step = 1.0e-6_wp real(wp), allocatable :: cn(:) @@ -484,14 +483,6 @@ subroutine test_dadr(error, mol, model) & dqlocdL(3, 3, mol%nat), dadr(3, mol%nat, mol%nat + 1), dadL(3, 3, mol%nat + 1), & & atrace(3, mol%nat), numtrace(3, mol%nat), numgrad(3, mol%nat, mol%nat + 1), qvec(mol%nat)) - ! Set tolerance higher if testing eeqbc model - select type(model) - type is(eeqbc_model) - thr2_local = 3.0_wp * thr2 - class default - thr2_local = thr2 - end select - call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) ! Obtain the vector of charges @@ -543,8 +534,7 @@ subroutine test_dadr(error, mol, model) dadr(:, iat, iat) = atrace(:, iat) + dadr(:, iat, iat) end do - ! higher tolerance for numerical gradient - if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2_local)) then + if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2)) then call test_failed(error, "Derivative of the A matrix does not match") print'(a)', "dadr:" print'(3es21.14)', dadr From 08509cd55dd0415dab179d4ff029a4587d1571c1 Mon Sep 17 00:00:00 2001 From: Leopold Seidler Date: Mon, 9 Mar 2026 16:58:27 +0100 Subject: [PATCH 2/3] fix format --- src/multicharge/model/eeqbc.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/multicharge/model/eeqbc.f90 b/src/multicharge/model/eeqbc.f90 index 251ec11d..9644dd9d 100644 --- a/src/multicharge/model/eeqbc.f90 +++ b/src/multicharge/model/eeqbc.f90 @@ -811,11 +811,11 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & real(wp), intent(out) :: dadL(:, :, :) real(wp), intent(out) :: atrace(:, :) - integer :: iat, jat, izp, jzp, img, kat - real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn, rvdw, wsw, dgam - real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3) - real(wp) :: dgamdL(3, 3), capi, capj - real(wp), allocatable :: dgamdr(:, :), dtrans(:, :) + integer :: iat, jat, izp, jzp, img, kat + real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn, rvdw, wsw, dgam + real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3) + real(wp) :: dgamdL(3, 3), capi, capj + real(wp), allocatable :: dgamdr(:, :), dtrans(:, :) ! Thread-private arrays for reduction real(wp), allocatable :: atrace_local(:, :) From ac0198007f4e4bb4b13be85fc435c85488e64f76 Mon Sep 17 00:00:00 2001 From: Leopold Seidler Date: Mon, 9 Mar 2026 17:23:13 +0100 Subject: [PATCH 3/3] vectorization + stripped unused vars --- src/multicharge/model/eeqbc.f90 | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/src/multicharge/model/eeqbc.f90 b/src/multicharge/model/eeqbc.f90 index 9644dd9d..5930baf9 100644 --- a/src/multicharge/model/eeqbc.f90 +++ b/src/multicharge/model/eeqbc.f90 @@ -735,14 +735,10 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & dtmp = 2.0_wp * exp(-arg) / (sqrtpi) atrace_local(:, iat) = +dtmp * qvec(jat) * dgamdr(:, iat) * cmat(jat, iat) + atrace_local(:, iat) atrace_local(:, jat) = +dtmp * qvec(iat) * dgamdr(:, jat) * cmat(iat, jat) + atrace_local(:, jat) - do kat = 1, mol%nat - if (kat /= iat) then - dadr_local(:, kat, iat) = +dtmp * qvec(jat) * dgamdr(:, kat) * cmat(jat, iat) + dadr_local(:, kat, iat) - end if - if (kat /= jat) then - dadr_local(:, kat, jat) = +dtmp * qvec(iat) * dgamdr(:, kat) * cmat(iat, jat) + dadr_local(:, kat, jat) - end if - end do + dadr_local(:, :, iat) = +dtmp * qvec(jat) * dgamdr(:, :) * cmat(jat, iat) + dadr_local(:, :, iat) + dadr_local(:, iat, iat) = -dtmp * qvec(jat) * dgamdr(:, iat) * cmat(jat, iat) + dadr_local(:, iat, iat) + dadr_local(:, :, jat) = +dtmp * qvec(iat) * dgamdr(:, :) * cmat(iat, jat) + dadr_local(:, :, jat) + dadr_local(:, jat, jat) = -dtmp * qvec(iat) * dgamdr(:, jat) * cmat(iat, jat) + dadr_local(:, jat, jat) dadL_local(:, :, iat) = +dtmp * qvec(jat) * dgamdL(:, :) * cmat(jat, iat) + dadL_local(:, :, iat) dadL_local(:, :, jat) = +dtmp * qvec(iat) * dgamdL(:, :) * cmat(iat, jat) + dadL_local(:, :, jat) @@ -811,7 +807,7 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & real(wp), intent(out) :: dadL(:, :, :) real(wp), intent(out) :: atrace(:, :) - integer :: iat, jat, izp, jzp, img, kat + integer :: iat, jat, izp, jzp, img real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn, rvdw, wsw, dgam real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3) real(wp) :: dgamdL(3, 3), capi, capj @@ -832,7 +828,7 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & !$omp parallel default(none) & !$omp shared(self, mol, cn, qloc, qvec, wsc, dadr, dadL, atrace) & !$omp shared (cmat, dcdr, dcdL, dcndr, dcndL, dqlocdr, dqlocdL, dtrans) & - !$omp private(iat, izp, jat, jzp, kat, img, gam, vec, r2, dtmp, norm_cn, arg, rvdw) & + !$omp private(iat, izp, jat, jzp, img, gam, vec, r2, dtmp, norm_cn, arg, rvdw) & !$omp private(radi, radj, dradi, dradj, capi, capj, dgamdr, dgamdL, dG, dS, wsw) & !$omp private(dgam, dadr_local, dadL_local, atrace_local) allocate(atrace_local, source=atrace) @@ -883,14 +879,10 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, & ! Effective charge width derivative atrace_local(:, iat) = -dgam * qvec(jat) * dgamdr(:, iat) + atrace_local(:, iat) atrace_local(:, jat) = -dgam * qvec(iat) * dgamdr(:, jat) + atrace_local(:, jat) - do kat = 1, mol%nat - if (kat /= iat) then - dadr_local(:, kat, iat) = -dgam * qvec(jat) * dgamdr(:, kat) + dadr_local(:, kat, iat) - end if - if (kat /= jat) then - dadr_local(:, kat, jat) = -dgam * qvec(iat) * dgamdr(:, kat) + dadr_local(:, kat, jat) - end if - end do + dadr_local(:, :, iat) = -dgam * qvec(jat) * dgamdr(:, :) + dadr_local(:, :, iat) + dadr_local(:, iat, iat) = +dgam * qvec(jat) * dgamdr(:, iat) + dadr_local(:, iat, iat) + dadr_local(:, :, jat) = -dgam * qvec(iat) * dgamdr(:, :) + dadr_local(:, :, jat) + dadr_local(:, jat, jat) = +dgam * qvec(iat) * dgamdr(:, jat) + dadr_local(:, jat, jat) dadL_local(:, :, iat) = -dgam * qvec(jat) * dgamdL(:, :) + dadL_local(:, :, iat) dadL_local(:, :, jat) = -dgam * qvec(iat) * dgamdL(:, :) + dadL_local(:, :, jat) @@ -1330,7 +1322,7 @@ subroutine get_dcpair_dir(kbc, rij, trans, rvdw, capi, capj, dgpair, dspair) real(wp), intent(out) :: dspair(3, 3) integer :: itr - real(wp) :: r1, arg, dtmp, dgtmp(3), dstmp(3, 3), vec(3) + real(wp) :: r1, dgtmp(3), dstmp(3, 3), vec(3) dgpair(:) = 0.0_wp dspair(:, :) = 0.0_wp