From 727f1e78e6aaf36e6a7bd8d82b8de8fd65b2ac9f Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert Date: Sun, 25 Oct 2026 09:59:11 +0100 Subject: [PATCH] Use stable CN weight calculation --- src/dftd3/model.f90 | 158 +++++++++------------------------- test/unit/test_model.f90 | 58 ++++++++++++- test/unit/test_regression.f90 | 17 +--- 3 files changed, 96 insertions(+), 137 deletions(-) diff --git a/src/dftd3/model.f90 b/src/dftd3/model.f90 index 0a2a02a0..6f077f7f 100644 --- a/src/dftd3/model.f90 +++ b/src/dftd3/model.f90 @@ -15,7 +15,6 @@ ! along with s-dftd3. If not, see . 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 @@ -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. !> @@ -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 diff --git a/test/unit/test_model.f90 b/test/unit/test_model.f90 index 2a1d6483..d3739e15 100644 --- a/test/unit/test_model.f90 +++ b/test/unit/test_model.f90 @@ -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 @@ -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 @@ -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 diff --git a/test/unit/test_regression.f90 b/test/unit/test_regression.f90 index 7834c784..49d6b7e1 100644 --- a/test/unit/test_regression.f90 +++ b/test/unit/test_regression.f90 @@ -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 @@ -63,11 +62,6 @@ 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) @@ -75,15 +69,6 @@ subroutine test_dftbplus_871(error) 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