1! 2! Copyright (C) 2004 PWSCF group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!----------------------------------------------------------------------- 10subroutine set_rho_core 11 !----------------------------------------------------------------------- 12 ! 13 ! input : all-electron wavefunctions + valence states 14 ! output: smoothed core charge for r > rcore 15 ! 16 use kinds, only : dp 17 use constants, only : pi 18 use io_global, only : stdout, ionode, ionode_id 19 use mp, only : mp_bcast 20 use mp_world, only : world_comm 21 use ld1inc, only : nlcc, grid, rhoc, aeccharge, psccharge, rcore, & 22 nwf, oc, rel, core_state, psi, file_core, new_core_ps,& 23 lpaw, lnc2paw 24 implicit none 25 26 real(DP) :: drho, const, br1, br2, & 27 eps1, eps2, br12, xc(8), a, b, eps12, totrho 28 real(DP), allocatable:: rhov(:) 29 real(DP), external :: int_0_inf_dr 30 integer :: i, ik, n, ns, ios 31 32 if (nlcc) then 33 write(stdout,'(/,5x,'' Computing core charge for nlcc: '')') 34 else 35 if (lpaw) write(stdout,'(/,5x,'' Computing core charge for PAW: '')') 36 end if 37 allocate (rhov(grid%mesh)) 38 ! 39 ! calculates core charge density 40 ! 41 do n=1,grid%mesh 42 rhov(n) = 0.0_dp 43 rhoc(n) = 0.0_dp 44 do ns=1,nwf 45 if (rel==2) then 46 if (core_state(ns)) then 47 rhoc(n)=rhoc(n)+oc(ns)*(psi(n,1,ns)**2+psi(n,2,ns)**2) 48 else 49 rhov(n)=rhov(n)+oc(ns)*(psi(n,1,ns)**2+psi(n,2,ns)**2) 50 endif 51 else 52 if (core_state(ns)) then 53 rhoc(n) = rhoc(n) + oc(ns)*psi(n,1,ns)**2 54 else 55 rhov(n) = rhov(n) + oc(ns)*psi(n,1,ns)**2 56 endif 57 endif 58 enddo 59 enddo 60 totrho = int_0_inf_dr(rhoc,grid,grid%mesh,2) 61 if (totrho<1.d-6.and.lpaw) then 62! 63! All valence charge for this atom (mainly for H) 64! 65 aeccharge(1:grid%mesh) = 0.0_DP 66 psccharge(1:grid%mesh) = 0.0_DP 67 goto 1100 68 endif 69 70! write(stdout,'("Integrated core charge",f15.10)') totrho 71 aeccharge(1:grid%mesh) = rhoc(1:grid%mesh) 72 ! 73 if (rcore > 0.0_dp) then 74 ! rcore read on input 75 do ik=1,grid%mesh 76 if (grid%r(ik) > rcore) go to 100 77 enddo 78 call infomsg('set_rho_core','rcore too big') 79 return 80 else 81 ! rcore determined by the condition rhoc(rcore) = 2*rhov(rcore) 82 do ik=1,grid%mesh 83 if (rhoc(ik) < 2.0 * rhov(ik)) go to 100 84 enddo 85 end if 86100 rcore=grid%r(ik) 87 drho = ( rhoc(ik+1)/grid%r2(ik+1) - rhoc(ik)/grid%r2(ik) ) / grid%dx / grid%r(ik) 88 ! 89 ! true_rho = rhoc(r)/r**2/4 pi 90 ! (factor 1/r from logarithmic mesh) 91 ! smoothened core charge density for nonlinear core correction: 92 ! rhoc(r) = core charge for r > rcore 93 ! rhoc(r) = r^2 a sin(b r)/r for r < rcore 94 ! 95 if (new_core_ps) then 96 call compute_phius(1,ik,aeccharge,rhoc,xc,0,' ') 97 else 98 if (drho > 0.0_dp) then 99 call infomsg('set_rho_core','d rho/ d r > 0') 100 return 101 endif 102 const= grid%r(ik)*drho / ( rhoc(ik)/grid%r2(ik) ) + 1.0_dp 103 if (const > 0.0_dp) then 104 br1 = 0.00001_dp 105 br2 = pi/2.0_dp-0.00001_dp 106 else 107 br1 = pi/2.0_dp+0.00001_dp 108 br2 = pi 109 end if 110 do n=1, 15 111 eps1 = br1 - const*tan(br1) 112 eps2 = br2 - const*tan(br2) 113 br12 = (br1+br2)/2.0_dp 114 eps12 = br12 - const*tan(br12) 115 if(eps1*eps12 < 0.0_dp) then 116 br2 = br12 117 else if(eps12*eps2 < 0.0_dp) then 118 br1 = br12 119 else 120 call errore('set_rho_core','error in bisection',n) 121 end if 122 end do 123 b = br12/grid%r(ik) 124 a = ( rhoc(ik)/grid%r2(ik) ) * grid%r(ik)/sin(br12) 125 do n=1,ik 126 rhoc(n) = a*sin(b*grid%r(n))/grid%r(n) * grid%r2(n) 127 end do 128 endif 129 if (lpaw) then 130 if (lnc2paw) then 131 ! Mimic NC calculation. If NLCC, the pseudized core charge. 132 if (nlcc) then 133 aeccharge(1:grid%mesh) = rhoc(1:grid%mesh) 134 ! Here one could set another pseudized ccharge, for 135 ! example with a larger matching radius. Right now, 136 ! just take the same as the AE (ie NC) one: 137 psccharge(1:grid%mesh) = rhoc(1:grid%mesh) 138 else 139 ! Reference NC calculation does not have core charge. 140 aeccharge(1:grid%mesh) = 0._dp 141 psccharge(1:grid%mesh) = 0._dp 142 end if 143 else 144! aeccharge(1:grid%mesh) = rhoc(1:grid%mesh) 145 if (nlcc) then 146 psccharge(1:grid%mesh) = rhoc(1:grid%mesh) 147 else 148 psccharge(1:grid%mesh) = 0._dp 149 end if 150 end if 151 end if 152 write(stdout,'(/,5x,'' r > '',f4.2,'' : true rho core'')') grid%r(ik) 153 if (new_core_ps) then 154 write(stdout,'(6x,"Core charge pseudized with two Bessel functions")') 155 else 156 write(stdout,110) grid%r(ik), a, b 157 endif 158110 format (5x, ' r < ',f4.2,' : rho core = a sin(br)/r', & 159 ' a=',f7.2,' b=',f7.2/) 1601100 continue 161 if (file_core .ne. ' ') then 162 write(stdout,'(6x, "***Writing file ",a, " ***")') trim(file_core) 163 if (ionode) & 164 open(unit=26,file=file_core, status='unknown', iostat=ios, err=300 ) 165300 call mp_bcast(ios, ionode_id, world_comm) 166 call errore('set_rho_core','opening file '//file_core,abs(ios)) 167 if (ionode) then 168 if (totrho>1.d-6) then 169 do n=1,grid%mesh 170 write(26,'(4f20.10)') grid%r(n),rhoc(n),rhov(n),aeccharge(n) 171 enddo 172 else 173 do n=1,grid%mesh 174 write(26,'(2f20.10)') grid%r(n),rhov(n) 175 enddo 176 endif 177 close(26) 178 endif 179 endif 180 totrho = int_0_inf_dr(rhoc,grid,grid%mesh,2) 181 write(stdout,'(6x,''Integrated core pseudo-charge : '',f6.2)') totrho 182 if (.not.nlcc) rhoc(1:grid%mesh) = 0.0_dp 183 deallocate (rhov) 184 return 185end subroutine set_rho_core 186