From fa01f374131dd21869c6d2c98a02c5d467bc3596 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Thu, 24 Jul 2025 18:09:30 +0200 Subject: [PATCH] Bug fix for initial guess with ECP potentials (#4343) --- src/atom_kind_orbitals.F | 2 +- src/atom_operators.F | 17 +++++++++++++++-- tests/QS/regtest-ecp-2/TEST_FILES.toml | 4 ++-- tests/QS/regtest-ecp/TEST_FILES.toml | 4 ++-- 4 files changed, 20 insertions(+), 7 deletions(-) diff --git a/src/atom_kind_orbitals.F b/src/atom_kind_orbitals.F index f5abf42859..9231df0a3e 100644 --- a/src/atom_kind_orbitals.F +++ b/src/atom_kind_orbitals.F @@ -1282,7 +1282,7 @@ CONTAINS elec_conf=ppeconf) ecp_atompot%zion = zeff ecp_atompot%econf(0:3) = ppeconf(0:3) - CALL get_potential(sgp_potential, lmax=ecp_atompot%lmax) + CALL get_potential(sgp_potential, sl_lmax=ecp_atompot%lmax) ! local CALL get_potential(sgp_potential, nloc=ecp_atompot%nloc, nrloc=ecp_atompot%nrloc, & aloc=ecp_atompot%aloc, bloc=ecp_atompot%bloc) diff --git a/src/atom_operators.F b/src/atom_operators.F index 796af37a9a..b1b08d59f1 100644 --- a/src/atom_operators.F +++ b/src/atom_operators.F @@ -600,7 +600,16 @@ CONTAINS ! add ecp_pseudo using numerical representation of basis IF (potential%ppot_type == ecp_pseudo) THEN ! scale 1/r potential - integrals%core = -potential%ecp_pot%zion*integrals%core + IF (ANY(potential%ecp_pot%nrloc(1:potential%ecp_pot%nloc) == 1)) THEN + alpha = 0.0_dp + DO k = 1, potential%ecp_pot%nloc + n = potential%ecp_pot%nrloc(k) + IF (n == 1) alpha = alpha + potential%ecp_pot%aloc(k) + END DO + integrals%core = (alpha - potential%ecp_pot%zion)*integrals%core + ELSE + integrals%core = -potential%ecp_pot%zion*integrals%core + END IF ! local potential m = basis%grid%nr rad => basis%grid%rad @@ -609,7 +618,11 @@ CONTAINS DO k = 1, potential%ecp_pot%nloc n = potential%ecp_pot%nrloc(k) alpha = potential%ecp_pot%bloc(k) - cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*EXP(-alpha*rad**2) + IF (n == 1) THEN + cpot(:) = cpot + potential%ecp_pot%aloc(k)/rad*(EXP(-alpha*rad**2) - 1.0_dp) + ELSE + cpot(:) = cpot + potential%ecp_pot%aloc(k)*rad**(n - 2)*EXP(-alpha*rad**2) + END IF END DO CALL numpot_matrix(integrals%core, cpot, basis, 0) ! non local pseudopotential diff --git a/tests/QS/regtest-ecp-2/TEST_FILES.toml b/tests/QS/regtest-ecp-2/TEST_FILES.toml index 317958c3b0..a95f7349da 100644 --- a/tests/QS/regtest-ecp-2/TEST_FILES.toml +++ b/tests/QS/regtest-ecp-2/TEST_FILES.toml @@ -1,3 +1,3 @@ "ICl_lanl2dz_gpw.inp" = [] -"Rn_stuttgart_gapw.inp" = [{matcher="M031", tol=1.0E-08, ref=3.56470376649E+03}] -"SbH3_def2_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-241.595963788403253}] +"Rn_stuttgart_gapw.inp" = [{matcher="M031", tol=1.0E-08, ref=3.58412821663E+03}] +"SbH3_def2_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-241.60143626280274}] diff --git a/tests/QS/regtest-ecp/TEST_FILES.toml b/tests/QS/regtest-ecp/TEST_FILES.toml index a00e7a171d..77e7b87d8e 100644 --- a/tests/QS/regtest-ecp/TEST_FILES.toml +++ b/tests/QS/regtest-ecp/TEST_FILES.toml @@ -1,5 +1,5 @@ -"ICl_lanl2dz_gpw.inp" = [{matcher="M011", tol=1.0E-11, ref=-23.515911094537550}] -"Rn_stuttgart_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-25.632304810373267}] +"ICl_lanl2dz_gpw.inp" = [{matcher="M011", tol=1.0E-11, ref=-26.20823873393772}] +"Rn_stuttgart_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-29.587239373669085}] "SbH3_def2_gapw.inp" = [{matcher="M011", tol=1.0E-11, ref=-241.414769303399083}] "HCl_ccECP.inp" = [{matcher="M011", tol=1.0E-12, ref=-15.457709431781252}] "HCl_force.inp" = []