diff --git a/pseudo-and-pao/LDA/test_default_PAOs.sh b/pseudo-and-pao/LDA/test_default_PAOs.sh new file mode 100755 index 00000000..4c898a8e --- /dev/null +++ b/pseudo-and-pao/LDA/test_default_PAOs.sh @@ -0,0 +1,22 @@ +#!/bin/bash +echo "Testing default PAO generation" +declare -i success=0 +declare -i failure=0 +declare -i total=0 +for ele in Ag As Ba Br Cd Cr F Ge Hf In Kr Mn Na Ni P Po Re Ru Sc Sn Tc Tl Xe Zr Al Au Be C Cl Cs Fe H Hg Ir Li Mo Nb O Pb Pt Rh S Se Sr Te V Y Ar B Bi Ca Co Cu Ga He I K Mg N Ne Os Pd Rb Rn Sb Si Ta Ti W Zn +do + let total++ + cd $ele + ../../../bin/MakeIonFiles > output 2>&1 + if [ `tail -1 output | grep -c Finish` -eq 1 ]; then + let success++ + rm output "$ele"CQ.ion input.log + else + let failure++ + echo "Possible failure in default PAO generation for "$ele". Please report to developers." + rm input.log + fi + cd .. +done +echo "Successes: "$success"/"$total +echo "Failures: "$failure"/"$total diff --git a/pseudo-and-pao/PBE/test_default_PAOs.sh b/pseudo-and-pao/PBE/test_default_PAOs.sh new file mode 100755 index 00000000..37a6d58d --- /dev/null +++ b/pseudo-and-pao/PBE/test_default_PAOs.sh @@ -0,0 +1,22 @@ +#!/bin/bash +echo "Testing default PAO generation" +declare -i success=0 +declare -i failure=0 +declare -i total=0 +for ele in Ag As Ba Br Cd Cr F Ge Hf In Kr Lu Mo Nb O Pb Pt Rh S Se Sr Te V Y Al Au Be C Cl Cs Fe H Hg Ir La Mg N Ne Os Pd Rb Rn Sb Si Ta Ti W Zn Ar B Bi Ca Co Cu Ga He I K Li Mn Na Ni P Po Re Ru Sc Sn Tc Tl Xe Zr +do + let total++ + cd $ele + ../../../bin/MakeIonFiles > output 2>&1 + if [ `tail -1 output | grep -c Finish` -eq 1 ]; then + let success++ + rm output "$ele"CQ.ion input.log + else + let failure++ + echo "Possible failure in default PAO generation for "$ele". Please report to developers." + rm input.log + fi + cd .. +done +echo "Successes: "$success"/"$total +echo "Failures: "$failure"/"$total diff --git a/pseudo-and-pao/PBEsol/test_default_PAOs.sh b/pseudo-and-pao/PBEsol/test_default_PAOs.sh new file mode 100755 index 00000000..37a6d58d --- /dev/null +++ b/pseudo-and-pao/PBEsol/test_default_PAOs.sh @@ -0,0 +1,22 @@ +#!/bin/bash +echo "Testing default PAO generation" +declare -i success=0 +declare -i failure=0 +declare -i total=0 +for ele in Ag As Ba Br Cd Cr F Ge Hf In Kr Lu Mo Nb O Pb Pt Rh S Se Sr Te V Y Al Au Be C Cl Cs Fe H Hg Ir La Mg N Ne Os Pd Rb Rn Sb Si Ta Ti W Zn Ar B Bi Ca Co Cu Ga He I K Li Mn Na Ni P Po Re Ru Sc Sn Tc Tl Xe Zr +do + let total++ + cd $ele + ../../../bin/MakeIonFiles > output 2>&1 + if [ `tail -1 output | grep -c Finish` -eq 1 ]; then + let success++ + rm output "$ele"CQ.ion input.log + else + let failure++ + echo "Possible failure in default PAO generation for "$ele". Please report to developers." + rm input.log + fi + cd .. +done +echo "Successes: "$success"/"$total +echo "Failures: "$failure"/"$total diff --git a/tools/BasisGeneration/read_module.f90 b/tools/BasisGeneration/read_module.f90 index 8843ac53..120824f1 100644 --- a/tools/BasisGeneration/read_module.f90 +++ b/tools/BasisGeneration/read_module.f90 @@ -420,9 +420,7 @@ contains ! ! Li, Be (npao set to one because l=0) if(paos%l(paos%n_shells)==1 .AND. paos%inner(paos%n_shells)==0) paos%npao(paos%n_shells) = 1 - ! We need sign change if s orbital is outermost and there is a semi-core p orbital - if(paos%l(paos%n_shells)==1 .AND. paos%polarised_shell==val%n_occ & - .AND. paos%inner(paos%n_shells)>0) paos%pol_pf = -one + ! Remove pol_pf here as it is now incorporated into find_polarisation DRB 2021/09/02 ! Use Vl if perturbing a p orbital with a valence or semi-core d orbital if(paos%l(paos%n_shells)==2) then do i=1,val%n_occ diff --git a/tools/BasisGeneration/schro_module.f90 b/tools/BasisGeneration/schro_module.f90 index 39cc83d6..2fb42d67 100644 --- a/tools/BasisGeneration/schro_module.f90 +++ b/tools/BasisGeneration/schro_module.f90 @@ -369,26 +369,27 @@ contains ! if(paos%flag_perturb_polarise) then i_shell = paos%n_shells - ell = paos%l(i_shell) + ! NB we pass l-1 here as the angular momentum of the shell being perturbed + ell = paos%l(i_shell)-1 + ! Note that npao is set to the value of n for the shell being perturbed en = paos%npao(i_shell) do zeta = 1,paos%nzeta(i_shell) allocate(paos%psi(zeta,i_shell)%f(nmesh)) paos%psi(zeta,i_shell)%f = zero - ! NB we pass l-1 here as we need to perturb the shell below - ell = paos%l(i_shell)-1 - en = paos%npao(i_shell) if(iprint>2) then write(*,fmt='(2x,"Perturbative polarisation")') write(*,fmt='(2x,"Species ",i2," n=",i2," l=",i2," zeta=",i2, " Rc=",f4.1," bohr")') & i_species, en, ell, zeta, paos%cutoff(zeta,i_shell) write(*,fmt='(4x,"Prefactor scaled by ",f5.1)') paos%pol_pf + write(*,fmt='(4x,"Energy: ",f12.5)') paos%energy(zeta,paos%polarised_shell) + write(*,fmt='(4x,"Polarising shell: ",i2)') paos%polarised_shell end if if(zeta>1.AND.paos%flag_zetas==1) then call find_split_norm(en,ell+1,paos%cutoff(zeta,i_shell),& paos%psi(zeta,i_shell)%f,paos%psi(1,i_shell)%f) else ! Set radius of polarisation PAO - paos%cutoff(zeta,i_shell) = paos%cutoff(zeta,paos%polarised_shell) !i_shell-1) + paos%cutoff(zeta,i_shell) = paos%cutoff(zeta,paos%polarised_shell) call find_polarisation(i_species,en,ell,paos%cutoff(zeta,i_shell),& paos%psi(zeta,paos%polarised_shell)%f,paos%psi(zeta,i_shell)%f,& paos%energy(zeta,paos%polarised_shell),vha,vxc,paos%pol_pf) @@ -1217,8 +1218,14 @@ contains l_half_sq = real(ell+1,double) + half l_half_sq = l_half_sq*l_half_sq l_l_plus_one = real((ell+2)*(ell+1),double) + ! NB The values of n and l we have here are for the perturbed function, so this is correct n_nodes = en - ell - 1 if(iprint>2) write(*,fmt='(2x,"For this polarisation function, we require ",i2," nodes")') n_nodes + ! This is needed when we are perturbing a valence state with a semi-core state + if(n_nodes>0 .and. pf_sign>zero .and. psi_l(1)5) write(*,fmt='(4x,"Changing sign of function as required")') + end if zval = pseudo(i_species)%z dx_sq_over_twelve = alpha*alpha/twelve alpha_sq_over_four = alpha*alpha/four @@ -1309,10 +1316,12 @@ contains end if end do if(loop>=n_loop) then - do i=1,nmax - write(50,*) rr(i),psi_pol(i),-rr(i)*rr(i)*psi_l(i) - end do - call flush(50) + if(iprint>5) then + do i=1,nmax + write(50,*) rr(i),psi_pol(i),-rr(i)*rr(i)*psi_l(i) + end do + call flush(50) + end if call cq_abort("ERROR in perturbative polarisation routines - prefactor not found") else if(iprint>2) write(*,fmt='("Finished perturbative polarisation search with prefactor ",f13.5)') prefac