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 module atmfuncs 9 10C This file contains a set of routines which provide all the information 11C about the basis set, pseudopotential, atomic mass, etc... of all the 12C chemical species present in the calculation. 13 14C The routines contained in this file can only be called after they 15C are initialized by calling the subroutine 'atom' for all the 16C different chemical species in the calculation: 17 18 use precision 19 use sys, only: die 20 use atm_types 21 use radial, only: rad_get, rad_func 22 use spher_harm, only: rlylm 23 24 implicit none 25! 26 character(len=79) message 27 integer, parameter :: max_l = 5 28 integer, parameter :: max_ilm = (max_l+1)*(max_l+1) 29 30 real(dp), parameter :: tiny20=1.e-20_dp 31 real(dp), parameter :: tiny12=1.e-12_dp 32 33 private :: chk, max_l, max_ilm, message 34 35 public :: nofis, nkbfis, izofis, massfis 36 public :: rcore, rchlocal, rcut, chcore_sub, epskb, uion 37 public :: atmpopfio, psch, zvalfis, floating, psover 38 public :: lofio, symfio, cnfigfio, zetafio, mofio 39 public :: labelfis, lomaxfis, nztfl, rphiatm, lmxkbfis 40 public :: phiatm, all_phi 41 public :: pol ! Added JMS Dec.2009 42 43 public :: orb_gindex, kbproj_gindex, vna_gindex, dftu_gindex 44 private 45 46 contains 47 48 subroutine chk(name,is) 49 character(len=*), intent(in) :: name 50 51 integer, intent(in) :: is 52 53 if ((is.lt.1).or.(is.gt.nspecies)) then 54 write(message,'(2a,i3,a,i3)') 55 $ name, ": Wrong species", is, ". Have", nspecies 56 call die(message) 57 endif 58 end subroutine chk 59! 60! 61 function floating(is) 62! Returns .true. if the species is really a "fake" one, intended 63! to provide some floating orbitals. 64 logical floating 65 integer, intent(in) :: is 66 67 floating = izofis(is) .lt. 0 68 end function floating 69 70 FUNCTION IZOFIS( IS ) 71 integer :: izofis ! Atomic number 72 integer, intent(in) :: is ! Species index 73 74 call chk('izofis',is) 75 76 izofis = species(is)%z 77 78 end function izofis 79 80 FUNCTION ZVALFIS( IS ) 81 real(dp) :: zvalfis ! Valence charge 82 integer, intent(in) :: is ! Species index 83 84 call chk('zvalfis',is) 85 86 zvalfis= species(is)%zval 87 end function zvalfis 88! 89 FUNCTION LABELFIS (IS) 90 character(len=20) :: labelfis ! Atomic label 91 integer, intent(in) :: is ! Species index 92 93 call chk('labelfis',is) 94 labelfis= species(is)%label 95 end function labelfis 96 97 FUNCTION LMXKBFIS (IS) 98 integer :: lmxkbfis ! Maximum ang mom of the KB projectors 99 integer, intent(in) :: is ! Species index 100 101 call chk('lmxkbfis',is) 102 lmxkbfis= species(is)%lmax_projs 103 end function lmxkbfis 104! 105 FUNCTION LOMAXFIS (IS) 106 integer :: lomaxfis ! Maximum ang mom of the Basis Functions 107 integer, intent(in) :: is ! Species index 108 109 call chk('lomaxfis',is) 110 111 lomaxfis = species(is)%lmax_basis 112 end function lomaxfis 113! 114 FUNCTION MASSFIS(IS) 115 real(dp) :: massfis ! Mass 116 integer, intent(in) :: is ! Species index 117 118 call chk('massfis',is) 119 massfis=species(is)%mass 120 end function massfis 121! 122 FUNCTION NKBFIS(IS) 123 integer :: nkbfis ! Total number of KB projectors 124 integer, intent(in) :: is ! Species index 125 126 call chk('nkbfis',is) 127 nkbfis = species(is)%nprojs 128 end function nkbfis 129! 130 131 FUNCTION NOFIS(IS) 132 integer :: nofis ! Total number of Basis functions 133 integer, intent(in) :: is ! Species index 134 135 call chk('nofis',is) 136 nofis = species(is)%norbs 137 end function nofis 138 139 FUNCTION UION ( IS ) 140 real(dp) uion 141 integer, intent(in) :: is ! Species index 142 call chk('uion',is) 143 uion = species(is)%self_energy 144 end function uion 145 146 FUNCTION RCORE(is) 147 real(dp) rcore 148 integer, intent(in) :: is ! Species index 149 150C Returns cutoff radius of the pseudo-core charge density for the non-linear 151C core corrections for xc potential. 152C Distances in Bohr 153 154 call chk('rcore',is) 155 rcore = species(is)%core%cutoff 156 157 end function rcore 158 159 FUNCTION RCHLOCAL(is) 160 real(dp) rchlocal 161 integer, intent(in) :: is ! Species index 162 163C Returns cutoff radius of the Vlocal charge density 164C Distances in Bohr 165 166 call chk('rchlocal',is) 167 rchlocal = species(is)%Chlocal%cutoff 168 169 end function rchlocal 170 171! AMENOFIS 172! 173!---- Global index helpers------------------------------ 174 175 FUNCTION orb_gindex (IS,IO) 176 integer orb_gindex 177 integer, intent(in) :: is ! Species index 178 integer, intent(in) :: io ! Orbital index (within atom) 179 180C Returns the global index of a basis orbital 181 182 call chk('orb_gindex',is) 183 if ( (io .gt. species(is)%norbs) .or. 184 $ (io .lt. 1)) call die("orb_gindex: Wrong io") 185 186 orb_gindex = species(is)%orb_gindex(io) 187 end function orb_gindex 188 189 FUNCTION kbproj_gindex (IS,IO) 190 integer kbproj_gindex 191 integer, intent(in) :: is ! Species index 192 integer, intent(in) :: io ! KBproj index 193 ! (within atom, <0, for compatibility) 194 195C Returns the global index of a KB projector 196 197 integer :: ko 198 199 call chk('kbproj_gindex',is) 200 ko = -io 201 202 if ( (ko .gt. species(is)%nprojs) .or. 203 $ (ko .lt. 1)) then 204 call die("kbproj_gindex: Wrong io") 205 endif 206 207 kbproj_gindex = species(is)%pj_gindex(ko) 208 end function kbproj_gindex 209 210 FUNCTION dftu_gindex (IS,IO) 211 integer dftu_gindex 212 integer, intent(in) :: is ! Species index 213 integer, intent(in) :: io ! Orbital index (within atom) 214 215C Returns the global index of a DFT+U projector 216 217 call chk('dftu_gindex',is) 218 if ( (io .gt. species(is)%nprojsdftu) .or. 219 $ (io .lt. 1)) call die("dftu_gindex: Wrong io") 220 221 dftu_gindex = species(is)%pjdftu_gindex(io) 222 end function dftu_gindex 223 224 FUNCTION vna_gindex (IS) 225 integer vna_gindex 226 integer, intent(in) :: is ! Species index 227 228C Returns the global index for a Vna function 229 230 call chk('vna_gindex',is) 231 vna_gindex = species(is)%vna_gindex 232 end function vna_gindex 233!------------------------------------------------ 234 235 FUNCTION ATMPOPFIO (IS,IO) 236 real(dp) atmpopfio 237 integer, intent(in) :: is ! Species index 238 integer, intent(in) :: io ! Orbital index (within atom) 239 240C Returns the population of the atomic basis orbitals in the atomic 241C ground state configuration. 242 243 call chk('atmpopfio',is) 244 if ( (io .gt. species(is)%norbs) .or. 245 $ (io .lt. 1)) call die("atmpopfio: Wrong io") 246 247 atmpopfio = species(is)%orb_pop(io) 248 end function atmpopfio 249 250 FUNCTION CNFIGFIO(IS,IO) 251 integer cnfigfio 252 integer, intent(in) :: is ! Species index 253 integer, intent(in) :: io ! Orbital index (within atom) 254 255C Returns the valence-shell configuration in the atomic ground state 256C (i.e. the principal quatum number for orbitals of angular momentum l) 257 258C INTEGER CNFIGFIO: Principal quantum number of the shell to what 259C the orbital belongs ( for polarization orbitals 260C the quantum number corresponds to the shell which 261C is polarized by the orbital io) 262 263 264 call chk('cnfigfio',is) 265 if ( (io .gt. species(is)%norbs) .or. 266 $ (io .lt. 1)) call die("cnfigfio: Wrong io") 267 268 cnfigfio = species(is)%orb_n(io) 269 270 end function cnfigfio 271 272 FUNCTION LOFIO (IS,IO) 273 integer lofio 274 integer, intent(in) :: is ! Species index 275 integer, intent(in) :: io ! Orbital index (within atom) 276 277C Returns total angular momentum quantum number of a given atomic basis 278C basis orbital or Kleynman-Bylander projector. 279 280C INTEGER IO : Orbital index (within atom) 281C IO > 0 => Basis orbitals 282C IO < 0 => Kleynman-Bylander projectors 283C IO = 0 => Local pseudopotential 284C************************OUTPUT***************************************** 285C INTEGER LOFIO : Quantum number L of orbital or KB projector 286 type(species_info), pointer :: spp 287 288 call chk('lofio',is) 289 290 spp => species(is) 291 if (io.gt.0) then 292 if (io.gt.spp%norbs) call die("lofio: No such orbital") 293 lofio = spp%orb_l(io) 294 else if (io.lt.0) then 295 if (-io.gt.spp%nprojs) call die("lofio: No such projector") 296 lofio = spp%pj_l(-io) 297 else 298 lofio = 0 299 endif 300 301 end function lofio 302 303 FUNCTION MOFIO (IS,IO) 304 integer mofio 305 integer, intent(in) :: is ! Species index 306 integer, intent(in) :: io ! Orbital index (within atom) 307 308C Returns m quantum number of a given atomic basis 309C basis orbital or Kleynman-Bylander projector. 310 311C INTEGER IO : Orbital index (within atom) 312C IO > 0 => Basis orbitals 313C IO < 0 => Kleynman-Bylander projectors 314C IO = 0 => Local pseudopotential 315C************************OUTPUT***************************************** 316C INTEGER MOFIO : Quantum number m of orbital or KB projector 317 type(species_info), pointer :: spp 318 319 call chk('mofio',is) 320 321 spp => species(is) 322 if (io.gt.0) then 323 if (io.gt.spp%norbs) call die("mofio: No such orbital") 324 mofio = spp%orb_m(io) 325 else if (io.lt.0) then 326 if (-io.gt.spp%nprojs) call die("mofio: No such projector") 327 mofio = spp%pj_m(-io) 328 else 329 mofio = 0 330 endif 331 332 end function mofio 333 334 FUNCTION ZETAFIO (IS,IO) 335 integer zetafio 336 integer, intent(in) :: is ! Species index 337 integer, intent(in) :: io ! Orbital index (within atom) 338 339C Returns zeta number of a 340C basis orbital 341 342C INTEGER IO : Orbital index (within atom) 343C IO > 0 => Basis orbitals 344C************************OUTPUT***************************************** 345C INTEGER ZETAFIO : Zeta number of orbital 346 type(species_info), pointer :: spp 347 348 call chk('mofio',is) 349 350 spp => species(is) 351 if (io.gt.0) then 352 if (io.gt.spp%norbs) call die("zetafio: No such orbital") 353 zetafio = spp%orbnl_z(spp%orb_index(io)) 354 else 355 call die('zetafio only deals with orbitals') 356 endif 357 358 end function zetafio 359 360 function rcut(is,io) 361 real(dp) rcut 362 integer, intent(in) :: is ! Species index 363 integer, intent(in) :: io ! Orbital index (within atom) 364 ! io> => basis orbitals 365 ! io<0 => KB projectors 366 ! io=0 : Local screened pseudopotential 367 368C Returns cutoff radius of Kleynman-Bylander projectors and 369C atomic basis orbitals. 370C Distances in Bohr 371 type(species_info), pointer :: spp 372 373 call chk('rcut',is) 374 375 spp => species(is) 376 if (io.gt.0) then 377 if (io.gt.spp%norbs) call die("rcut: No such orbital") 378 379 rcut = spp%orbnl(spp%orb_index(io))%cutoff 380 else if (io.lt.0) then 381 if (-io.gt.spp%nprojs) call die("rcut: No such projector") 382 rcut = spp%pjnl(spp%pj_index(-io))%cutoff 383 else 384 rcut = spp%vna%cutoff 385 endif 386 387 end function rcut 388! 389 390! 391 FUNCTION SYMFIO (IS,IO) 392 character(len=20) symfio 393 integer, intent(in) :: is ! Species index 394 integer, intent(in) :: io ! Orbital index (within atom) 395 396C Returns a label describing the symmetry of the 397C basis orbital or Kleynman-Bylander projector. 398C INTEGER IO : Orbital index (within atom) 399C IO > 0 => Basis orbitals 400C IO < 0 => Kleynman-Bylander projectors 401 402C INTEGER SYMFIO : Symmetry of the orbital or KB projector 403C 2) Returns 's' for IO = 0 404 405 integer ilm, i, lorb, morb 406 integer, parameter :: lmax_sym=4 407 408 character(len=11) sym_label((lmax_sym+1)*(lmax_sym+1)) 409 410 data sym_label(1) 411 . / 's' / 412 data (sym_label(i),i=2,4) 413 . / 'py', 'pz', 'px' / 414 data (sym_label(i),i=5,9) 415 . / 'dxy', 'dyz', 'dz2', 'dxz', 'dx2-y2' / 416 data (sym_label(i),i=10,16) 417 . / 'fy(3x2-y2)', 'fxyz', 'fz2y', 'fz3', 418 . 'fz2x', 'fz(x2-y2)', 'fx(x2-3y2)' / 419 data (sym_label(i),i=17,25) 420 . / 'gxy(x2-y2)', 'gzy(3x2-y2)', 'gz2xy', 'gz3y', 'gz4', 421 . 'gz3x', 'gz2(x2-y2)', 'gzx(x2-3y2)', 'gx4+y4' / 422 type(species_info), pointer :: spp 423 424 call chk('rcut',is) 425 426 spp => species(is) 427 if (io.gt.0) then 428 if (io.gt.spp%norbs) call die("symfio: No such orbital") 429 else if (io.lt.0) then 430 if (-io.gt.spp%nprojs) call die("symfio: No such projector") 431 else 432 symfio = 's' 433 endif 434 435 lorb=lofio(is,io) 436 morb=mofio(is,io) 437 438 if(lorb.gt.lmax_sym ) then 439 symfio=' ' 440 else 441 ilm=lorb*lorb+lorb+morb+1 442 if(pol(is,io)) then 443 symfio='P'//sym_label(ilm) 444 else 445 symfio=sym_label(ilm) 446 endif 447 endif 448 449 end function symfio 450! 451! End of FIOs ---------------------------------------------------- 452! 453 FUNCTION POL (IS,IO) 454 logical pol 455 integer, intent(in) :: is ! Species index 456 integer, intent(in) :: io ! Orbital index (within atom) 457 ! io>0 => basis orbitals 458 459C If true, the orbital IO is a perturbative polarization orbital 460 type(species_info), pointer :: spp 461 462 spp => species(is) 463 if ( (io .gt. species(is)%norbs) .or. 464 $ (io .le. 0)) call die("pol: Wrong io") 465 466 pol = spp%orbnl_ispol(spp%orb_index(io)) 467 468 end function pol 469 470 FUNCTION EPSKB (IS,IO) 471 real(dp) epskb 472 integer, intent(in) :: is ! Species index 473 integer, intent(in) :: io ! KB proyector index (within atom) 474 ! May be positive or negative 475 ! (only ABS(IO) is used). 476 477C Returns the energies epsKB_l of the Kleynman-Bylander projectors: 478C <Psi|V_KB|Psi'> = <Psi|V_local|Psi'> + 479C Sum_lm( epsKB_l * <Psi|Phi_lm> * <Phi_lm|Psi'> ) 480C where Phi_lm is returned by subroutine PHIATM. 481C Energy in Rydbergs. 482 type(species_info), pointer :: spp 483 484 integer ik 485 486 spp => species(is) 487 ik = abs(io) 488 if ((ik.gt.spp%nprojs) .or. 489 $ (ik .lt. 1) ) call die("epskb: No such projector") 490 epskb = spp%pjnl_ekb(spp%pj_index(ik)) 491 492 end function epskb 493 494!-------------------------------------------------------------------- 495 subroutine vna_sub(is,r,v,grv) 496 integer, intent(in) :: is ! Species index 497 real(dp), intent(in) :: r(3) ! Point vector, relative to atom 498 real(dp), intent(out) :: v ! Value of local pseudopotential 499 real(dp), intent(out) :: grv(3) ! Gradient of local pseudopotential 500 501C Returns local part of neutral-atom Kleynman-Bylander pseudopotential. 502C Distances in Bohr, Energies in Rydbergs 503C 2) Returns exactly zero when |R| > RCUT(IS,0) 504 type(rad_func), pointer :: func 505 506 real(dp) rmod, dvdr 507 508 call chk('vna_sub',is) 509 510 v = 0.0_dp 511 grv(1:3) = 0.0_dp 512 513 if (floating(is)) return 514 515 func => species(is)%vna 516 517 rmod = sqrt(sum(r*r)) 518 if (rmod .gt. func%cutoff) return 519 520 call rad_get(func,rmod,v,dvdr) 521 rmod = rmod + tiny20 522 grv(1:3) = dvdr * r(1:3)/rmod 523 524 end subroutine vna_sub 525 526 subroutine psch(is,r,ch,grch) 527 integer, intent(in) :: is ! Species index 528 real(dp), intent(in) :: r(3) ! Point vector, relative to atom 529 real(dp), intent(out) :: ch ! Local pseudopot. charge dens. 530 real(dp), intent(out) :: grch(3) ! Gradient of local ps. ch. dens. 531 532C Returns 'local-pseudotential charge density'. 533C Distances in Bohr, Energies in Rydbergs 534C Density in electrons/Bohr**3 535C 2) Returns exactly zero when |R| > Rchloc 536 type(rad_func), pointer :: func 537 538 real(dp) :: rmod, dchdr 539 540 call chk('psch',is) 541 542 ch = 0.0_dp 543 grch(1:3) = 0.0_dp 544 545 if (floating(is)) return 546 547 func => species(is)%chlocal 548 rmod = sqrt(sum(r*r)) 549 if (rmod .gt. func%cutoff) return 550 551 call rad_get(func,rmod,ch,dchdr) 552 rmod = rmod + tiny20 553 grch(1:3) = dchdr * r(1:3)/rmod 554 555 end subroutine psch 556 557 subroutine chcore_sub(is,r,ch,grch) 558 integer, intent(in) :: is ! Species index 559 real(dp), intent(in) :: r(3) ! Point vector, relative to atom 560 real(dp), intent(out) :: ch ! Value of pseudo-core charge dens. 561 real(dp), intent(out) :: grch(3) ! Gradient of pseudo-core ch. dens. 562 563C Returns returns pseudo-core charge density for non-linear core correction 564C in the xc potential. 565C Distances in Bohr, Energies in Rydbergs, Density in electrons/Bohr**3 566C 2) Returns exactly zero when |R| > Rcore 567 type(rad_func), pointer :: func 568 569 real(dp) rmod, dchdr 570 571 call chk('chcore_sub',is) 572 573 ch = 0.0_dp 574 grch(1:3) = 0.0_dp 575 576 if (floating(is)) return 577 578 func => species(is)%core 579 rmod = sqrt(sum(r*r)) 580 rmod = rmod + tiny20 ! Moved here. JMS, Dec.2012 581 if (rmod .gt. func%cutoff) return 582 583 call rad_get(func,rmod,ch,dchdr) 584! rmod=rmod+tiny20 ! Removed. JMS, Dec.2012 585 grch(1:3) = dchdr * r(1:3)/rmod 586 587 end subroutine chcore_sub 588 589 subroutine phiatm(is,io,r,phi,grphi) 590 integer, intent(in) :: is ! Species index 591 integer, intent(in) :: io ! Orbital index (within atom) 592! IO > 0 => Basis orbitals 593! IO = 0 => Local screened pseudopotential 594! IO < 0 => Kleynman-Bylander projectors 595 real(dp), intent(in) :: r(3) ! Point vector, relative to atom 596 real(dp), intent(out) :: phi ! Basis orbital, KB projector, or 597 ! local pseudopotential 598 real(dp), intent(out) :: grphi(3)! Gradient of BO, KB proj, or Loc ps 599 600C Returns Kleynman-Bylander local pseudopotential, nonlocal projectors, 601C and atomic basis orbitals (and their gradients). 602C Distances in Bohr 603C 1) Each projector and basis function has a well defined total 604C angular momentum (quantum number l). 605C 2) Basis functions are normalized and mutually orthogonal 606C 3) Projection functions are normalized and mutually orthogonal 607C 4) Normalization of KB projectors |Phi_lm> is such that 608C <Psi|V_KB|Psi'> = <Psi|V_local|Psi'> + 609C Sum_lm( epsKB_l * <Psi|Phi_lm> * <Phi_lm|Psi'> ) 610C where epsKB_l is returned by function EPSKB 611C 5) Prints a message and stops when no data exits for IS and/or IO 612C 6) Returns exactly zero when |R| > RCUT(IS,IO) 613C 7) PHIATM with IO = 0 is strictly equivalent to VNA_SUB 614 type(species_info), pointer :: spp 615 type(rad_func), pointer :: func 616 617 real(dp) rmod, phir, dphidr 618 real(dp) rly(max_ilm), grly(3,max_ilm) 619 integer l, m, ik, ilm 620 621 phi = 0.0_dp 622 grphi(1:3) = 0.0_dp 623 624 spp => species(is) 625 if (io.gt.0) then 626 if (io.gt.spp%norbs) call die("phiatm: No such orbital") 627 func => spp%orbnl(spp%orb_index(io)) 628 l = spp%orb_l(io) 629 m = spp%orb_m(io) 630 else if (io.lt.0) then 631 if (floating(is)) return 632 ik = -io 633 if (ik.gt.spp%nprojs) call die("phiatm: No such projector") 634 func => spp%pjnl(spp%pj_index(ik)) 635 l = spp%pj_l(ik) 636 m = spp%pj_m(ik) 637 else ! io=0 638 if (floating(is)) return 639 func => spp%vna 640 l = 0 641 m = 0 642 endif 643 644 rmod = sqrt(sum(r*r)) + tiny20 645 if(rmod.gt.func%cutoff-tiny12) return 646 647 call rad_get(func,rmod,phir,dphidr) 648 649 if (io.eq.0) then 650 phi=phir 651 grphi(1:3)=dphidr*r(1:3)/rmod 652 else 653 654 ilm = l*l + l + m + 1 655 call rlylm( l, r, rly, grly ) 656 phi = phir * rly(ilm) 657 grphi(1)=dphidr*rly(ilm)*r(1)/rmod+phir*grly(1,ilm) 658 grphi(2)=dphidr*rly(ilm)*r(2)/rmod+phir*grly(2,ilm) 659 grphi(3)=dphidr*rly(ilm)*r(3)/rmod+phir*grly(3,ilm) 660 661 endif 662 663 end subroutine phiatm 664 665 666 subroutine rphiatm(is,io,r,phi,dphidr) 667 integer, intent(in) :: is ! Species index 668 integer, intent(in) :: io ! Orbital index (within atom) 669! IO > 0 => Basis orbitals 670! IO = 0 => Local screened pseudopotential 671! IO < 0 => Kleynman-Bylander projectors 672 real(dp), intent(in) :: r ! Radial distance, relative to atom 673 real(dp), intent(out) :: phi ! Basis orbital, KB projector, or 674 ! local pseudopotential 675 real(dp), intent(out) :: dphidr ! Radial derivative of BO, 676 ! KB proj, or Loc pseudopot. 677 678C Returns the radial component of 679C Kleynman-Bylander local pseudopotential, nonlocal projectors, 680C and atomic basis orbitals (and their radial drivatives) 681C Distances in Bohr 682C 1) Each projector and basis function has a well defined total 683C angular momentum (quantum number l). 684C 2) Basis functions are normalized and mutually orthogonal 685C 3) Projection functions are normalized and mutually orthogonal 686C 4) Normalization of KB projectors |Phi_lm> is such that 687C <Psi|V_KB|Psi'> = <Psi|V_local|Psi'> + 688C Sum_lm( epsKB_l * <Psi|Phi_lm> * <Phi_lm|Psi'> ) 689C where epsKB_l is returned by function EPSKB 690C 6) Returns exactly zero when |R| > RCUT(IS,IO) 691C 7) RPHIATM with ITYPE = 0 is strictly equivalent to VNA_SUB 692 type(species_info), pointer :: spp 693 type(rad_func), pointer :: func 694 695 real(dp) rmod, phir 696 integer l, m, ik 697 698 phi = 0.0_dp 699 dphidr = 0._dp 700 701 spp => species(is) 702 if (io.gt.0) then 703 if (io.gt.spp%norbs) call die("rphiatm: No such orbital") 704 func => spp%orbnl(spp%orb_index(io)) 705 l = spp%orb_l(io) 706 m = spp%orb_m(io) 707 else if (io.lt.0) then 708 if (floating(is)) return 709 ik = -io 710 if (ik.gt.spp%nprojs) call die("rphiatm: No such projector") 711 func => spp%pjnl(spp%pj_index(ik)) 712 l = spp%pj_l(ik) 713 m = spp%pj_m(ik) 714 else 715 if (floating(is)) return 716 func => spp%vna 717 l = 0 718 m = 0 719 endif 720 721 rmod = r + tiny20 722 if(rmod.gt.func%cutoff-tiny12) return 723 724 call rad_get(func,rmod,phir,dphidr) 725 726 if (l.eq.0) then 727 phi=phir 728 elseif (l.eq.1) then 729 phi=phir*r 730 dphidr=dphidr*r 731 dphidr=dphidr+phir 732 else 733 phi=phir*r**l 734 dphidr=dphidr * r**l 735 dphidr=dphidr + l * phir * r**(l-1) 736 endif 737 738 end subroutine rphiatm 739 740 741 subroutine all_phi( is, it, r, maxnphi, nphi, phi, grphi ) 742 integer, intent(in) :: is ! Species index 743 integer, intent(in) :: it ! Orbital-type switch: 744 ! IT > 0 => Basis orbitals 745 ! IT < 0 => KB projectors 746 real(dp), intent(in) :: r(3) ! Point vector, relative to atom 747 integer, intent(in) :: maxnphi ! Maximum number of phi's 748 integer, intent(out) :: nphi ! Number of phi's 749 real(dp), intent(out) :: phi(maxnphi) ! Basis orbital, KB projector, or 750 ! local pseudopotential 751 real(dp), optional, intent(out) :: grphi(3,maxnphi) ! Gradient of phi 752 753C Returns Kleynman-Bylander local pseudopotential, nonlocal projectors, 754C and atomic basis orbitals (and their gradients). 755C Same as phiatm but returns all orbitals or KB projectors of the atom 756C Written by D.Sanchez-Portal and J.M.Soler. Jan. 2000 757 758C Distances in Bohr 759C 1) Each projector and basis function has a well defined total 760C angular momentum (quantum number l). 761C 2) Basis functions are normalized and mutually orthogonal 762C 3) Projection functions are normalized and mutually orthogonal 763C 4) Normalization of KB projectors |Phi_lm> is such that 764C <Psi|V_KB|Psi'> = <Psi|V_local|Psi'> + 765C Sum_lm( epsKB_l * <Psi|Phi_lm> * <Phi_lm|Psi'> ) 766C where epsKB_l is returned by function EPSKB 767C 5) Prints a message and stops when no data exits for IS 768C 6) Returns exactly zero when |R| > RCUT(IS,IO) 769C 8) If arrays phi or grphi are too small, returns with the required 770C value of nphi 771 type(species_info), pointer :: spp 772 773 integer i, jlm, l, lmax, m, maxlm 774 double precision rmod, phir, dphidr 775 real(dp) rly(max_ilm), grly(3,max_ilm) 776 777 integer :: ilm(maxnphi) 778 double precision :: rmax(maxnphi) 779 logical :: within(maxnphi) 780 781 call chk('all_phi',is) 782 spp => species(is) 783 784! Find number of orbitals 785 if (it.gt.0) then 786 nphi=spp%norbs 787 elseif (it.lt.0) then 788 nphi=spp%nprojs 789 else 790 call die("all_phi: Please use phiatm to get Vna...") 791 endif 792 793 if (nphi.gt.maxnphi) call die('all_phi: maxphi too small') 794 795 if (it.gt.0) then 796 do i = 1, nphi 797 l = spp%orb_l(i) 798 m = spp%orb_m(i) 799 ilm(i) = l*(l+1)+m+1 800 rmax(i) = spp%orbnl(spp%orb_index(i))%cutoff 801 enddo 802 else 803 do i = 1, nphi 804 rmax(i) = spp%pjnl(spp%pj_index(i))%cutoff 805 l = spp%pj_l(i) 806 m = spp%pj_m(i) 807 ilm(i) = l*(l+1)+m+1 808 enddo 809 endif 810 811! Initialize orbital values 812 phi(1:nphi) = 0._dp 813 if (present(grphi)) grphi(:,1:nphi) = 0._dp 814 815 if ((it.lt.0) .and. floating(is)) return 816 817! Find for which orbitals rmod < rmax and test for quick return 818 rmod = sqrt(sum(r*r)) + tiny20 819 within(1:nphi) = ( rmax(1:nphi) > rmod ) 820 if (.not.any(within(1:nphi))) return 821 822! Find spherical harmonics 823 maxlm = maxval( ilm(1:nphi), mask=within(1:nphi) ) 824 lmax=nint(sqrt(real(maxlm,dp)))-1 825 call rlylm(lmax,r,rly,grly) 826 827! Find values 828 829 i_loop: do i=1,nphi 830 831! Check if rmod > rmax 832 if (.not.within(i)) cycle i_loop 833 834! Find radial part 835 if (it.gt.0) then 836 call rad_get(spp%orbnl(spp%orb_index(i)),rmod,phir,dphidr) 837 else 838 call rad_get(spp%pjnl(spp%pj_index(i)),rmod,phir,dphidr) 839 end if 840 841! Multiply radial and angular parts 842 jlm = ilm(i) 843 phi(i) = phir * rly(jlm) 844 if (present(grphi)) 845 . grphi(:,i) = dphidr * rly(jlm) * r(:) / rmod + 846 . phir * grly(:,jlm) 847 848 enddo i_loop 849 850 end subroutine all_phi 851! 852! 853! This routine takes two species arguments 854! 855 subroutine psover(is1,is2,r,energ,dedr) 856 integer, intent(in) :: is1, is2 ! Species indexes 857 real(dp), intent(in) :: r ! Distance between atoms 858 real(dp), intent(out) :: energ ! Value of the correction 859 ! interaction energy 860 real(dp), intent(out) :: dedr ! Radial derivative of the correction 861 862C Returns electrostatic correction to the ions interaction energy 863C due to the overlap of the two 'local pseudopotential charge densities' 864C Distances in Bohr, Energies in Rydbergs 865C 2) Returns exactly zero when |R| > Rchloc 866 type(rad_func), pointer :: func 867 868 integer ismx, ismn, indx 869 real(dp) r_local 870 871 call chk('psover',is1) 872 call chk('psover',is2) 873 874 energ=0.0_dp 875 dedr=0.0_dp 876 877 if (floating(is1) .or. floating(is2)) return 878 879 ismx=max(is1,is2) 880 ismn=min(is1,is2) 881 indx=((ismx-1)*ismx)/2 + ismn 882 func => elec_corr(indx) 883 884 if ( r .gt. func%cutoff - tiny12 ) return 885 886 call rad_get(func,r,energ,dedr) 887 r_local = r+tiny20 888 energ=2.0_dp*energ/r_local 889 dedr=(-energ + 2.0_dp*dedr)/r_local 890 891 end subroutine psover 892 893! 894! Deprecated 895! 896 FUNCTION NZTFL (IS,L) 897 integer nztfl 898 integer, intent(in) :: is ! Species index 899 integer, intent(in) :: l ! Angular momentum of the basis funcs 900C Returns the number of different basis functions 901C with the same angular momentum and for a given species 902 type(species_info), pointer :: spp 903 904 integer i 905 906 call chk('nztfl',is) 907 spp => species(is) 908 909 nztfl = 0 910 do i = 1, spp%norbs 911 if (spp%orb_l(i).eq.l) nztfl = nztfl+1 912 enddo 913 914 end function nztfl 915 916 FUNCTION NKBL_FUNC (IS,L) 917 integer nkbl_func 918 integer, intent(in) :: is ! Species index 919 integer, intent(in) :: l ! Angular momentum of the basis funcs 920 921C Returns the number of different KB projectors 922C with the same angular momentum and for a given species 923 type(species_info), pointer :: spp 924 925 integer i 926 927 call chk('nkbl_func',is) 928 929 spp => species(is) 930 931 nkbl_func = 0 932 do i = 1, spp%nprojs 933 if (spp%pj_l(i).eq.l) nkbl_func = nkbl_func+1 934 enddo 935 936 end function nkbl_func 937 938 end module atmfuncs 939 940 941 942 943 944