Skip to content

Fix dadr for 0D and 3D #60

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

Draft
lmseidler wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Fix dadr for 0D and 3D #60

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
75d0893
fixed dadr for 0d and 3d
lmseidler Mar 9, 2026
08509cd
fix format
lmseidler Mar 9, 2026
ac01980
vectorization + stripped unused vars
lmseidler Mar 9, 2026
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
36 changes: 24 additions & 12 deletions src/multicharge/model/eeqbc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(:, :)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Comment thread
lmseidler marked this conversation as resolved.
Outdated
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)

Expand Down Expand Up @@ -805,7 +811,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
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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Comment thread
lmseidler marked this conversation as resolved.
Outdated
dadL_local(:, :, iat) = -dgam * qvec(jat) * dgamdL(:, :) + dadL_local(:, :, iat)
dadL_local(:, :, jat) = -dgam * qvec(iat) * dgamdL(:, :) + dadL_local(:, :, jat)

Expand Down
11 changes: 1 addition & 10 deletions test/unit/test_model.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(:)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
12 changes: 1 addition & 11 deletions test/unit/test_pbc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(:)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down