1 2* 3* $Id$ 4* 5 6* *************************** 7* * * 8* * generate_ELF * 9* * * 10* *************************** 11 12* This routine returns the electron localization function (ELF) defined 13* by Becke and Edgecombe, J. Chem. Phys., vol 92, page 5397-5403, 1990. 14* 15* ELF = (1+ (D/D0)^2)^-1 16* 17* where 18* 19* D0 = (3/5)*(6*pi^2)^(2/3) * dn^(5/3) 20* D = Sum(i) |del*psi_i|^2 - 0.25*|del*dn|^2/dn 21* 22 subroutine generate_ELF(npack1,ne,psi,dn,elf) 23 implicit none 24#include "errquit.fh" 25 integer npack1,ne 26 complex*16 psi(npack1,ne) 27 real*8 dn(*) 28 real*8 elf(*) 29 30#include "bafdecls.fh" 31 32* **** local variables **** 33 logical value 34 integer i,k,npack0,nfft3d,n2ft3d 35 integer dng(2),dng2(2),dns(2),dnsp(2) 36 integer nx,ny,nz 37 real*8 tt,vv,uu,scaluu 38 real*8 scal1,scal2 39 40* **** external functions **** 41 real*8 lattice_omega 42 external lattice_omega 43 44 call D3dB_nx(1,nx) 45 call D3dB_ny(1,ny) 46 call D3dB_nz(1,nz) 47 scal1 = 1.0d0/dble(nx*ny*nz) 48 scal2 = 1.0d0/lattice_omega() 49 50 call Pack_npack(0,npack0) 51 call D3dB_nfft3d(1,nfft3d) 52 call D3dB_n2ft3d(1,n2ft3d) 53 54 call dcopy(n2ft3d,0.0d0,0,elf,1) 55 56 value = BA_push_get(mt_dbl,n2ft3d,'dns',dns(2),dns(1)) 57 value = value.and. 58 > BA_push_get(mt_dbl,n2ft3d,'dnsp',dnsp(2),dnsp(1)) 59 value = value.and. 60 > BA_push_get(mt_dcpl,nfft3d,'dng',dng(2),dng(1)) 61 value = value.and. 62 > BA_push_get(mt_dcpl,nfft3d,'dng2',dng2(2),dng2(1)) 63 if (.not.value) call errquit('generate_ELF:pushing stack',0, 64 & MA_ERR) 65 call Ursenbach_smoother(n2ft3d,dn, 66 > dbl_mb(dns(1)), 67 > dbl_mb(dnsp(1))) 68 do i=1,n2ft3d 69 dbl_mb(dns(1)+i-1) = log(dn(i)) 70 end do 71 72* **** elf <-- elf + |del*rho|^2/rho **** 73 do i=1,n2ft3d 74 dbl_mb(dns(1)+i-1) = log(dn(i)) 75 end do 76 !call D3dB_r_SMul(1,scal1,dbl_mb(dns(1)),dcpl_mb(dng(1))) 77 call D3dB_r_SMul(1,scal1,dbl_mb(int(dn(1))),dcpl_mb(dng(1))) 78 call D3dB_r_SMul(1,scal1,dbl_mb(dns(1)),dcpl_mb(dng2(1))) 79 call D3dB_r_Zero_Ends(1,dcpl_mb(dng(1))) 80 call D3dB_r_Zero_Ends(1,dcpl_mb(dng2(1))) 81 call D3dB_rc_fft3f(1,dcpl_mb(dng(1))) 82 call D3dB_rc_fft3f(1,dcpl_mb(dng2(1))) 83 call Pack_c_pack(0,dcpl_mb(dng(1))) 84 call Pack_c_pack(0,dcpl_mb(dng2(1))) 85 call gradient_energy_density2(0,npack0, 86 > dcpl_mb(dng(1)), 87 > dcpl_mb(dng2(2)), 88 > n2ft3d,elf) 89 !call D3dB_rr_Mul(1,dbl_mb(dnsp(1)),elf,elf) 90 !call D3dB_rr_Mul(1,dbl_mb(dnsp(1)),elf,elf) 91 92* **** elf <-- -0.25*elf **** 93 !call D3dB_rr_Divide(1,elf,dn,elf) 94c call D3dB_r_SMul(1,(-0.25d0),elf,elf) 95 call D3dB_r_SMul1(1,(-0.25d0),elf) 96 97* **** elf <-- elf + Sum(i) |del*psi_i|^2 **** 98 do i=1,ne 99 call dcopy(n2ft3d,0.0d0,0,dcpl_mb(dng(1)),1) 100 call Pack_c_SMul(1,dsqrt(scal2),psi(1,i),dcpl_mb(dng(1))) 101 call gradient_energy_density(1,npack1,dcpl_mb(dng(1)), 102 > n2ft3d,elf) 103 end do 104 105* **** elf <-- (1+(elf/D0)^2)^(-1) **** 106 scaluu = 4.0d0*datan(1.0d0) 107 scaluu = (3.0d0/5.0d0)*(6.0d0*scaluu**2)**(2.0d0/3.0d0) 108 do k=1,n2ft3d 109 tt = elf(k) 110 vv = dn(k) 111 vv = scaluu*(vv**(5.0d0/3.0d0)) 112 !uu = 1.0d0+(tt/vv)**2 113 !elf(k) = 1.0d0/uu 114 uu = vv*vv + tt*tt 115 uu = 1.0d0/uu 116 elf(k) = vv*vv*uu 117 end do 118 119 10 continue 120 121 value = BA_pop_stack(dng2(2)) 122 value = value.and.BA_pop_stack(dng(2)) 123 value = value.and.BA_pop_stack(dnsp(2)) 124 value = value.and.BA_pop_stack(dns(2)) 125 if (.not.value) call errquit('generate_ELF:popping stack',1, 126 & MA_ERR) 127 128 return 129 end 130 131 132* *********************************** 133* * * 134* * gradient_energy_density * 135* * * 136* *********************************** 137* 138* This routine calculates 139* 140* df = df + |del*f|^2 141* 142 subroutine gradient_energy_density(nb,npack,f, 143 > n2ft3d,df) 144 implicit none 145#include "errquit.fh" 146 integer nb,npack 147 complex*16 f(*) 148 integer n2ft3d 149 real*8 df(*) 150 151#include "bafdecls.fh" 152 153* **** local variables **** 154 logical value 155 integer nfft3d 156 integer fx(2),fy(2),fz(2) 157 integer Gx(2),Gy(2),Gz(2) 158 159* **** external functions **** 160 integer G_indx 161 external G_indx 162 163 164 nfft3d = n2ft3d/2 165 value = BA_push_get(mt_dbl,n2ft3d, 166 > 'fx',fx(2),fx(1)) 167 value = value.and.BA_push_get(mt_dbl,n2ft3d, 168 > 'fy',fy(2),fy(1)) 169 value = value.and.BA_push_get(mt_dbl,n2ft3d, 170 > 'fz',fz(2),fz(1)) 171 value = value.and.BA_push_get(mt_dbl,nfft3d, 172 > 'Gx',Gx(2),Gx(1)) 173 value = value.and.BA_push_get(mt_dbl,nfft3d, 174 > 'Gy',Gy(2),Gy(1)) 175 value = value.and.BA_push_get(mt_dbl,nfft3d, 176 > 'Gz',Gz(2),Gz(1)) 177 if (.not.value) 178 > call errquit('gradient_energy_density:push stack',0, MA_ERR) 179 180* **** define Gx,Gy,Gz in packed space **** 181 call dcopy(nfft3d,dbl_mb(G_indx(1)),1,dbl_mb(Gx(1)),1) 182 call dcopy(nfft3d,dbl_mb(G_indx(2)),1,dbl_mb(Gy(1)),1) 183 call dcopy(nfft3d,dbl_mb(G_indx(3)),1,dbl_mb(Gz(1)),1) 184 call Pack_t_pack(nb,dbl_mb(Gx(1))) 185 call Pack_t_pack(nb,dbl_mb(Gy(1))) 186 call Pack_t_pack(nb,dbl_mb(Gz(1))) 187 188* ***** calculate: **** 189* **** fx(G)=i*Gx(G)*f(G) **** 190* **** fy(G)=i*Gy(G)*f(G) **** 191* **** fz(G)=i*Gz(G)*f(G) **** 192 call Pack_itc_Mul(nb,dbl_mb(Gx(1)),f,dbl_mb(fx(1))) 193 call Pack_itc_Mul(nb,dbl_mb(Gy(1)),f,dbl_mb(fy(1))) 194 call Pack_itc_Mul(nb,dbl_mb(Gz(1)),f,dbl_mb(fz(1))) 195 196 value = BA_pop_stack(Gz(2)) 197 value = value.and.BA_pop_stack(Gy(2)) 198 value = value.and.BA_pop_stack(Gx(2)) 199 if (.not.value) 200 > call errquit('gradient_energy_density:pop stack',1, MA_ERR) 201 202* **** Fourier transform fx,fy,fz to real space **** 203 call Pack_c_unpack(nb,dbl_mb(fx(1))) 204 call Pack_c_unpack(nb,dbl_mb(fy(1))) 205 call Pack_c_unpack(nb,dbl_mb(fz(1))) 206 call D3dB_cr_fft3b(1,dbl_mb(fx(1))) 207 call D3dB_cr_fft3b(1,dbl_mb(fy(1))) 208 call D3dB_cr_fft3b(1,dbl_mb(fz(1))) 209 210* **** calculate df = df + fx*fx + fy*fy + fz*fz **** 211c call D3dB_rr_Sqr(1,dbl_mb(fx(1)),dbl_mb(fx(1))) 212c call D3dB_rr_Sqr(1,dbl_mb(fy(1)),dbl_mb(fy(1))) 213c call D3dB_rr_Sqr(1,dbl_mb(fz(1)),dbl_mb(fz(1))) 214c call D3dB_rr_Sum(1,dbl_mb(fx(1)),df,df) 215c call D3dB_rr_Sum(1,dbl_mb(fy(1)),df,df) 216c call D3dB_rr_Sum(1,dbl_mb(fz(1)),df,df) 217 call D3dB_rr_Sqr1(1,dbl_mb(fx(1))) 218 call D3dB_rr_Sqr1(1,dbl_mb(fy(1))) 219 call D3dB_rr_Sqr1(1,dbl_mb(fz(1))) 220 call D3dB_rr_Sum2(1,dbl_mb(fx(1)),df) 221 call D3dB_rr_Sum2(1,dbl_mb(fy(1)),df) 222 call D3dB_rr_Sum2(1,dbl_mb(fz(1)),df) 223 224 value = BA_pop_stack(fz(2)) 225 value = value.and.BA_pop_stack(fy(2)) 226 value = value.and.BA_pop_stack(fx(2)) 227 if (.not.value) 228 > call errquit('gradient_energy_density:pop stack',2, MA_ERR) 229 230 return 231 end 232 233* *********************************** 234* * * 235* * gradient_energy_density2 * 236* * * 237* *********************************** 238* 239* This routine calculates 240* 241* df = df + del*f * del*h 242* 243 subroutine gradient_energy_density2(nb,npack,f,h, 244 > n2ft3d,df) 245 implicit none 246#include "errquit.fh" 247 integer nb,npack 248 complex*16 f(*),h(*) 249 integer n2ft3d 250 real*8 df(*) 251 252#include "bafdecls.fh" 253 254* **** local variables **** 255 logical value 256 integer nfft3d 257 integer fx(2),fy(2),fz(2) 258 integer hx(2),hy(2),hz(2) 259 integer Gx(2),Gy(2),Gz(2) 260 261* **** external functions **** 262 integer G_indx 263 external G_indx 264 265 266 nfft3d = n2ft3d/2 267 value = BA_push_get(mt_dbl,n2ft3d, 268 > 'hx',hx(2),hx(1)) 269 value = value.and.BA_push_get(mt_dbl,n2ft3d, 270 > 'hy',hy(2),hy(1)) 271 value = value.and.BA_push_get(mt_dbl,n2ft3d, 272 > 'hz',hz(2),hz(1)) 273 value = value.and.BA_push_get(mt_dbl,n2ft3d, 274 > 'fx',fx(2),fx(1)) 275 value = value.and.BA_push_get(mt_dbl,n2ft3d, 276 > 'fy',fy(2),fy(1)) 277 value = value.and.BA_push_get(mt_dbl,n2ft3d, 278 > 'fz',fz(2),fz(1)) 279 280 value = value.and.BA_push_get(mt_dbl,nfft3d, 281 > 'Gx',Gx(2),Gx(1)) 282 value = value.and.BA_push_get(mt_dbl,nfft3d, 283 > 'Gy',Gy(2),Gy(1)) 284 value = value.and.BA_push_get(mt_dbl,nfft3d, 285 > 'Gz',Gz(2),Gz(1)) 286 if (.not.value) 287 > call errquit('gradient_energy_density:push stack',0, MA_ERR) 288 289* **** define Gx,Gy,Gz in packed space **** 290 call dcopy(nfft3d,dbl_mb(G_indx(1)),1,dbl_mb(Gx(1)),1) 291 call dcopy(nfft3d,dbl_mb(G_indx(2)),1,dbl_mb(Gy(1)),1) 292 call dcopy(nfft3d,dbl_mb(G_indx(3)),1,dbl_mb(Gz(1)),1) 293 call Pack_t_pack(nb,dbl_mb(Gx(1))) 294 call Pack_t_pack(nb,dbl_mb(Gy(1))) 295 call Pack_t_pack(nb,dbl_mb(Gz(1))) 296 297* ***** calculate: **** 298* **** fx(G)=i*Gx(G)*f(G) **** 299* **** fy(G)=i*Gy(G)*f(G) **** 300* **** fz(G)=i*Gz(G)*f(G) **** 301 call Pack_itc_Mul(nb,dbl_mb(Gx(1)),f,dbl_mb(fx(1))) 302 call Pack_itc_Mul(nb,dbl_mb(Gy(1)),f,dbl_mb(fy(1))) 303 call Pack_itc_Mul(nb,dbl_mb(Gz(1)),f,dbl_mb(fz(1))) 304 call Pack_itc_Mul(nb,dbl_mb(Gx(1)),h,dbl_mb(hx(1))) 305 call Pack_itc_Mul(nb,dbl_mb(Gy(1)),h,dbl_mb(hy(1))) 306 call Pack_itc_Mul(nb,dbl_mb(Gz(1)),h,dbl_mb(hz(1))) 307 308 value = BA_pop_stack(Gz(2)) 309 value = value.and.BA_pop_stack(Gy(2)) 310 value = value.and.BA_pop_stack(Gx(2)) 311 if (.not.value) 312 > call errquit('gradient_energy_density2:pop stack',1, MA_ERR) 313 314* **** Fourier transform fx,fy,fz to real space **** 315 call Pack_c_unpack(nb,dbl_mb(fx(1))) 316 call Pack_c_unpack(nb,dbl_mb(fy(1))) 317 call Pack_c_unpack(nb,dbl_mb(fz(1))) 318 call D3dB_cr_fft3b(1,dbl_mb(fx(1))) 319 call D3dB_cr_fft3b(1,dbl_mb(fy(1))) 320 call D3dB_cr_fft3b(1,dbl_mb(fz(1))) 321 322 call Pack_c_unpack(nb,dbl_mb(hx(1))) 323 call Pack_c_unpack(nb,dbl_mb(hy(1))) 324 call Pack_c_unpack(nb,dbl_mb(hz(1))) 325 call D3dB_cr_fft3b(1,dbl_mb(hx(1))) 326 call D3dB_cr_fft3b(1,dbl_mb(hy(1))) 327 call D3dB_cr_fft3b(1,dbl_mb(hz(1))) 328 329* **** calculate df = df + hx*fx + hy*fy + hz*fz **** 330c call D3dB_rr_Mul(1,dbl_mb(hx(1)),dbl_mb(fx(1)),dbl_mb(fx(1))) 331c call D3dB_rr_Mul(1,dbl_mb(hy(1)),dbl_mb(fy(1)),dbl_mb(fy(1))) 332c call D3dB_rr_Mul(1,dbl_mb(hz(1)),dbl_mb(fz(1)),dbl_mb(fz(1))) 333c call D3dB_rr_Sum(1,dbl_mb(fx(1)),df,df) 334c call D3dB_rr_Sum(1,dbl_mb(fy(1)),df,df) 335c call D3dB_rr_Sum(1,dbl_mb(fz(1)),df,df) 336 call D3dB_rr_Mul2(1,dbl_mb(hx(1)),dbl_mb(fx(1))) 337 call D3dB_rr_Mul2(1,dbl_mb(hy(1)),dbl_mb(fy(1))) 338 call D3dB_rr_Mul2(1,dbl_mb(hz(1)),dbl_mb(fz(1))) 339 call D3dB_rr_Sum2(1,dbl_mb(fx(1)),df) 340 call D3dB_rr_Sum2(1,dbl_mb(fy(1)),df) 341 call D3dB_rr_Sum2(1,dbl_mb(fz(1)),df) 342 343 value = BA_pop_stack(fz(2)) 344 value = value.and.BA_pop_stack(fy(2)) 345 value = value.and.BA_pop_stack(fx(2)) 346 value = value.and.BA_pop_stack(hz(2)) 347 value = value.and.BA_pop_stack(hy(2)) 348 value = value.and.BA_pop_stack(hx(2)) 349 if (.not.value) 350 > call errquit('gradient_energy_density2:pop stack',2, MA_ERR) 351 352 return 353 end 354 355 356 357 358* ********************************** 359* * * 360* * generate_dmatrix_column * 361* * * 362* ********************************** 363 364 subroutine generate_dmatrix_column(ms,xyz, 365 > ispin,ne,n2ft3d,psi_r,rho) 366 implicit none 367 integer ms 368 real*8 xyz(3) 369 integer ispin,ne(2),n2ft3d 370 real*8 psi_r(n2ft3d,ne(1)+ne(2)) 371 real*8 rho(*) 372 373#include "bafdecls.fh" 374#include "errquit.fh" 375 376* **** local variables **** 377 logical value 378 integer i,k,jj,jjindex(100),rgrid(2),nx,ny,nz 379 real*8 summ(100),tsum,msum,x,y,z,dv,dx,dy,dz,rr,rrmin 380 381* **** external functions **** 382 real*8 lattice_omega 383 external lattice_omega 384 385 if (.not.BA_push_get(mt_dbl,3*n2ft3d,'rgrd',rgrid(2),rgrid(1))) 386 > call errquit('pspw_dplot_loop:push stack',0,MA_ERR) 387 call lattice_r_grid(dbl_mb(rgrid(1))) 388 389 call D3dB_nx(1,nx) 390 call D3dB_ny(1,ny) 391 call D3dB_nz(1,nz) 392 dv = lattice_omega()/dble(nx*ny*nz) 393 394c *** find nearest i *** 395 rrmin = 9.9d99 396 jj = 0 397 do i=1,n2ft3d 398 x = dbl_mb(rgrid(1)+3*(i-1)) 399 y = dbl_mb(rgrid(1)+3*(i-1)+1) 400 z = dbl_mb(rgrid(1)+3*(i-1)+2) 401 dx = x-xyz(1) 402 dy = y-xyz(2) 403 dz = z-xyz(3) 404 rr = dx**2 + dy**2 + dz**2 405 if (rr.lt.rrmin) then 406 rrmin = rr 407 jj = i 408 end if 409 end do 410 411c jj = j 412c msum = 0.0d0 413c do i=1,n2ft3d 414c tsum = 0.0d0 415c do k=1,ne(ms) 416c tsum = tsum + psi_r(i,k+(ms-1)*ne(1))**2 417c end do 418c if (tsum.gt.msum) then 419c jj = i 420c msum = tsum 421c end if 422c end do 423 write(*,*) "jj=",jj,rrmin 424 x = dbl_mb(rgrid(1)+3*(jj-1)) 425 y = dbl_mb(rgrid(1)+3*(jj-1)+1) 426 z = dbl_mb(rgrid(1)+3*(jj-1)+2) 427 428 do k=1,ne(ms) 429 summ(k) = psi_r(jj,k+(ms-1)*ne(1)) 430 end do 431 call dcopy(n2ft3d,0.0d0,0,rho,1) 432 do i=1,n2ft3d 433 do k=1,ne(ms) 434 rho(i) = rho(i) + psi_r(i,k+(ms-1)*ne(1))*summ(k) 435 end do 436 end do 437 do i=1,n2ft3d 438 tsum = tsum + rho(i) 439 end do 440 tsum = tsum * dv 441 write(*,*) "<rho(:,j)> = ",jj,x,y,z,tsum 442 443 if (.not.BA_pop_stack(rgrid(2))) 444 > call errquit('pspw_dplot_loop:pop stack',0,MA_ERR) 445 return 446 end 447 448 449