Skip to content

Use stable CN weight calculation #246

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
awvwgk merged 1 commit into from
Oct 25, 2026
Merged

Use stable CN weight calculation #246

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
158 changes: 39 additions & 119 deletions src/dftd3/model.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@
! along with s-dftd3. If not, see <https://www.gnu.org/licenses/>.

module dftd3_model
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan
use dftd3_citation, only : citation_type, get_citation, doi_fourier_d3
use dftd3_data, only : get_r4r2_val, get_vdw_rad
use dftd3_fourier_decomposition, only : d3_lowrank_c6, d3_lowrank_config, new_lowrank_c6
Expand Down Expand Up @@ -218,120 +217,49 @@ subroutine weight_references(self, mol, cn, gwvec, gwdcn, gwd2cn)
!> second derivative of the weighting function w.r.t. the coordination number
real(wp), intent(out), optional :: gwd2cn(:, :)

integer :: iat, izp, iref, mref
real(wp) :: norm, dnorm, d2norm, gw, expw, expd, expdd, gwk, dgwk, d2gwk
real(wp) :: wf2, dcn
real(wp), allocatable :: gwt(:)

mref = maxval(self%ref)

if (present(gwdcn)) then
gwvec(:, :) = 0.0_wp
gwdcn(:, :) = 0.0_wp
if (present(gwd2cn)) gwd2cn(:, :) = 0.0_wp
wf2 = 2 * self%wf

!$omp parallel default(none) &
!$omp shared(gwvec, gwdcn, gwd2cn, mol, self, cn, wf2, mref) &
!$omp private(iat, izp, iref, norm, dnorm, d2norm, gw, expw, expd, expdd, &
!$omp& gwk, dgwk, d2gwk, dcn, gwt)
allocate(gwt(mref))
!$omp do schedule(runtime)
do iat = 1, mol%nat
izp = mol%id(iat)
norm = 0.0_wp
dnorm = 0.0_wp
d2norm = 0.0_wp
do iref = 1, self%ref(izp)
gw = weight_cn(self%wf, cn(iat), self%cn(iref, izp))
gwt(iref) = gw
dcn = self%cn(iref, izp) - cn(iat)
norm = norm + gw
dnorm = dnorm + wf2 * dcn * gw
d2norm = d2norm + (wf2*wf2*dcn*dcn - wf2) * gw
end do
norm = 1.0_wp / norm
do iref = 1, self%ref(izp)
expw = gwt(iref)
dcn = self%cn(iref, izp) - cn(iat)
expd = wf2 * dcn * expw
expdd = (wf2*wf2*dcn*dcn - wf2) * expw
gwk = expw * norm
if (is_exceptional(gwk)) then
if (maxval(self%cn(:self%ref(izp), izp)) == self%cn(iref, izp)) then
gwk = 1.0_wp
else
gwk = 0.0_wp
end if
end if
gwvec(iref, iat) = gwk

dgwk = expd * norm - expw * dnorm * norm**2
if (is_exceptional(dgwk)) then
dgwk = 0.0_wp
end if
gwdcn(iref, iat) = dgwk

if (present(gwd2cn)) then
d2gwk = expdd * norm - 2.0_wp * expd * dnorm * norm**2 &
& - expw * d2norm * norm**2 + 2.0_wp * expw * dnorm**2 * norm**3
if (is_exceptional(d2gwk)) then
d2gwk = 0.0_wp
end if
gwd2cn(iref, iat) = d2gwk
end if
end do
end do
!$omp end do
deallocate(gwt)
!$omp end parallel

else

gwvec(:, :) = 0.0_wp

!$omp parallel default(none) &
!$omp shared(gwvec, mol, self, cn, mref) &
!$omp private(iat, izp, iref, norm, gw, expw, gwk, gwt)
allocate(gwt(mref))
!$omp do schedule(runtime)
do iat = 1, mol%nat
izp = mol%id(iat)
norm = 0.0_wp
do iref = 1, self%ref(izp)
gw = weight_cn(self%wf, cn(iat), self%cn(iref, izp))
gwt(iref) = gw
norm = norm + gw
end do
norm = 1.0_wp / norm
do iref = 1, self%ref(izp)
gwk = gwt(iref) * norm
if (is_exceptional(gwk)) then
if (maxval(self%cn(:self%ref(izp), izp)) == self%cn(iref, izp)) then
gwk = 1.0_wp
else
gwk = 0.0_wp
end if
end if
gwvec(iref, iat) = gwk
end do
end do
!$omp end do
deallocate(gwt)
!$omp end parallel
end if
integer :: iat, izp, nref
real(wp) :: dmin, norm, mean, var
real(wp), allocatable :: dcn(:), gw(:), slope(:)

gwvec(:, :) = 0.0_wp
if (present(gwdcn)) gwdcn(:, :) = 0.0_wp
if (present(gwd2cn)) gwd2cn(:, :) = 0.0_wp

!$omp parallel default(none) &
!$omp shared(gwvec, gwdcn, gwd2cn, mol, self, cn) &
!$omp private(iat, izp, nref, dmin, norm, mean, var, dcn, gw, slope)
allocate(dcn(maxval(self%ref)), gw(maxval(self%ref)), slope(maxval(self%ref)))
!$omp do schedule(runtime)
do iat = 1, mol%nat
izp = mol%id(iat)
nref = self%ref(izp)
dcn(:nref) = self%cn(:nref, izp) - cn(iat)
! Shift by the closest reference to keep the largest exponent at zero
! and the normalization at least one (log-sum-exp), as done in D3S
dmin = minval(dcn(:nref)**2)
gw(:nref) = exp(-self%wf * (dcn(:nref)**2 - dmin))
norm = sum(gw(:nref))
gw(:nref) = gw(:nref) / norm
gwvec(:nref, iat) = gw(:nref)

if (present(gwdcn)) then
slope(:nref) = 2 * self%wf * dcn(:nref)
mean = sum(gw(:nref) * slope(:nref))
slope(:nref) = slope(:nref) - mean
gwdcn(:nref, iat) = gw(:nref) * slope(:nref)
if (present(gwd2cn)) then
var = sum(gw(:nref) * slope(:nref)**2)
gwd2cn(:nref, iat) = gw(:nref) * (slope(:nref)**2 - var)
end if
end if
end do
!$omp end do
deallocate(dcn, gw, slope)
!$omp end parallel

end subroutine weight_references


!> Check whether we are dealing with an exceptional value, NaN or Inf
elemental function is_exceptional(val)
real(wp), intent(in) :: val
logical :: is_exceptional
is_exceptional = ieee_is_nan(val) .or. abs(val) > huge(val)
end function is_exceptional


!> Calculate atomic dispersion coefficients and their derivatives w.r.t.
!> the coordination number.
!>
Expand Down Expand Up @@ -491,12 +419,4 @@ subroutine get_atomic_c6(self, mol, gwvec, gwdcn, c6, dc6dcn, gwd2cn, d2c6dcn2,
end subroutine get_atomic_c6


elemental function weight_cn(wf,cn,cnref) result(cngw)
real(wp),intent(in) :: wf, cn, cnref
real(wp) :: cngw
intrinsic :: exp
cngw = exp ( -wf * ( cn - cnref )**2 )
end function weight_cn


end module dftd3_model
58 changes: 56 additions & 2 deletions test/unit/test_model.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ module test_model
use mctc_env, only : wp
use mctc_env_testing, only : new_unittest, unittest_type, error_type, check, &
& test_failed
use mctc_io_structure, only : structure_type
use mctc_io_structure, only : structure_type, new
use mctc_ncoord, only : ncoord_type, new_ncoord, cn_count
use mstore, only : get_structure
implicit none
Expand Down Expand Up @@ -53,7 +53,8 @@ subroutine collect_model(testsuite)
& new_unittest("gw-amf3", test_gw_amf3), &
& new_unittest("ghost-atoms", test_ghost_atoms), &
& new_unittest("dgw-mb04", test_dgw_mb04), &
& new_unittest("dgw-mb05", test_dgw_mb05) &
& new_unittest("dgw-mb05", test_dgw_mb05), &
& new_unittest("gw-stable", test_gw_stable) &
& ]

end subroutine collect_model
Expand Down Expand Up @@ -336,6 +337,59 @@ subroutine test_dgw_mb05(error)
end subroutine test_dgw_mb05


!> Weights and derivatives must stay finite and normalized for CNs far away
!> from all references, where the unshifted Gaussians underflow
subroutine test_gw_stable(error)

!> Error handling
type(error_type), allocatable, intent(out) :: error

type(structure_type) :: mol
type(d3_model) :: d3
integer :: ic, iat, nref
real(wp), parameter :: step = 1.0e-5_wp
real(wp), parameter :: cns(*) = [-100.0_wp, -13.5_wp, 0.0_wp, 2.5_wp, 14.0_wp, 17.4_wp, 17.5_wp, &
& 30.0_wp, 100.0_wp, 1.0e4_wp]
real(wp) :: cn(3), xyz(3, 3)
real(wp), allocatable :: gwvec(:, :), gwdcn(:, :), gwd2cn(:, :)
real(wp), allocatable :: gwr(:, :), gwl(:, :), dgwr(:, :), dgwl(:, :)

xyz(:, :) = 0.0_wp
call new(mol, [1, 6, 73], xyz)
call new_d3_model(d3, mol)
allocate(gwvec(maxval(d3%ref), mol%nat), source=0.0_wp)
allocate(gwdcn, gwd2cn, gwr, gwl, dgwr, dgwl, mold=gwvec)

do ic = 1, size(cns)
cn(:) = cns(ic)
call d3%weight_references(mol, cn, gwvec, gwdcn, gwd2cn)
call d3%weight_references(mol, cn + step, gwr, dgwr)
call d3%weight_references(mol, cn - step, gwl, dgwl)

do iat = 1, mol%nat
nref = d3%ref(mol%id(iat))
call check(error, sum(gwvec(:nref, iat)), 1.0_wp, thr=thr)
if (allocated(error)) return
call check(error, minval(gwvec(:nref, iat)) >= 0.0_wp)
if (allocated(error)) return
! far away from all references only the closest one survives
if (abs(cns(ic)) > 50.0_wp) then
call check(error, gwvec(minloc(abs(d3%cn(:nref, mol%id(iat)) - cns(ic)), 1), iat), &
& 1.0_wp, thr=thr)
if (allocated(error)) return
end if
end do

call check(error, norm2(gwdcn - (gwr - gwl) / (2*step)), 0.0_wp, thr=thr2)
if (allocated(error)) return
call check(error, norm2(gwd2cn - (dgwr - dgwl) / (2*step)), 0.0_wp, thr=10*thr2)
if (allocated(error)) return
end do

end subroutine test_gw_stable



subroutine test_ghost_atoms(error)

!> Error handling
Expand Down
17 changes: 1 addition & 16 deletions test/unit/test_regression.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,7 @@ module test_regression
use dftd3, only : get_dispersion, realspace_cutoff, rational_damping_param, &
& d3_model, new_d3_model
use mctc_env, only : wp
use mctc_env_testing, only : new_unittest, unittest_type, error_type, check, &
& test_failed
use mctc_env_testing, only : new_unittest, unittest_type, error_type, check
use mctc_io, only : structure_type, new
implicit none
private
Expand Down Expand Up @@ -63,27 +62,13 @@ subroutine test_dftbplus_871(error)
call new_d3_model(d3, mol)
call get_dispersion(mol, d3, param, realspace_cutoff(), energy, gradient, sigma)

if (is_exceptional(norm2(gradient))) then
call test_failed(error, "Exceptional value propagated to the gradient")
return
end if

call check(error, energy, -8.9964956516666404E-5_wp, thr=thr)
if (allocated(error)) return
call check(error, norm2(gradient), 8.9844177273972976E-9_wp, thr=thr)

end subroutine test_dftbplus_871


!> Check whether we are dealing with an exceptional value, NaN or Inf
elemental function is_exceptional(val)
use, intrinsic :: ieee_arithmetic, only : ieee_is_nan
real(wp), intent(in) :: val
logical :: is_exceptional
is_exceptional = ieee_is_nan(val) .or. abs(val) > huge(val)
end function is_exceptional


subroutine structure_dftbplus_871(self)
type(structure_type), intent(out) :: self
integer, parameter :: nat = 90
Expand Down
Loading