1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 subroutine phirphi_opt( nua, na, nuo, no, scell, xa, maxnh, lasto, 9 & lastkb, iphorb, iphKB, isa, numh, 10 & listhptr, listh, dk, matrix, Sp ) 11C ********************************************************************* 12C If matrix='R' 13C Finds the matrix elements of the dk*r between basis 14C orbitals, where dk is a given vector and r is the position operator, 15C or the momentum operator plus a correction due to the non-local 16C potential in the case of an infinite solid. 17C 18C S_a,b(R_a,R_b) = <phi_a(r-R_a-t_a)|dk*r|phi_b(r-R_b-t_b)> 19C 20C 21C Being R_a and R_b lattice vectors, and t_a and t_b the coordiantes 22C of atoms a and b in the unit cell 23C 24C If matrix='P' 25C Finds the matrix elements of the -i*dk*p (i.e. minus dk dot gradient ) 26C between basis orbitals, phi_a(r-R_a) and phi_b(r-R_b), 27C where dk is a given vector and p is the momentum operator, 28C a correction must be included due to the use of non-local pseudoptentials 29C VNL=|f(r-R_c)><f(r'-R_c)| 30C 31C S_a,b(R_a,R_b) = -i*2.0d0*<phi_a(r-R_a-t_a)|dk*p|phi_b(r-R_b-t_b)>+ 32C <phi_a(r-R_a-t_a)|f(r-R_c)><f(r'-R_c)|dk*(r'-R_c)|phi_b(r'-R_b-t_b)> 33C -<phi_a(r-R_a-t_a)|dk*(r-R_c)|f(r-R_c)><f(r'-R_c)|phi_b(r'-R_b-t_b)> 34C 35C The factor of two in front of dk*p is related to the use of Ry and 36C Bohr instead of Ha and Ry. The factor will be canceled out 37C afterwards when we divide by (E_i-E-j), the energy denominator. 38C Being R_a and R_b lattice vectors, and t_a and t_b the coordiantes 39C of atoms a and b in the unit cell 40C 41 42C Energies in Ry. Lengths in Bohr. 43C Written by DSP August 1999 ( Based in routine overfsm and nlefsm) 44C Restyled for f90 version by JDG. June 2004 45C **************************** INPUT ********************************** 46C integer nua : Number of atoms in unit cell 47C integer na : Number of atoms in supercell 48C integer nuo : Number of basis orbitals in unit cell 49C integer no : Number of basis orbitals in super cell 50C real*8 scell(3,3) : Supercell vectors SCELL(IXYZ,IVECT) 51C real*8 xa(3,na) : Atomic positions in cartesian coordinates 52C integer maxnh : First dimension of listh 53C integer lasto(0:na) : Last orbital index of each atom 54C integer lastkb(0:na) : Last KB porjector index of each atom 55C integer iphorb(no) : Orbital index of each orbital in its atom 56C integer iphKB(nokb) : KB proj. index of each KB proj. in its atom 57C integer isa(na) : Species index of each atom 58C integer numh(nuo) : Number of nonzero elements of each row 59C of the overlap matrix 60C integer listh(maxnh) : Column indexes of the nonzero elements 61C of each row of the overlap matrix 62C real*8 dk(3) : Vector in k-space 63C character*1 matrix : 'R' for molecules or atoms, the matrix 64C elements of the position operator are 65C calculated 66C 'P' for solids, the matrix elements of 67C the momentum operator are calculated. 68C **************************** OUTPUT ********************************* 69C real*8 Sp(maxnh) : Sparse overlap matrix 70C ********************************************************************* 71 72 use precision, only : dp 73 use parallel, only : Node, Nodes 74 use atmparams, only : lmx2, nzetmx, nsemx 75 use atmfuncs, only : epskb, lofio, mofio, rcut, rphiatm 76 use atmfuncs, only : orb_gindex, kbproj_gindex 77 use atm_types, only : nspecies 78 use parallelsubs, only : GlobalToLocalOrb, LocalToGlobalOrb 79 use alloc, only : re_alloc, de_alloc 80 use sys, only : die 81 use neighbour, only : jna=>jan, xij, r2ij 82 use neighbour, only : mneighb, reset_neighbour_arrays 83 use m_new_matel, only : new_matel 84 85 implicit none 86 87C Passed variables 88 integer :: maxnh, na, no, nua, nuo 89 integer :: iphorb(no), isa(na), lasto(0:na), lastkb(0:na), 90 & listh(maxnh), listhptr(nuo), numh(nuo), iphKB(*) 91 real(dp) :: scell(3,3), dk(3), Sp(maxnh), xa(3,na) 92 character :: matrix*1 93 94C Internal variables 95C maxkba = maximum number of KB projectors of one atom 96C maxno = maximum number of basis orbitals overlapping a KB projector 97 98 integer 99 & ia, iio, ind, io, ioa, is, ix, 100 & j, ja, jn, jo, joa, js, nnia, norb, ka 101 102 real(dp) 103 & grSij(3), rij, Sij, xinv(3), sum 104 105 integer 106 & ikb, ina, ino, jno, ko, koa, ks, ig, jg, kg, 107 & nkb, nna, nno, ilm1, ilm2, npoints, 108 & ir 109 110 real(dp) 111 & epsk, grSki(3), 112 & rki, rmax, rmaxkb, rmaxo, 113 & Sik, Sjk, Sikr, Sjkr, 114 & dintg2(3), dint1, dint2, dintgmod2, 115 & dintg1(3), dintgmod1, 116 & phi1, phi2, dphi1dr, dphi2dr, Sir0, r 117 118 integer, pointer, save :: iano(:) 119 integer, pointer, save :: iono(:) 120 integer, save :: maxkba = 25 121 integer, save :: maxno = 1000 122!N logical, pointer, save :: calculated(:,:,:) 123 logical, pointer, save :: listed(:) 124 logical, pointer, save :: needed(:) 125 logical :: within 126 real(dp), save :: dx = 0.01d0 127!N real(dp), pointer, save :: Pij(:,:,:) 128 real(dp), pointer, save :: Si(:) 129 real(dp), pointer, save :: Ski(:,:,:) 130 real(dp), save :: tiny = 1.0d-9 131 real(dp), pointer, save :: Vi(:) 132 133C Start timer 134 call timer('phirphiopt',1) 135 136C Check input matrix 137 if(matrix.ne.'P'.and.matrix.ne.'R') 138 & call die('phirphi_opt: matrix only can take values R or P') 139 140C Nullify pointers 141 nullify( listed, needed, Si, Vi, iano, iono, Ski ) 142!N nullify( Pij, calculated ) 143 144C Allocate arrays 145 norb = lmx2*nzetmx*nsemx 146 call re_alloc( listed, 1, no, 'listed', 'phirphi_opt' ) 147 call re_alloc( needed, 1, no, 'needed', 'phirphi_opt' ) 148 call re_alloc( Si, 1, no, 'Si', 'phirphi_opt' ) 149 call re_alloc( Vi, 1, no, 'Vi', 'phirphi_opt' ) 150 call re_alloc( iano, 1, maxno, 'iano', 'phirphi_opt' ) 151 call re_alloc( iono, 1, maxno, 'iono', 'phirphi_opt' ) 152 call re_alloc( Ski, 1, 2, 1, maxkba, 1, maxno, 153 & 'Ski', 'phirphi_opt' ) 154!N call re_alloc( Pij, 1, norb, 1, norb, 1, nspecies, 155!N & 'Pij', 'phirphi_opt' ) 156!N call re_alloc( calculated, 1, norb, 1, norb, 1, nspecies, 157!N & 'calculated', 'phirphi_opt' ) 158 159C Initialise quantities 160 do io = 1,maxnh 161 Sp(io) = 0.0d0 162 enddo 163!N do ka = 1,nspecies 164!N do io = 1,norb 165!N do jo = 1,norb 166!N calculated(jo,io,ka) = .false. 167!N enddo 168!N enddo 169!N enddo 170 171C Find maximum range 172 rmaxo = 0.0d0 173 rmaxkb = 0.0d0 174 do ia = 1,na 175 is = isa(ia) 176 do ikb = lastkb(ia-1)+1,lastkb(ia) 177 ioa = iphKB(ikb) 178 rmaxkb = max( rmaxkb, rcut(is,ioa) ) 179 enddo 180 do io = lasto(ia-1)+1,lasto(ia) 181 ioa = iphorb(io) 182 rmaxo = max( rmaxo, rcut(is,ioa) ) 183 enddo 184 enddo 185 rmax = rmaxo + rmaxkb 186 187C Correction due to the non-locality of the potential have to be 188C calculated for the momentum operator matrix elements 189 190 if (matrix.eq.'P') then 191 192C Initialize arrayd Vi only once 193 no = lasto(na) 194 do jo = 1,no 195 Vi(jo) = 0.0d0 196 listed(jo) = .false. 197 needed(jo) = .false. 198 enddo 199 200C Find out which orbitals are needed on this node 201 do iio = 1,nuo 202 call LocalToGlobalOrb(iio,Node,Nodes,io) 203 needed(io) = .true. 204 do j = 1,numh(iio) 205 ind = listhptr(iio) + j 206 jo = listh(ind) 207 needed(jo) = .true. 208 enddo 209 enddo 210 211C Loop on atoms with KB projectors 212 do ka = 1,na 213 ks = isa(ka) 214 nkb = lastkb(ka) - lastkb(ka-1) 215 if (nkb.gt.maxkba) then 216 maxkba = nkb + 10 217 call re_alloc( Ski, 1, 2, 1, maxkba, 1, maxno, copy=.true., 218 & name='Ski', routine='phirphi_opt' ) 219 endif 220 221C Find neighbour atoms 222 call mneighb( scell, rmax, na, xa, ka, 0, nna ) 223 224C Find neighbour orbitals 225 nno = 0 226 do ina = 1,nna 227 ia = jna(ina) 228 229 is = isa(ia) 230 rki = sqrt(r2ij(ina)) 231 do io = lasto(ia-1)+1,lasto(ia) 232 if (needed(io)) then 233 call GlobalToLocalOrb(io,Node,Nodes,iio) 234 ioa = iphorb(io) 235 ig = orb_gindex(is,ioa) 236 237C Find if orbital is within range 238 within = .false. 239 do ko = lastkb(ka-1)+1,lastkb(ka) 240 koa = iphKB(ko) 241 if ( rki .lt. rcut(is,ioa)+rcut(ks,koa) ) 242 & within = .true. 243 enddo 244 245C Find overlap between neighbour orbitals and KB projectors 246 if (within) then 247 nno = nno + 1 248 if (nno.gt.maxno) then 249 maxno = nno + 500 250 call re_alloc( iano, 1, maxno, copy=.true., 251 & name='iano', routine='phirphi_opt' ) 252 call re_alloc( iono, 1, maxno, copy=.true., 253 & name='iono', routine='phirphi_opt' ) 254 call re_alloc( Ski, 1, 2, 1, maxkba, 1, maxno, 255 & copy=.true., name='Ski', 256 & routine='phirphi_opt' ) 257 endif 258 iono(nno) = io 259 iano(nno) = ia 260 ikb = 0 261 do ko = lastkb(ka-1)+1,lastkb(ka) 262 ikb = ikb + 1 263 koa = iphKB(ko) 264 kg = kbproj_gindex(ks,koa) 265 do ix = 1,3 266 xinv(ix) = - xij(ix,ina) 267 enddo 268 call new_MATEL('S', ig, kg, xinv, 269 & Ski(1,ikb,nno), grSki) 270 271 sum = 0.0d0 272 if (abs(dk(1)).gt.tiny) then 273 call new_MATEL('X', ig, kg, xinv, 274 & Sik, grSki) 275 sum = sum + Sik*dk(1) 276 endif 277 if (abs(dk(2)).gt.tiny) then 278 call new_MATEL('Y', ig, kg, xinv, 279 & Sik, grSki) 280 sum = sum + Sik*dk(2) 281 endif 282 if (abs(dk(3)).gt.tiny) then 283 call new_MATEL('Z', ig, kg, xinv, 284 & Sik, grSki) 285 sum = sum + Sik*dk(3) 286 endif 287 Ski(2,ikb,nno) = sum 288 enddo 289 endif 290 endif 291 enddo 292 293 enddo 294 295C Loop on neighbour orbitals 296 do ino = 1,nno 297 io = iono(ino) 298 call GlobalToLocalOrb(io,Node,Nodes,iio) 299 if (iio .gt. 0) then 300 ia = iano(ino) 301 if (ia .le. nua) then 302 303C Scatter filter of desired matrix elements 304 do j = 1,numh(iio) 305 ind = listhptr(iio) + j 306 jo = listh(ind) 307 listed(jo) = .true. 308 enddo 309 310C Find matrix elements with other neighbour orbitals 311 do jno = 1,nno 312 jo = iono(jno) 313 if (listed(jo)) then 314 315C Loop on KB projectors 316 ikb = 0 317 do ko = lastkb(ka-1)+1,lastkb(ka) 318 ikb = ikb + 1 319 koa = iphKB(ko) 320 epsk = epskb(ks,koa) 321 Sik = Ski(1,ikb,ino) 322 Sikr= Ski(2,ikb,ino) 323 Sjk = Ski(1,ikb,jno) 324 Sjkr= Ski(2,ikb,jno) 325 Vi(jo) = Vi(jo) + epsk * (Sik * Sjkr - Sikr * Sjk) 326 enddo 327 328 endif 329 enddo 330 331C Pick up contributions to H and restore Di and Vi 332 do j = 1,numh(iio) 333 ind = listhptr(iio) + j 334 jo = listh(ind) 335 Sp(ind) = Sp(ind) + Vi(jo) 336 Vi(jo) = 0.0d0 337 listed(jo) = .false. 338 enddo 339 340 endif 341 endif 342 enddo 343 enddo 344 345 endif 346 347C Initialize neighb subroutine 348 call mneighb( scell, 2.0d0*rmaxo, na, xa, 0, 0, nnia ) 349 350 do jo = 1,no 351 Si(jo) = 0.0d0 352 enddo 353 354 do ia = 1,nua 355 356 is = isa(ia) 357 call mneighb( scell, 2.0d0*rmaxo, na, xa, ia, 0, nnia ) 358 359 do io = lasto(ia-1)+1,lasto(ia) 360 call GlobalToLocalOrb(io,Node,Nodes,iio) 361 if (iio .gt. 0) then 362 ioa = iphorb(io) 363 ig = orb_gindex(is,ioa) 364 do jn = 1,nnia 365 do ix = 1,3 366 xinv(ix) = - xij(ix,jn) 367 enddo 368 ja = jna(jn) 369 rij = sqrt( r2ij(jn) ) 370 do jo = lasto(ja-1)+1,lasto(ja) 371 joa = iphorb(jo) 372 js = isa(ja) 373 jg = orb_gindex(js,joa) 374 375 if (rcut(is,ioa)+rcut(js,joa) .gt. rij) then 376 377 if (matrix.eq.'R') then 378 379 call new_MATEL('X', jg, ig, xinv, 380 & Sij, grSij ) 381 Si(jo) = Sij*dk(1) 382 383 call new_MATEL('Y', jg, ig, xinv, 384 & Sij, grSij ) 385 Si(jo) = Si(jo) + Sij*dk(2) 386 387 call new_MATEL('Z', jg, ig, xinv, 388 & Sij, grSij ) 389 Si(jo) = Si(jo) + Sij*dk(3) 390 391 call new_MATEL('S', ig, jg, xij(1:3,jn), 392 & Sij, grSij ) 393 Si(jo) = Si(jo) + Sij*( 394 & xa(1,ia)*dk(1) 395 & + xa(2,ia)*dk(2) 396 & + xa(3,ia)*dk(3)) 397 398 else 399 400 if (rij.lt.tiny) then 401C Perform the direct computation of the matrix element of the momentum 402C within the same atom 403!N if ( .not.calculated(joa,ioa,is) ) then 404 ilm1 = lofio(is,ioa)**2 + lofio(is,ioa) + 405 & mofio(is,ioa) + 1 406 ilm2 = lofio(is,joa)**2 + lofio(is,joa) + 407 & mofio(is,joa) + 1 408 call intgry(ilm1,ilm2,dintg2) 409 call intyyr(ilm1,ilm2,dintg1) 410 dintgmod1 = dintg1(1)**2 + dintg1(2)**2 + 411 & dintg1(3)**2 412 dintgmod2 = dintg2(1)**2 + dintg2(2)**2 + 413 & dintg2(3)**2 414 Sir0 = 0.0d0 415 if ((dintgmod2.gt.tiny).or.(dintgmod1.gt.tiny)) 416 & then 417 dint1 = 0.0d0 418 dint2 = 0.0d0 419 npoints = int(max(rcut(is,ioa),rcut(is,joa)) 420 & /dx) + 2 421 do ir = 1,npoints 422 r = dx*(ir-1) 423 call rphiatm(is,ioa,r,phi1,dphi1dr) 424 call rphiatm(is,joa,r,phi2,dphi2dr) 425 dint1 = dint1 + dx*phi1*dphi2dr*r**2 426 dint2 = dint2 + dx*phi1*phi2*r 427 enddo 428C The factor of two because we use Ry for the Hamiltonian 429 Sir0 = 430 & -2.0d0*(dk(1)*(dint1*dintg1(1)+dint2*dintg2(1))+ 431 & dk(2)*(dint1*dintg1(2)+dint2*dintg2(2))+ 432 & dk(3)*(dint1*dintg1(3)+dint2*dintg2(3))) 433 endif 434!N Pij(ioa,joa,is) = - Sir0 435!N Pij(joa,ioa,is) = Sir0 436 Si(jo) = Sir0 437!N calculated(ioa,joa,is) = .true. 438!N calculated(joa,ioa,is) = .true. 439!N endif 440 441 else 442C Matrix elements between different atoms are taken from the 443C gradient of the overlap 444 call new_MATEL('S', ig, jg, xij(1:3,jn), 445 & Sij, grSij ) 446C The factor of two because we use Ry for the Hamiltonian 447 Si(jo) = 448 & 2.0d0*(grSij(1)*dk(1) 449 & + grSij(2)*dk(2) 450 & + grSij(3)*dk(3)) 451 452 endif 453 454 endif 455 endif 456 enddo 457 enddo 458 do j = 1,numh(iio) 459 ind = listhptr(iio) + j 460 jo = listh(ind) 461 Sp(ind) = Sp(ind) + Si(jo) 462 Si(jo) = 0.0d0 463 enddo 464 endif 465 enddo 466 enddo 467 468C Free local memory 469! call new_MATEL( 'S', 0, 0, 0, 0, xij, Sij, grSij ) 470 call reset_neighbour_arrays( ) 471!N call de_alloc( calculated, 'calculated', 'phirphi_opt' ) 472!N call de_alloc( Pij, 'Pij', 'phirphi_opt' ) 473 call de_alloc( Ski, 'Ski', 'phirphi_opt' ) 474 call de_alloc( iono, 'iono', 'phirphi_opt' ) 475 call de_alloc( iano, 'iano', 'phirphi_opt' ) 476 call de_alloc( Vi, 'Vi', 'phirphi_opt' ) 477 call de_alloc( Si, 'Si', 'phirphi_opt' ) 478 call de_alloc( needed, 'needed', 'phirphi_opt' ) 479 call de_alloc( listed, 'listed', 'phirphi_opt' ) 480 481C Stop timer 482 call timer('phirphiopt',2) 483 484 return 485 end 486 487 488 subroutine Intgry( ilm1, ilm2, dintg ) 489C ******************************************************************* 490C Returns a vector with the value of the integral of a spherical 491C harmonic with the gradient of another spherical harmonic. 492C written by DSP. August 1999 (from subroutine ylmexp by J. Soler) 493C ************************* INPUT *********************************** 494C integer ilm1 : first spherical harmonic index: 495C ilm1=L1*L1+L1+M1+1. 496C integer ilm2 : second spherical harmonic index 497C ************************* OUTPUT ********************************** 498C real*8 dintg(3) : Vector with the value of the 499C (angular) integral of the product 500C Yl1m1*grad(Yl2m2) 501C ******************************************************************* 502C external rlylm(lmax,rvec,y,grady) : Returns spherical harmonics, 503C times r**l 504C ******************************************************************* 505 506 use precision, only: dp 507 use alloc 508 use spher_harm, only : gauleg, rlylm 509 510 implicit none 511 512C Passed variables 513 integer ilm1, ilm2 514 real(dp) dintg(3) 515 516C Internal variables 517 integer 518 & lmax, l2, ix, isp, iz, nsp, im, 519 & maxlm, maxsp 520 real(dp) 521 & phi, pi, theta, dl, r(3) 522 523 integer, save :: maxl = 8 524 logical, save :: frstme = .true. 525 real(dp), pointer, save :: gry(:,:) 526 real(dp), pointer, save :: w(:) 527 real(dp), pointer, save :: wsp(:) 528 real(dp), pointer, save :: x(:,:) 529 real(dp), pointer, save :: y(:) 530 real(dp), pointer, save :: z(:) 531 532 external 533 & dot 534 535C Nullify pointers and initialise arrays on first call 536 if (frstme) then 537 nullify(gry,w,wsp,x,y,z) 538 maxlm = (maxl+1)*(maxl+1) 539 maxsp = (maxl+1)*(2*maxl+1) 540 call re_alloc(gry,1,3,0,maxlm,name='gry',routine='intgry') 541 call re_alloc(w,1,maxl+1,name='w',routine='intgry') 542 call re_alloc(wsp,1,maxsp,name='wsp',routine='intgry') 543 call re_alloc(x,1,3,1,maxsp,name='x',routine='intgry') 544 call re_alloc(y,0,maxlm,name='y',routine='intgry') 545 call re_alloc(z,1,maxl+1,name='z',routine='intgry') 546 frstme = .false. 547 endif 548 549C Find special points and weights for gaussian quadrature 550 if (max(ilm1,ilm2).eq.1) then 551 lmax = 1 552 else 553 dl = dble(max(ilm1,ilm2)) - 1.0d0 + 1.0d-2 554 lmax = int(dsqrt(dl)) + 1 555 endif 556 557C Check array dimensions 558 if (lmax.gt.maxl) then 559 maxl = lmax 560 maxlm = (maxl+1)*(maxl+1) 561 maxsp = (maxl+1)*(2*maxl+1) 562 call re_alloc(gry,1,3,0,maxlm,name='gry',routine='intgry') 563 call re_alloc(w,1,maxl+1,name='w',routine='intgry') 564 call re_alloc(wsp,1,maxsp,name='wsp',routine='intgry') 565 call re_alloc(x,1,3,1,maxsp,name='x',routine='intgry') 566 call re_alloc(y,0,maxlm,name='y',routine='intgry') 567 call re_alloc(z,1,maxl+1,name='z',routine='intgry') 568 endif 569C 570 call gauleg( -1.0d0, 1.0d0, z, w, lmax+1 ) 571 pi = 4.0d0*atan(1.0d0) 572 nsp = 0 573 do iz = 1,lmax+1 574 theta = acos( z(iz) ) 575 do im = 0,2*lmax 576 nsp = nsp + 1 577 phi = im*2.0d0*pi/(2*lmax+1) 578 x(1,nsp) = sin(theta)*cos(phi) 579 x(2,nsp) = sin(theta)*sin(phi) 580 x(3,nsp) = cos(theta) 581 wsp(nsp) = w(iz)*(2.0d0*pi) / (2*lmax+1) 582 enddo 583 enddo 584 585C Find the integrals 586 if (ilm2.eq.1) then 587 l2 = 0 588 else 589 dl = dble(ilm2) - 1.0d0 + 1.0d-2 590 l2 = int(dsqrt(dl)) 591 endif 592 do ix = 1,3 593 dintg(ix) = 0.0d0 594 enddo 595 if (l2.eq.0) return 596 do isp = 1,nsp 597 r(1) = x(1,isp) 598 r(2) = x(2,isp) 599 r(3) = x(3,isp) 600 call rlylm( lmax, r, y, gry ) 601 do ix = 1,3 602 dintg(ix) = dintg(ix) + wsp(isp)*y(ilm1-1)* 603 & (gry(ix,ilm2-1) - l2*y(ilm2-1)*x(ix,isp)) 604 enddo 605 enddo 606 607 end 608 609 610 subroutine Intyyr( ilm1, ilm2, dintg) 611C ******************************************************************* 612C Returns a vector with the value of the integral of the product of spherical 613C two spherical harmonics and the radial unit vector. 614C ************************* INPUT *********************************** 615C integer ilm1 : first spherical harmonic index: 616C ilm1=L1*L1+L1+M1+1. 617C integer ilm2 : second spherical harmonic index 618C ************************* OUTPUT ********************************** 619C real*8 dintg(3) : Vector with the value of the 620C (angular) integral of the product 621C Yl1m1*Yl2m2*(x/r,y/r,z/r) 622C ******************************************************************* 623C external rlylm(lmax,rvec,y,grady) : Returns spherical harmonics, 624C times r**l 625C ******************************************************************* 626 627 use precision, only: dp 628 use alloc 629 use spher_harm, only : gauleg, rlylm 630 631 implicit none 632 633C Passed variables 634 integer ilm1, ilm2 635 real(dp) dintg(3) 636 637C Internal variables 638 integer 639 & lmax,ix, isp, iz, nsp, im, maxlm, maxsp 640 real(dp) 641 & phi, pi, theta, dl, r(3) 642 643 integer, save :: maxl = 8 644 logical, save :: frstme = .true. 645 real(dp), pointer, save :: gry(:,:) 646 real(dp), pointer, save :: w(:) 647 real(dp), pointer, save :: wsp(:) 648 real(dp), pointer, save :: x(:,:) 649 real(dp), pointer, save :: y(:) 650 real(dp), pointer, save :: z(:) 651 652 external 653 & dot 654 655C Nullify pointers and initialise arrays on first call 656 if (frstme) then 657 nullify(gry,w,wsp,x,y,z) 658 maxlm = (maxl+1)*(maxl+1) 659 maxsp = (maxl+1)*(2*maxl+1) 660 call re_alloc(gry,1,3,0,maxlm,name='gry',routine='intyyr') 661 call re_alloc(w,1,maxl+1,name='w',routine='intyyr') 662 call re_alloc(wsp,1,maxsp,name='wsp',routine='intyyr') 663 call re_alloc(x,1,3,1,maxsp,name='x',routine='intyyr') 664 call re_alloc(y,0,maxlm,name='y',routine='intyyr') 665 call re_alloc(z,1,maxl+1,name='z',routine='intyyr') 666 frstme = .false. 667 endif 668 669C Find special points and weights for Gaussian quadrature 670 if (max(ilm1,ilm2).eq.1) then 671 lmax = 1 672 else 673 dl = dble(max(ilm1,ilm2)) - 1.0d0 + 1.0d-2 674 lmax = int(dsqrt(dl)) + 1 675 endif 676 677C Check array dimensions 678 if (lmax.gt.maxl) then 679 maxl = lmax 680 maxlm = (maxl+1)*(maxl+1) 681 maxsp = (maxl+1)*(2*maxl+1) 682 call re_alloc(gry,1,3,0,maxlm,name='gry',routine='intgry') 683 call re_alloc(w,1,maxl+1,name='w',routine='intgry') 684 call re_alloc(wsp,1,maxsp,name='wsp',routine='intgry') 685 call re_alloc(x,1,3,1,maxsp,name='x',routine='intgry') 686 call re_alloc(y,0,maxlm,name='y',routine='intgry') 687 call re_alloc(z,1,maxl+1,name='z',routine='intgry') 688 endif 689C 690 call gauleg( -1.0d0, 1.0d0, z, w, lmax+1 ) 691 pi = 4.0d0*atan(1.0d0) 692 nsp = 0 693 do iz = 1,lmax+1 694 theta = acos(z(iz)) 695 do im = 0,2*lmax 696 nsp = nsp + 1 697 phi = im * 2.0d0*pi/(2*lmax+1) 698 x(1,nsp) = sin(theta)*cos(phi) 699 x(2,nsp) = sin(theta)*sin(phi) 700 x(3,nsp) = cos(theta) 701 wsp(nsp) = w(iz)*(2.0d0*pi)/(2*lmax+1) 702 enddo 703 enddo 704 705C Find the integrals 706 do ix = 1,3 707 dintg(ix) = 0.0d0 708 enddo 709 do isp = 1,nsp 710 r(1) = x(1,isp) 711 r(2) = x(2,isp) 712 r(3) = x(3,isp) 713 call rlylm( lmax, r, y, gry ) 714 do ix = 1,3 715 dintg(ix) = dintg(ix) + 716 & wsp(isp)*y(ilm1-1)*y(ilm2-1)*x(ix,isp) 717 enddo 718 enddo 719 720 end 721