Skip to content

Commit

Permalink
Add partial core corrections and tidy code
Browse files Browse the repository at this point in the history 
Note that partial core corrections (pcc) are often referred to as
non-linear core corrections (NLCC) in the literature
  • Loading branch information
@davidbowler
davidbowler committed Jun 24, 2024
1 parent 87d7751 commit 7ff3df2
Show file tree
Hide file tree
Showing 2 changed files with 91 additions and 26 deletions.
39 changes: 32 additions & 7 deletions tools/BasisGeneration/read_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -515,7 +515,7 @@ subroutine read_hgh_input(i_species)
character(len=80) :: line
logical :: flag_core_done = .false.
logical, dimension(3,0:3) :: flag_min
real(double) :: dummy, dummy2, highest_energy, root_two, proj, rr_lp, pj, pjp
real(double) :: dummy, dummy2, highest_energy, root_two, proj, rr_lp, pj, pjp, r_core, r_core_2, c_core
real(double) :: rl_base, rl_sqrt, rr, rr_l, rr_rl, rr_rl2, rr_rl4, rr_rl6, charge
real(double), dimension(:,:), allocatable :: gamma_fac
real(double), dimension(:,:), allocatable :: hnl
Expand Down Expand Up @@ -659,12 +659,22 @@ subroutine read_hgh_input(i_species)
write(*,fmt='("n, l and occupancy: ",i1," ",i1,f6.2)') val%n(i), val%l(i), val%occ(i)
end do
val%n_occ = n_occ
! Test for PCC
ios = 0
r_core = zero
c_core = zero
n_read = 0 ! Compatibility with CP2K files; not used
read(lun,*,iostat=ios) r_core, n_read, c_core
if(ios==0) then
pseudo(i_species)%flag_pcc = .true.
write(*,fmt='("This pseudopotential includes partial core corrections")')
end if
call io_close(lun)
!
! Grid size
!
ngrid = log(45.0_double/(beta/pseudo(i_species)%z))/log(1.012_double) ! Following Hamann
write(*,*) 'Number of grid points ',ngrid
if(iprint>2) write(*,fmt='("Number of grid points ",i5)') ngrid
call allocate_vkb(ngrid,i_species)
! Assign pseudo-n value for nodes
do i = 1,n_shells
Expand All @@ -681,15 +691,13 @@ subroutine read_hgh_input(i_species)
end do
if(iprint>3) write(*,fmt='(i2," valence shells, with ",i2," occupied")') n_shells, n_occ
root_two = sqrt(two)
! Read density or set to zero
! Read density from charge.dat or set to zero if file not present
open(unit=40,file='charge.dat',status='old',iostat=ios)
if(ios==0) then
!charge = zero
do i=1,ngrid
read(40,*) rr,local_and_vkb%charge(i)
charge = charge + alpha*rr*rr*rr*local_and_vkb%charge(i)
! if(abs(rr-(beta/pseudo(i_species)%z)*exp(alpha*real(i-1,double)))>1e-5_double) &
! write(*,*) 'Mesh error: ',rr,(beta/pseudo(i_species)%z)*exp(alpha*real(i-1,double))
end do
close(40)
if(iprint>4) write(*,*) 'Charge read in integrates to : ',charge
Expand All @@ -713,6 +721,19 @@ subroutine read_hgh_input(i_species)
local_and_vkb%r_vkb = local_and_vkb%rr(ngrid)
local_and_vkb%ngrid_vkb = ngrid
!
! Create PCC
!
if(pseudo(i_species)%flag_pcc) then
allocate(local_and_vkb%pcc(ngrid))
local_and_vkb%pcc = zero
! There is a parameter n_read above which is currently unused
r_core_2 = r_core*r_core
do i=1,ngrid
rr = local_and_vkb%rr(i)
local_and_vkb%pcc(i) = c_core*exp(-half*rr*rr/r_core_2) ! May need to remove factor of 4pi
end do
end if
!
! Create non-local projectors
!
! i from 1 to 3
Expand All @@ -734,21 +755,25 @@ subroutine read_hgh_input(i_species)
n_r_proj_max = 1
flag_min = .true.
do i=1,ngrid
if(local_and_vkb%rr(i)>half) then
do ell = 0,max_l
if(local_and_vkb%n_proj(ell)>0) then
rr_rl = local_and_vkb%rr(i)/hgh_data(i_species)%r(ell)
do j = 1,local_and_vkb%n_proj(ell)!3 ! Fix later to account for number of projectors per l
do j = 1,local_and_vkb%n_proj(ell)
rr_l = local_and_vkb%rr(i)**(ell + 2*(j-1))
proj = root_two*rr_l*exp(-0.5*rr_rl*rr_rl)/gamma_fac(j,ell)
if(abs(proj)<1e-8.and.flag_min(j,ell)) then
if(local_and_vkb%rr(i)>local_and_vkb%core_radius(ell)) local_and_vkb%core_radius(ell) = local_and_vkb%rr(i)
if(local_and_vkb%rr(i)>local_and_vkb%rr(n_r_proj_max)) n_r_proj_max = i
flag_min(j,ell) = .false.
write(*,'("l=",i1," core radius ",f6.3," bohr")') ell, local_and_vkb%core_radius(ell)
end if
end do
end if
end do
end if
end do
do ell = 0,max_l
write(*,'("l=",i1," core radius ",f6.3," bohr")') ell, local_and_vkb%core_radius(ell)
end do
! Now calculate projectors
local_and_vkb%r_vkb = local_and_vkb%rr(n_r_proj_max)
Expand Down
78 changes: 59 additions & 19 deletions tools/BasisGeneration/schro_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,12 @@ subroutine find_unconfined_valence_states(i_species,vha,vxc)
iter = 0
if(iprint>2) write(*,fmt='(/2x,"Finding unconfined energies for valence states")')
resid = one
!
! SCF iteration needed for HGH where we do not have an atomic density supplied
! I have found that if we start from no charge density then it's better to allow
! the charge to gradually change towards the correct ionic charge rather than
! rescaling after the first iteration
!
do while(resid>1e-12_double.and.iter<maxiter)
iter = iter+1
newcharge = zero
Expand All @@ -124,11 +130,11 @@ subroutine find_unconfined_valence_states(i_species,vha,vxc)
en = val%npao(i_shell)
large_energy = val%en_ps(i_shell)
if(ps_format==hgh) then
if(resid<0.01_double) then
large_energy = val%en_ps(i_shell)
else
!if(resid<0.01_double) then
! large_energy = val%en_ps(i_shell)
!else
large_energy = zero!half*val%en_ps(i_shell)
end if
!end if
end if
!if(abs(large_energy)<RD_ERR.and.i_shell>1) large_energy = val%en_pao(i_shell-1)
call find_eigenstate_and_energy_vkb(i_species,en,ell,radius_large, psi,large_energy,vha,vxc)
Expand All @@ -149,11 +155,12 @@ subroutine find_unconfined_valence_states(i_species,vha,vxc)
resid = resid + drdi(i)*rr_squared(i)*(local_and_vkb%charge(i)-newcharge(i))**2
check = check + drdi(i)*rr_squared(i)*local_and_vkb%charge(i)
end do
! Simple linear mixing
local_and_vkb%charge = alpha_scf*newcharge + (one-alpha_scf)*local_and_vkb%charge
! We can use these lines to write out the charge for future solvers
!open(unit=70,file='charge_out.dat',position="append")
!do i=1,nmesh
! write(70,*) rr(i), newcharge(i)!local_and_vkb%charge(i)
! write(70,*) rr(i), newcharge(i)
!end do
!close(70)
! Integrate
Expand All @@ -167,17 +174,7 @@ subroutine find_unconfined_valence_states(i_species,vha,vxc)
! This rescales charge to have full valence value, but can be unstable
!if(abs(check-total_charge)>RD_ERR) local_and_vkb%charge = local_and_vkb%charge*total_charge/check
if(iprint>2) write(*,fmt='("Iteration ",i4," residual ",e14.6)') iter,resid
if(iprint>0) then
write(*,fmt='(/2x,"Unconfined valence state energies (Ha)")')
write(*,fmt='(2x," n l AE energy PAO energy")')
do i_shell = 1, val%n_occ
ell = val%l(i_shell)
en = val%n(i_shell)
write(*,fmt='(2x,2i3,2f18.10)') en, ell, val%en_ps(i_shell), &
val%en_pao(i_shell)
if(ps_format==hgh) val%en_ps(i_shell) = val%en_pao(i_shell)
end do
else if(ps_format==hgh) then
if(ps_format==hgh) then
do i_shell = 1, val%n_occ
val%en_ps(i_shell) = val%en_pao(i_shell)
end do
Expand All @@ -187,6 +184,25 @@ subroutine find_unconfined_valence_states(i_species,vha,vxc)
call get_vxc(nmesh,rr,local_and_vkb%charge,vxc)
if(pseudo(i_species)%flag_pcc) local_and_vkb%charge = local_and_vkb%charge - local_and_vkb%pcc
end do
if(iprint>0) then
write(*,fmt='(/2x,"Unconfined valence state energies (Ha)")')
if(ps_format==hgh) then
write(*,fmt='(2x," n l PAO energy")')
do i_shell = 1, val%n_occ
ell = val%l(i_shell)
en = val%n(i_shell)
write(*,fmt='(2x,2i3,f18.10)') en, ell,val%en_pao(i_shell)
end do
else
write(*,fmt='(2x," n l AE energy PAO energy")')
do i_shell = 1, val%n_occ
ell = val%l(i_shell)
en = val%n(i_shell)
write(*,fmt='(2x,2i3,2f18.10)') en, ell, val%en_ps(i_shell), &
val%en_pao(i_shell)
end do
end if
end if
if(iter>maxiter) call cq_abort("Exceeded iterations in SCF")
! We can use these lines to write out the charge for future solvers
open(unit=70,file='charge_out.dat')
Expand Down Expand Up @@ -716,7 +732,7 @@ subroutine find_default_cutoffs(i_species,vha,vxc)
deltaE_large_radius, deltaE_large_radius*HaToeV
small_cutoff(i_shell) = large_cutoff(i_shell)
end if
write(*,*) '# Radii: ',large_cutoff(i_shell),small_cutoff(i_shell)
if(iprint>3) write(*,*) '# Radii: ',large_cutoff(i_shell),small_cutoff(i_shell)
end do
! Create cutoffs based on defaults chosen by user
if(paos%flag_cutoff==pao_cutoff_energies.OR.paos%flag_cutoff==pao_cutoff_default) then ! Same energy for all l/n shells
Expand Down Expand Up @@ -912,7 +928,7 @@ subroutine find_eigenstate_and_energy_vkb(i_species,en,ell,Rc,psi,energy,vha,vxc
real(double), OPTIONAL :: width, prefac

! Local variables
real(double) :: g_temp, dy_L, dy_R
real(double) :: g_temp, dy_L, dy_R, ipsi_in, ipsi_out
real(double), dimension(:), allocatable :: f, potential
integer :: classical_tp, i, n_crossings, n_nodes, n_loop, loop, nmax, n_kink, n_nodes_lower, n_nodes_upper, n_kink_vkb
real(double) :: l_half_sq, dx_sq_over_twelve, fac, norm, d_energy, e_lower, e_upper, df_cusp, cusp_psi, tol
Expand Down Expand Up @@ -1041,7 +1057,7 @@ subroutine find_eigenstate_and_energy_vkb(i_species,en,ell,Rc,psi,energy,vha,vxc
psi(n_kink:nmax) = psi(n_kink:nmax)*fac
xin = psi(n_kink)*f(n_kink)- psi(n_kink+1)*f(n_kink+1)
gin = psi(n_kink)
gsgin = psi(n_kink)*psi(n_kink)*drdi_squared(n_kink)
gsgin = psi(n_kink)*psi(n_kink)*drdi_squared(n_kink)
! Remember that psi is y in numerov - don't forget factor of root(r)
! Normalise
norm = zero
Expand Down Expand Up @@ -1081,6 +1097,30 @@ subroutine find_eigenstate_and_energy_vkb(i_species,en,ell,Rc,psi,energy,vha,vxc
!d_energy = cusp_psi
if(iprint>4) write(*,fmt='(2x,"Energy shift (alt): ",f18.10)') cusp_psi
if(iprint>5) write(*,fmt='(2x,"Number of nodes: ",i4)') n_crossings
!! Integrate to kink in both directions - this is another estimate of the
!! energy change required but is rather approximate because of the method
!! used to calculate dpsi/dr
!ipsi_out = zero
!do i=1,n_kink
! ipsi_out = ipsi_out + psi(i)*psi(i)*drdi_squared(i)
!end do
!ipsi_out = ipsi_out/(psi(n_kink)*psi(n_kink))
!ipsi_in = zero
!do i=n_kink,nmax
! ipsi_in = ipsi_in + psi(i)*psi(i)*drdi_squared(i)
!end do
!ipsi_in = ipsi_in/(psi(n_kink)*psi(n_kink))
!!dy_L = (psi(n_kink)*sqrt(drdi(n_kink))/rr(n_kink)-psi(n_kink-1)*sqrt(drdi(n_kink-1))/rr(n_kink-1)) &
!! /(rr(n_kink)-rr(n_kink-1))
!!dy_R = (psi(n_kink+1)*sqrt(drdi(n_kink+1))/rr(n_kink+1)-psi(n_kink)*sqrt(drdi(n_kink))/rr(n_kink)) &
!! /(rr(n_kink+1)-rr(n_kink))
!dy_L = (psi(n_kink)*drdi(n_kink)-psi(n_kink-1)*drdi(n_kink-1)) &
! /(rr(n_kink)-rr(n_kink-1))
!dy_R = (psi(n_kink+1)*drdi(n_kink+1)-psi(n_kink)*drdi(n_kink)) &
! /(rr(n_kink+1)-rr(n_kink))
!d_energy = -(dy_L/psi(n_kink) - dy_R/psi(n_kink))/(ipsi_in + ipsi_out)
!write(*,*) 'New d_energy is ',d_energy
! Write out psi here?
if ( n_crossings /= n_nodes) then
if ( n_crossings > n_nodes ) then
e_upper = energy
Expand Down

0 comments on commit 7ff3df2

Please sign in to comment.