Skip to content

Make method_init_density='pp' work for psp8 and add a generic atomic-density fallback #1260

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

Open
MZKC wants to merge 2 commits into
base: develop-2.0.0
Choose a base branch
from
Open

Make method_init_density='pp' work for psp8 and add a generic atomic-density fallback #1260

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c73acfe
feat: make method_init_density='pp' work for psp8 and add a generic f…
Jun 18, 2026
b29617d
fix: error out when the psp8 radial grid exceeds nrmax instead of sil…
MZKC Jun 24, 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
45 changes: 43 additions & 2 deletions src/atom/pp/input_pp.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,17 @@ subroutine input_pp(pp,hx,hy,hz)

if ( flag_beta_proj_is_given ) then
flag_potential_is_given=.false.
if(method_init_density/='read_dns_cube' .and. method_init_density/='wf' .and. ps_format(ik)/='UPF') then
stop "radial density is not available (method_init_density=pp...)"
if(method_init_density/='read_dns_cube' .and. method_init_density/='wf') then
! method_init_density='pp' needs a radial atomic density. It is read
! from the file for UPF (<PP_RHOATOM>) and psp8 (valence density). If
! the pseudopotential provides none (rho_pp_tbl is still all zero),
! fall back to a generic atomic density built from the valence charge
! so that 'pp' works for any pseudopotential instead of aborting.
if( maxval(abs(pp%rho_pp_tbl(:,ik))) < 1d-30 ) then
call set_generic_atomic_density(pp,ik)
write(*,'(a,i3,a)') " method_init_density=pp: pseudopotential of element ",ik, &
& " has no tabulated atomic density; using a generic atomic density."
end if
end if
end if
if ( any(pp%vpp_so/=0.0d0) ) flag_so=.true.
Expand Down Expand Up @@ -285,6 +294,26 @@ subroutine input_pp(pp,hx,hy,hz)
return
end subroutine input_pp
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130
subroutine set_generic_atomic_density(pp,ik)
use structures,only : s_pp_info
implicit none
type(s_pp_info),intent(inout) :: pp
integer,intent(in) :: ik
integer :: i
real(8) :: rs, r
! Generic spherical atomic valence density used as an initial-density guess
! (method_init_density='pp') when the pseudopotential provides no tabulated
! atomic density. Exponential model rho(r) = N*exp(-r/rs) normalized to the
! number of valence electrons (pp%zps), stored as 4*pi*r^2*rho:
! rho_pp_tbl(i) = (zps/(2*rs^3)) * r^2 * exp(-r/rs), with Int rho_pp_tbl dr = zps.
rs = 1.5d0 ! exponential scale length [a.u.]; this is only an initial guess
pp%rho_pp_tbl(:,ik) = 0.0d0
do i=1,pp%mr(ik)
r = pp%rad(i,ik)
pp%rho_pp_tbl(i,ik) = (dble(pp%zps(ik))/(2.0d0*rs**3)) * r**2 * exp(-r/rs)
end do
end subroutine set_generic_atomic_density
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130
subroutine read_ps_ky(pp,rrc,ik,ps_file)
use structures,only : s_pp_info
use salmon_global,only : Lmax_ps
Expand Down Expand Up @@ -575,9 +604,21 @@ subroutine read_ps_abinitpsp8(pp,rrc,rhor_nlcc,flag_nlcc_element,ik,ps_file)
write(*,*) "sum(rho_nlcc)=",sum_rho_nlcc*pi4*dr
end if
! valence charge density
! The psp8 file stores the valence density as 4*pi*rho (same convention as its
! NLCC data). Store it as 4*pi*r^2*rho in rho_pp_tbl (the convention used by
! calc_density_pp / method_init_density='pp'), i.e. r^2 * (4*pi*rho) = r^2*rho_tmp.
! The radial grid must fit rho_pp_tbl; otherwise the valence density would be
! silently truncated (and left inconsistent with sum_rho_pp below), corrupting
! the method_init_density='pp' initial density. Fail explicitly instead.
if (pp%mr(ik)+1 > pp%nrmax) then
write(*,*) "Error: psp8 radial grid size",pp%mr(ik)+1,"exceeds pp%nrmax",pp%nrmax, &
"in read_ps_abinitpsp8 (element",ik,")"
stop "read_ps_abinitpsp8: radial grid larger than pp%nrmax"
end if
sum_rho_pp=0.0d0
do i=1,pp%mr(ik)+1
read(4,*) dummy_text, r_tmp, rho_tmp
pp%rho_pp_tbl(i,ik) = pp%rad(i,ik)**2 * rho_tmp
sum_rho_pp=sum_rho_pp+rho_tmp*r_tmp**2
end do
write(*,*) "sum(rho_pp)=",sum_rho_pp*dr
Expand Down