TABLE OF CONTENTS
ABINIT/integrho [ Functions ]
NAME
integrho
FUNCTION
This routine integrates the electron density inside the atomic surface already calculated - it reads the file *.surf The radial integration is always performed with splines and the two angular integrations with Gauss quadrature
COPYRIGHT
Copyright (C) 2002-2017 ABINIT group (PCasek,FF,XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
aim_dtset = the structured entity containing all input variables znucl_batom=the nuclear charge of the Bader atom
OUTPUT
(see side effects)
SIDE EFFECTS
This routine works primarily on the data contained in the aimfields and aimprom modules WARNING This file does not follow the ABINIT coding rules (yet)
PARENTS
drvaim
CHILDREN
bschg1,coeffs_gausslegint,spline,vgh_rho
SOURCE
40 #if defined HAVE_CONFIG_H 41 #include "config.h" 42 #endif 43 44 #include "abi_common.h" 45 46 subroutine integrho(aim_dtset,znucl_batom) 47 48 use defs_basis 49 use defs_aimfields 50 use defs_aimprom 51 use defs_parameters 52 use defs_abitypes 53 use m_profiling_abi 54 use m_splines 55 use m_errors 56 57 use m_numeric_tools, only : coeffs_gausslegint 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'integrho' 63 use interfaces_63_bader, except_this_one => integrho 64 !End of the abilint section 65 66 implicit none 67 68 !Arguments ------------------------------------ 69 !scalars 70 type(aim_dataset_type),intent(in) :: aim_dtset 71 72 !Local variables ------------------------------ 73 !scalars 74 integer :: batom,chs,iat,ii,inx,inxf,ipos,jj,kk,ll,nn,nph,nth 75 real(dp) :: chg,chgint,cintr,ct1,ct2,lder,nsphe,phimax,phimin,rder 76 real(dp) :: rsmax,rsmin,ss,stp,themax,themin,uu 77 real(dp) :: znucl_batom,zz 78 logical :: gaus,weit 79 !arrays 80 real(dp) :: grho(3),hrho(3,3),shift(3),unvec(3),vv(3) 81 real(dp),allocatable :: ncrho(:),nsp2(:),nsp3(:),nsp4(:),rdint(:,:),rr(:) 82 real(dp),allocatable :: vdd(:),vrho(:),wgrs(:,:),work(:) 83 84 ! ********************************************************************* 85 86 gaus=.true. 87 weit=.true. 88 89 write(std_out,*) 'npt = ',aim_dtset%npt 90 91 rewind(unts) 92 read(unts,*) batom,shift ! Warning : batom is read, instead of coming from aim_dtset 93 read(unts,*) nth,themin,themax ! Warning : these numbers are read, instead of coming from aim_dtset 94 read(unts,*) nph,phimin,phimax ! Warning : these numbers are read, instead of coming from aim_dtset 95 96 write(std_out,*) 'NTH NPH ',nth,nph 97 98 ABI_ALLOCATE(wgrs,(nth,nph)) 99 ABI_ALLOCATE(rdint,(nth,nph)) 100 101 do ii=1,nth 102 do jj=1,nph 103 if (weit) then 104 read(unts,*) th(ii),ph(jj),rs(ii,jj),wgrs(ii,jj) 105 else 106 read(unts,*) th(ii),ph(jj),rs(ii,jj) 107 end if 108 end do 109 end do 110 read(unts,*) rsmin,rsmax 111 112 113 if (gaus) then 114 ct1=cos(themin) 115 ct2=cos(themax) 116 call coeffs_gausslegint(ct1,ct2,cth,wcth,nth) 117 call coeffs_gausslegint(phimin,phimax,ph,wph,nph) 118 end if 119 120 do ii=1,nth 121 do jj=1,nph 122 if (.not.weit) then 123 if (gaus) then 124 wgrs(ii,jj)=wcth(ii)*wph(jj) 125 else 126 wgrs(ii,jj)=1._dp 127 end if 128 end if 129 end do 130 end do 131 132 133 do ii=1,nth 134 do jj=1,nph 135 if (rs(ii,jj) < rsmin) rsmin=rs(ii,jj) 136 end do 137 end do 138 139 140 !INTEGRATION OF THE CORE DENSITY 141 142 nn=typat(batom) 143 kk=ndat(nn) 144 145 146 !spherical integration of the core density in the sphere 147 !of the minimal Bader radius 148 149 !COEF. FOR SPHERICAL INTEGRATION 150 151 ABI_ALLOCATE(nsp2,(kk)) 152 ABI_ALLOCATE(nsp3,(kk)) 153 ABI_ALLOCATE(nsp4,(kk)) 154 ABI_ALLOCATE(ncrho,(kk)) 155 156 do ii=1,kk 157 ncrho(ii)=crho(ii,nn)*4._dp*pi*rrad(ii,nn)*rrad(ii,nn) 158 nsp3(ii)=4._dp*pi*(2._dp*crho(ii,nn)+2._dp*rrad(ii,nn)*sp2(ii,nn)+& 159 & rrad(ii,nn)*rrad(ii,nn)*sp3(ii,nn)) 160 end do 161 162 if (rsmin < rrad(ndat(nn),nn)) then ! search index 163 inx=0 164 if (rsmin < rrad(1,nn)) then 165 MSG_ERROR('absurd') 166 elseif (rsmin > rrad(ndat(nn),nn)) then 167 inx=ndat(nn) 168 else 169 do while (rsmin >= rrad(inx+1,nn)) 170 inx=inx+1 171 end do 172 end if 173 else 174 inx=ndat(nn) 175 end if 176 177 cintr=4._dp/3._dp*pi*rrad(1,nn)**3*crho(1,nn) 178 179 !spline integration 180 181 do ii=1,inx-1 182 uu=rrad(ii+1,nn)-rrad(ii,nn) 183 cintr=cintr+(ncrho(ii)+ncrho(ii+1))*uu/2._dp-uu*uu*uu/2.4d1*(nsp3(ii)+nsp3(ii+1)) 184 end do 185 if (inx/=ndat(nn)) then 186 uu=rsmin-rrad(inx,nn) 187 zz=rrad(inx+1,nn)-rsmin 188 ss=rrad(inx+1,nn)-rrad(inx,nn) 189 cintr=cintr+ncrho(inx)/2._dp*(ss-zz*zz/ss)+ncrho(inx+1)/2._dp*uu*uu/ss+& 190 nsp3(inx)/1.2d1*(zz*zz*ss-zz*zz*zz*zz/2._dp/ss-ss*ss*ss/2._dp)+& 191 nsp3(inx+1)/1.2d1*(uu*uu*uu*uu/2._dp/ss-uu*uu*ss) 192 end if 193 194 195 !INTEGRATION OF THE REST OF THE CORE DENSITY 196 !(for gauss quadrature) 197 !For the Gauss quadrature it is added 198 !to the radial integrated valence density 199 200 rdint(:,:)=0._dp 201 nsphe=0._dp 202 do ii=1,nth 203 do jj=1,nph 204 if (inx==ndat(nn)) cycle 205 inxf=inx 206 if (rs(ii,jj) < rsmin) then 207 write(std_out,*) rs(ii,jj),rsmin 208 MSG_ERROR('in surface') 209 elseif (rs(ii,jj) > rrad(ndat(nn),nn)) then 210 inxf=ndat(nn) 211 else 212 do while (rs(ii,jj) >= rrad(inxf+1,nn)) 213 inxf=inxf+1 214 end do 215 end if 216 217 if (inxf==inx) then 218 uu=rrad(inx+1,nn)-rs(ii,jj) 219 zz=rrad(inx+1,nn)-rsmin 220 ss=rrad(inx+1,nn)-rrad(inx,nn) 221 222 rdint(ii,jj)=(ncrho(inx)/2._dp/ss-nsp3(inx)/1.2d1*ss)*(zz*zz-uu*uu)+& 223 nsp3(inx)/2.4d1/ss*(zz**4-uu**4) 224 uu=rs(ii,jj)-rrad(inx,nn) 225 zz=rsmin-rrad(inx,nn) 226 rdint(ii,jj)=rdint(ii,jj)+(uu*uu-zz*zz)*(ncrho(inx+1)/2._dp/ss-nsp3(inx+1)/1.2d1*ss)+& 227 nsp3(inx+1)/2.4d1/ss*(uu**4-zz**4) 228 else 229 uu=rrad(inx+1,nn)-rsmin 230 zz=rsmin-rrad(inx,nn) 231 232 rdint(ii,jj)=ncrho(inx)/2._dp/ss*uu*uu+ncrho(inx+1)/2._dp*(ss-zz*zz/ss)+& 233 nsp3(inx)/1.2d1*(uu**4/2._dp/ss-uu*uu*ss)+nsp3(inx+1)/1.2d1*(zz*zz*ss-ss**3/2._dp-zz**4/2._dp/ss) 234 if (inxf > inx+1) then 235 do kk=inx+1,inxf-1 236 uu=rrad(kk+1,nn)-rrad(kk,nn) 237 rdint(ii,jj)=rdint(ii,jj)+(ncrho(kk)+ncrho(kk+1))*uu/2._dp-uu*uu*uu/2.4d1*(nsp3(kk)+nsp3(kk+1)) 238 end do 239 end if 240 241 if (inxf/=ndat(nn)) then 242 uu=rs(ii,jj)-rrad(inxf,nn) 243 zz=rrad(inxf+1,nn)-rs(ii,jj) 244 ss=rrad(inxf+1,nn)-rrad(inxf,nn) 245 rdint(ii,jj)=rdint(ii,jj)+ncrho(inxf)/2._dp*(ss-zz*zz/ss)+ncrho(inxf+1)/2._dp*uu*uu/ss+& 246 nsp3(inxf)/1.2d1*(zz*zz*ss-zz*zz*zz*zz/2._dp/ss-ss*ss*ss/2._dp)+& 247 nsp3(inxf+1)/1.2d1*(uu*uu*uu*uu/2._dp/ss-uu*uu*ss) 248 end if 249 end if 250 rdint(ii,jj)=rdint(ii,jj)/4._dp/pi 251 nsphe=nsphe+rdint(ii,jj)*wgrs(ii,jj) 252 end do 253 end do 254 nsphe=nsphe*(pi/(themin-themax))*(two_pi/(phimax-phimin)) 255 256 write(untout,*) 257 write(untout,*) "CHARGE INTEGRATION" 258 write(untout,*) "==================" 259 write(untout,'(" Core density contribution: ",/,/," ",F16.8)') cintr+nsphe 260 261 write(std_out,*) ':INTECOR ', cintr+nsphe 262 263 ABI_DEALLOCATE(ncrho) 264 ABI_DEALLOCATE(nsp2) 265 ABI_DEALLOCATE(nsp3) 266 ABI_DEALLOCATE(nsp4) 267 268 !INTEGRATION OF THE VALENCE DENSITY 269 270 ABI_ALLOCATE(rr,(aim_dtset%npt+1)) 271 ABI_ALLOCATE(vrho,(aim_dtset%npt+1)) 272 ABI_ALLOCATE(vdd,(aim_dtset%npt+1)) 273 274 !in the case of the only irho appelation 275 276 nn=0 277 do ii=-3,3 278 do jj=-3,3 279 do kk=-3,3 280 nn=nn+1 281 atp(1,nn)=ii*1._dp 282 atp(2,nn)=jj*1._dp 283 atp(3,nn)=kk*1._dp 284 call bschg1(atp(:,nn),1) 285 if ((ii==0).and.(jj==0).and.(kk==0)) ipos=nn 286 end do 287 end do 288 end do 289 nnpos=nn 290 iat=batom 291 292 !XG020629 There is a problem with this routine 293 !(or vgh_rho), when one uses the PGI compiler : 294 !The following line is needed, otherwise, iat and ipos 295 !are set to 0 inside vgh_now. Why ???? 296 write(std_out,*)' integrho : iat,ipos=',iat,ipos 297 ! 298 299 nsphe=0._dp 300 ABI_ALLOCATE(work,(aim_dtset%npt+1)) 301 do ii=1,nth 302 do jj=1,nph 303 304 stp=rs(ii,jj)/aim_dtset%npt 305 unvec(1)=sin(th(ii))*cos(ph(jj)) 306 unvec(2)=sin(th(ii))*sin(ph(jj)) 307 unvec(3)=cos(th(ii)) 308 do kk=0,aim_dtset%npt 309 rr(kk+1)=kk*stp 310 vv(:)=xatm(:,batom)+kk*stp*unvec(:) 311 chs=-2 312 call vgh_rho(vv,chg,grho,hrho,uu,iat,ipos,chs) 313 vrho(kk+1)=chg*rr(kk+1)*rr(kk+1) 314 if (kk==aim_dtset%npt) then 315 rder=0._dp 316 do ll=1,3 317 rder=rder+grho(ll)*unvec(ll) 318 end do 319 rder=rder*rr(kk+1)*rr(kk+1)+2._dp*rr(kk+1)*chg 320 end if 321 end do 322 lder=0._dp 323 kk=aim_dtset%npt+1 324 call spline(rr,vrho,kk,lder,rder,vdd) 325 326 ! INTEGRATION 327 328 do kk=1,aim_dtset%npt 329 rdint(ii,jj)=rdint(ii,jj)+stp/2._dp*(vrho(kk)+vrho(kk+1))& 330 & -stp*stp*stp/24._dp*(vdd(kk)+vdd(kk+1)) 331 end do 332 nsphe=nsphe+rdint(ii,jj)*wgrs(ii,jj) 333 end do 334 end do 335 ABI_DEALLOCATE(work) 336 337 if (gaus.or.weit) then 338 nsphe=nsphe*(pi/(themin-themax))*(two_pi/(phimax-phimin)) 339 else 340 nsphe=nsphe/(nth*nph)*2.0*two_pi 341 end if 342 chgint=cintr+nsphe 343 344 write(untout,'(/," Different density contributions: Core (only spherical part) and the rest ",/,/," ",2F16.8)') & 345 & cintr, nsphe 346 write(untout,'(/,a,i4,a,f14.8)') ' For atom number ',batom,', the number of electrons in the Bader volume is ',chgint 347 write(untout,'(a,f15.7,a,f17.8)') ' The nuclear charge is',znucl_batom,', so that the Bader charge is ',znucl_batom-chgint 348 write(untout,*) 349 write(std_out,*) ':INTEPAR ', cintr, nsphe 350 write(std_out,*) ':RHOTOT ',batom,chgint 351 352 end subroutine integrho