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 old_atmfuncs 9 10C This module 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 14! AG: It is now a "legacy" module to interface the old and new data 15! structures. 16 17 use precision, only: dp 18 use sys 19 use atmparams, only: nzetmx, lmaxd, nsemx 20 use atmparams, only: maxos, nkbmx, ntbmax 21 use alloc, only: re_alloc, de_alloc 22 implicit none 23 24 integer, save, public :: nsmax 25 26 integer, public, pointer :: izsave(:) 27 integer, public, pointer :: nomax(:) 28 integer, public, pointer :: nkbmax(:) 29 30 integer, public, pointer :: lmxosave(:) 31 integer, public, pointer :: nkblsave(:,:) 32 integer, public, pointer :: npolorbsave(:,:,:) 33 integer, public, pointer :: nsemicsave(:,:) 34 integer, public, pointer :: nzetasave(:,:,:) 35 integer, public, pointer :: cnfigtb(:,:,:) 36 37 logical, public, pointer :: semicsave(:) 38 39 real(dp), public, pointer :: zvaltb(:) 40 integer, public, pointer :: lmxkbsave(:) 41 real(dp), public, pointer :: smasstb(:) 42 real(dp), public, pointer :: chargesave(:) 43 real(dp), public, pointer :: slfe(:) 44 real(dp), public, pointer :: lambdatb(:,:,:,:) 45 real(dp), public, pointer :: filtercuttb(:,:,:) 46 47 real(dp), public, pointer :: qtb(:,:) 48 49 real(dp), public, pointer :: table(:,:,:) 50 real(dp), public, pointer :: tabpol(:,:,:) 51 real(dp), public, pointer :: tab2(:,:,:) 52 real(dp), public, pointer :: tab2pol(:,:,:) 53 54 55 real(dp), public, pointer :: coretab(:,:,:) 56 real(dp), public, pointer :: chloctab(:,:,:) 57 real(dp), public, pointer :: vlocaltab(:,:,:) 58 real(dp), public, pointer :: rctb(:,:,:) 59 real(dp), public, pointer :: rcotb(:,:,:,:) 60 real(dp), public, pointer :: rcpoltb(:,:,:,:) 61 62 character(len=20), save, public, pointer :: label_save(:) 63 character(len=10), save, public, pointer :: basistype_save(:) 64 65! Public routines 66 public :: labelfis, izofis, zvalfis, nkbfis 67 public :: massfis, lomaxfis, nofis, lmxkbfis 68 public :: cnfigfio, lofio, mofio 69 public :: atmpopfio , epskb, rcut 70 public :: clear_tables, allocate_old_arrays 71 public :: deallocate_old_arrays 72 73 74 PRIVATE 75 76 CONTAINS !================================================ 77! 78 subroutine allocate_old_arrays() 79 80 !allocate(rcotb(nzetmx,0:lmaxd,nsemx,nsmax)) 81 nullify( rcotb ) 82 call re_alloc( rcotb, 1, nzetmx, 0, lmaxd, 1, nsemx, 1, nsmax, 83 & 'rcotb', 'old_atmfuncs' ) 84 85 !allocate(rcpoltb(nzetmx,0:lmaxd,nsemx,nsmax)) 86 nullify( rcpoltb ) 87 call re_alloc( rcpoltb, 1, nzetmx, 0, lmaxd, 1, nsemx, 1, nsmax, 88 & 'rcpoltb', 'old_atmfuncs' ) 89 !allocate(lambdatb(nzetmx,0:lmaxd,nsemx,nsmax)) 90 nullify( lambdatb ) 91 call re_alloc( lambdatb, 1, nzetmx, 0, lmaxd, 1, nsemx, 92 & 1, nsmax, 'lambdatb', 'old_atmfuncs' ) 93 !allocate(filtercuttb(0:lmaxd,nsemx,nsmax)) 94 nullify( filtercuttb ) 95 call re_alloc( filtercuttb, 0, lmaxd, 1, nsemx, 96 & 1, nsmax, 'filtercuttb', 'old_atmfuncs' ) 97 !allocate(qtb(maxos,nsmax)) 98 nullify( qtb ) 99 call re_alloc( qtb, 1, maxos, 1, nsmax, 100 & 'qtb', 'old_atmfuncs' ) 101 !allocate(slfe(nsmax)) 102 nullify( slfe ) 103 call re_alloc( slfe, 1, nsmax, 'slfe', 'old_atmfuncs' ) 104 !allocate(rctb(nkbmx,0:lmaxd,nsmax)) 105 nullify( rctb ) 106 call re_alloc( rctb, 1, nkbmx, 0, lmaxd, 1, nsmax, 107 & 'rctb', 'old_atmfuncs' ) 108 109 !allocate(smasstb(nsmax)) 110 nullify( smasstb ) 111 call re_alloc( smasstb, 1, nsmax, 'smasstb', 'old_atmfuncs' ) 112 !allocate(chargesave(nsmax)) 113 nullify( chargesave ) 114 call re_alloc( chargesave, 1, nsmax, 115 & 'chargesave', 'old_atmfuncs' ) 116! 117! Table: This is a hybrid 118! negative values of the second index correspond to KB projectors 119! positive values of the second index correspond to orbitals 120! A second index of zero (0) corresponds to the local potential 121! 122! The first index has ntbmax "real points" and two extra 123! entries for bookeeping 124! The total number of angular momentum entries is lmaxd+1 (since 125! s=0, p=1, etc) 126! 127! 128! allocate(table((ntbmax+2),-nkbmx*(lmaxd+1):nzetmx*nsemx* 129! (lmaxd+1),nsmax)) 130 nullify( table ) 131 call re_alloc( table, 1, ntbmax+2, 132 & -nkbmx*(lmaxd+1), nzetmx*nsemx*(lmaxd+1), 133 & 1, nsmax, 'table', 'old_atmfuncs' ) 134! allocate(tab2(ntbmax,-nkbmx*(lmaxd+1):nzetmx*nsemx* 135! (lmaxd+1),nsmax)) 136 nullify( tab2 ) 137 call re_alloc( tab2, 1, ntbmax, 138 & -nkbmx*(lmaxd+1), nzetmx*nsemx*(lmaxd+1), 139 & 1, nsmax, 'tab2', 'old_atmfuncs' ) 140! allocate(tabpol((ntbmax+2),nzetmx*nsemx*(lmaxd+1),nsmax)) 141 nullify( tabpol ) 142 call re_alloc( tabpol, 1, ntbmax+2, 143 & 1, nzetmx*nsemx*(lmaxd+1), 144 & 1, nsmax, 'tabpol', 'old_atmfuncs' ) 145! allocate(tab2pol(ntbmax,nzetmx*nsemx*(lmaxd+1),nsmax)) 146 nullify( tab2pol ) 147 call re_alloc( tab2pol, 1, ntbmax, 148 & 1, nzetmx*nsemx*(lmaxd+1), 149 & 1, nsmax, 'tab2pol', 'old_atmfuncs' ) 150 151! allocate(coretab(ntbmax+1,2,nsmax)) 152 nullify( coretab ) 153 call re_alloc( coretab, 1, ntbmax+1, 1, 2, 1, nsmax, 154 & 'coretab', 'old_atmfuncs' ) 155 156! allocate(chloctab((ntbmax+1),2,nsmax)) 157 nullify( chloctab ) 158 call re_alloc( chloctab, 1, ntbmax+1, 1, 2, 1, nsmax, 159 & 'chloctab', 'old_atmfuncs' ) 160! allocate(vlocaltab((ntbmax+1),2,nsmax)) 161 nullify( vlocaltab ) 162 call re_alloc( vlocaltab, 1, ntbmax+1, 1, 2, 1, nsmax, 163 & 'vlocaltab', 'old_atmfuncs' ) 164 165 !allocate(izsave(nsmax)) 166 nullify( izsave ) 167 call re_alloc( izsave, 1, nsmax, 'izsave', 'old_atmfuncs' ) 168 !allocate(lmxkbsave(nsmax)) 169 nullify( lmxkbsave ) 170 call re_alloc( lmxkbsave, 1, nsmax, 171 & 'lmxkbsave', 'old_atmfuncs' ) 172 !allocate(lmxosave(nsmax)) 173 nullify( lmxosave ) 174 call re_alloc( lmxosave, 1, nsmax, 'lmxosave', 'old_atmfuncs' ) 175 !allocate(npolorbsave(0:lmaxd,nsemx,nsmax)) 176 nullify( npolorbsave ) 177 call re_alloc( npolorbsave, 0, lmaxd, 1, nsemx, 1, nsmax, 178 & 'npolorbsave', 'old_atmfuncs' ) 179 !allocate(nsemicsave(0:lmaxd,nsmax)) 180 nullify( nsemicsave ) 181 call re_alloc( nsemicsave, 0, lmaxd, 1, nsmax, 182 & 'nsemicsave', 'old_atmfuncs' ) 183 !allocate(nzetasave(0:lmaxd,nsemx,nsmax)) 184 nullify( nzetasave ) 185 call re_alloc( nzetasave, 0, lmaxd, 1, nsemx, 1, nsmax, 186 & 'nzetasave', 'old_atmfuncs' ) 187 !allocate(nomax(nsmax)) 188 nullify( nomax ) 189 call re_alloc( nomax, 1, nsmax, 'nomax', 'old_atmfuncs' ) 190 !allocate(nkbmax(nsmax)) 191 nullify( nkbmax ) 192 call re_alloc( nkbmax, 1, nsmax, 'nkbmax', 'old_atmfuncs' ) 193 !allocate(zvaltb(nsmax)) 194 nullify( zvaltb ) 195 call re_alloc( zvaltb, 1, nsmax, 'zvaltb', 'old_atmfuncs' ) 196 !allocate(cnfigtb(0:lmaxd,nsemx,nsmax)) 197 nullify( cnfigtb ) 198 call re_alloc( cnfigtb, 0, lmaxd, 1, nsemx, 1, nsmax, 199 & 'cnfigtb', 'old_atmfuncs' ) 200 !allocate(nkblsave(0:lmaxd,nsmax)) 201 nullify( nkblsave ) 202 call re_alloc( nkblsave, 0, lmaxd, 1, nsmax, 203 & 'nkblsave', 'old_atmfuncs' ) 204! 205 nullify (label_save) 206 allocate(label_save(nsmax)) 207! call re_alloc(label_save,1,nsmax,"label_save", 208! $ routine= "allocate_old_arrays") 209 nullify (basistype_save) 210 allocate(basistype_save(nsmax)) 211! call re_alloc(basistype_save,1,nsmax,"basistype_save", 212! $ routine= "allocate_old_arrays") 213 nullify (semicsave) 214 call re_alloc( semicsave, 1, nsmax, 215 & 'semicsave', 'old_atmfuncs' ) 216 217 end subroutine allocate_old_arrays 218 219 subroutine deallocate_old_arrays() 220 221 call de_alloc( rcotb, 'rcotb', 'old_atmfuncs' ) 222 call de_alloc( rcpoltb, 'rcpoltb', 'old_atmfuncs' ) 223 call de_alloc( lambdatb, 'lambdatb', 'old_atmfuncs' ) 224 call de_alloc( filtercuttb, 'filtercuttb', 'old_atmfuncs' ) 225 call de_alloc( qtb, 'qtb', 'old_atmfuncs' ) 226 call de_alloc( slfe, 'slfe', 'old_atmfuncs' ) 227 call de_alloc( rctb, 'rctb', 'old_atmfuncs' ) 228 call de_alloc( smasstb, 'smasstb', 'old_atmfuncs' ) 229 call de_alloc( chargesave, 'chargesave', 'old_atmfuncs' ) 230 call de_alloc( table, 'table', 'old_atmfuncs' ) 231 call de_alloc( tab2, 'tab2', 'old_atmfuncs' ) 232 call de_alloc( tabpol, 'tabpol', 'old_atmfuncs' ) 233 call de_alloc( tab2pol, 'tab2pol', 'old_atmfuncs' ) 234 call de_alloc( coretab, 'coretab', 'old_atmfuncs' ) 235 call de_alloc( chloctab, 'chloctab', 'old_atmfuncs' ) 236 call de_alloc( vlocaltab, 'vlocaltab', 'old_atmfuncs' ) 237 call de_alloc( izsave, 'izsave', 'old_atmfuncs' ) 238 call de_alloc( lmxkbsave, 'lmxkbsave', 'old_atmfuncs' ) 239 call de_alloc( lmxosave, 'lmxosave', 'old_atmfuncs' ) 240 call de_alloc( npolorbsave, 'npolorbsave', 'old_atmfuncs' ) 241 call de_alloc( nsemicsave, 'nsemicsave', 'old_atmfuncs' ) 242 call de_alloc( nzetasave, 'nzetasave', 'old_atmfuncs' ) 243 call de_alloc( nomax, 'nomax', 'old_atmfuncs' ) 244 call de_alloc( nkbmax, 'nkbmax', 'old_atmfuncs' ) 245 call de_alloc( zvaltb, 'zvaltb', 'old_atmfuncs' ) 246 call de_alloc( cnfigtb, 'cnfigtb', 'old_atmfuncs' ) 247 call de_alloc( nkblsave, 'nkblsave', 'old_atmfuncs' ) 248 call de_alloc( semicsave, 'semicsave', 'old_atmfuncs' ) 249 deallocate( label_save ) 250! call de_alloc( label_save, 'label_save', 'old_atmfuncs' ) 251 deallocate( basistype_save ) 252! call de_alloc( basistype_save, 'basistype_save', 'old_atmfuncs' ) 253 254 end subroutine deallocate_old_arrays 255 256 subroutine clear_tables() 257 258 integer is 259 260 do is=1,nsmax 261 izsave(is)=0 262 lmxosave(is)=0 263 lmxkbsave(is)=0 264 label_save(is)=' ' 265 nkbmax(is)=0 266 nomax(is)=0 267 semicsave(is)=.false. 268 269 nsemicsave(:,is) = 0 270 nzetasave(:,:,is) = 0 271 rcotb(:,:,:,is) = 0.0_dp 272 lambdatb(:,:,:,is) = 0.0_dp 273 filtercuttb(:,:,is) = 0.0_dp 274 rcpoltb(:,:,:,is) = 0.0_dp 275 276 table(:,:,is) = 0.0_dp 277 tab2(:,:,is) = 0.0_dp 278 tabpol(:,:,is) = 0.0_dp 279 tab2pol(:,:,is) = 0.0_dp 280 281 qtb(1:maxos,is)=0.00_dp 282 283 enddo 284 end subroutine clear_tables 285 286 subroutine check_is(name,is) 287 character(len=*), intent(in) :: name 288 integer, intent(in) :: is 289 290 if ((is.lt.1).or.(is.gt.nsmax)) then 291 write(6,*) trim(name),': THERE ARE NO DATA FOR IS=',IS 292 write(6,*) trim(name),': ISMIN= 1, NSMAX= ',nsmax 293 call die() 294 endif 295 end subroutine check_is 296! 297! 298 FUNCTION IZOFIS( IS ) 299 integer :: izofis ! Atomic number 300 integer, intent(in) :: is ! Species index 301 302 call check_is('izofis',is) 303 izofis=izsave(is) 304 305 end function izofis 306! 307 FUNCTION ZVALFIS( IS ) 308 real(dp) :: zvalfis ! Valence charge 309 integer, intent(in) :: is ! Species index 310 311 call check_is('zvalfis',is) 312 313 zvalfis=zvaltb(is) 314 end function zvalfis 315! 316 FUNCTION LABELFIS (IS) 317 character(len=20) :: labelfis ! Atomic label 318 integer, intent(in) :: is ! Species index 319 320 call check_is('labelfis',is) 321 labelfis=label_save(is) 322 end function labelfis 323! 324 FUNCTION LMXKBFIS (IS) 325 integer :: lmxkbfis ! Maximum ang mom of the KB projectors 326 integer, intent(in) :: is ! Species index 327 328 call check_is('lmxkbfis',is) 329 lmxkbfis=lmxkbsave(is) 330 end function lmxkbfis 331! 332 FUNCTION LOMAXFIS (IS) 333 integer :: lomaxfis ! Maximum ang mom of the Basis Functions 334 integer, intent(in) :: is ! Species index 335 336 integer lmx, nsm 337 338 call check_is('lomaxfis',is) 339 340 lomaxfis=0 341 lmx=lmxosave(is) 342 do nsm=1,nsemicsave(lmx,is)+1 343 if(npolorbsave(lmx,nsm,is).gt.0) lomaxfis=lmx+1 344 enddo 345 346 lomaxfis=max(lomaxfis,lmx) 347 end function lomaxfis 348! 349 350 FUNCTION MASSFIS(IS) 351 real(dp) :: massfis ! Mass 352 integer, intent(in) :: is ! Species index 353 354 call check_is('massfis',is) 355 massfis=smasstb(is) 356 end function massfis 357! 358 FUNCTION NKBFIS(IS) 359 integer :: nkbfis ! Total number of KB projectors 360 integer, intent(in) :: is ! Species index 361 362 call check_is('nkbfis',is) 363 nkbfis=nkbmax(is) 364 end function nkbfis 365! 366 367 FUNCTION NOFIS(IS) 368 integer :: nofis ! Total number of Basis functions 369 integer, intent(in) :: is ! Species index 370 371 call check_is('nofis',is) 372 nofis=nomax(is) 373 end function nofis 374 375!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! AMENoFIS 376 377 FUNCTION ATMPOPFIO (IS,IO) 378 real(dp) atmpopfio 379 integer, intent(in) :: is ! Species index 380 integer, intent(in) :: io ! Orbital index (within atom) 381 382C Returns the population of the atomic basis orbitals in the atomic 383C ground state configuration. 384 385 call check_is('atmpopfio',is) 386 if((io.gt.nomax(is)).or.(io.lt.1)) then 387 write(6,*) 'ATMPOPFIO: THERE ARE NO DATA FOR IO=',IO 388 write(6,*) 'ATMPOPFIO: IOMIN= 1', ' IOMAX= ',nomax(is) 389 call die 390 endif 391 392 atmpopfio=qtb(io,is) 393 end function atmpopfio 394 395! 396! 397 FUNCTION CNFIGFIO(IS,IO) 398 integer cnfigfio 399 integer, intent(in) :: is ! Species index 400 integer, intent(in) :: io ! Orbital index (within atom) 401 402C INTEGER CNFIGFIO: Principal quantum number of the shell to which 403C the orbital belongs 404 405C (Formerly, for polarization orbitals 406C the quantum number corresponded to the shell which 407C is polarized by the orbital io. 408C This behavior for polarization orbitals is meaningless, 409C and has been replaced.) 410 411 integer l, norb, izeta, ipol, nsm 412 413 call check_is('cnfigfio',is) 414 if ((io.gt.nomax(is)).or.(io.lt.1)) then 415 write(6,*) 'CNFIGFIO: THERE ARE NO DATA FOR IO=',IO 416 write(6,*) 'CNFIGFIO: IOMIN= 1', 417 . ' IOMAX= ',nomax(is) 418 call die() 419 endif 420 421 ! This routine assumes that polarization orbitals are the last ones 422 ! in the list indexed by 'io' 423 424 norb=0 425 do 10 l=0,lmxosave(is) 426 do 8 nsm=1,nsemicsave(l,is)+1 427 do 5 izeta=1,nzetasave(l,nsm,is) 428 norb=norb+(2*l+1) 429 if(norb.ge.io) then 430 cnfigfio=cnfigtb(l,nsm,is) 431 return 432 endif 433 5 continue 434 8 continue 43510 continue 436 437 do 20 l=0,lmxosave(is) 438 do 18 nsm=1,nsemicsave(l,is)+1 439 do 15 ipol=1, npolorbsave(l,nsm,is) 440 norb=norb+(2*(l+1)+1) 441 if(norb.ge.io) then 442 ! Deal properly with polarization orbitals 443 ! Return the highest n for l+1 444 cnfigfio=maxval(cnfigtb(l+1,:,is)) 445 return 446 endif 44715 continue 44818 continue 44920 continue 450 write(6,*) 'CNFIGFIO: ERROR: ORBITAL INDEX IO=',IO 451 write(6,*) 'CNFIGFIO: ERROR: NOT FOUND' 452 call die() 453 454 end function cnfigfio 455! 456! 457 FUNCTION LOFIO (IS,IO) 458 integer lofio 459 integer, intent(in) :: is ! Species index 460 integer, intent(in) :: io ! Orbital index (within atom) 461 462C Returns total angular momentum quantum number of a given atomic basis 463C basis orbital or Kleynman-Bylander projector. 464 465C INTEGER IO : Orbital index (within atom) 466C IO > 0 => Basis orbitals 467C IO < 0 => Kleynman-Bylander projectors 468C************************OUTPUT***************************************** 469C INTEGER LOFIO : Quantum number L of orbital or KB projector 470 471 integer l, norb, izeta, ipol, nkb, nsm 472 473 call check_is('lofio',is) 474 if ((io.gt.nomax(is)).or.(io.lt.-nkbmax(is))) then 475 write(6,*) 'LOFIO: THERE ARE NO DATA FOR IO=',IO 476 write(6,*) 'LOFIO: IOMIN= ',-nkbmax(is), 477 . ' IOMAX= ',nomax(is) 478 CALL DIE 479 endif 480 481 if (io.gt.0) then 482 483 norb=0 484 do 10 l=0,lmxosave(is) 485 do 8 nsm=1,nsemicsave(l,is)+1 486 do 5 izeta=1,nzetasave(l,nsm,is) 487 norb=norb+(2*l+1) 488 if(norb.ge.io) goto 30 489 5 continue 490 8 continue 49110 continue 492 493 do 20 l=0,lmxosave(is) 494 do 18 nsm=1,nsemicsave(l,is)+1 495 do 15 ipol=1, npolorbsave(l,nsm,is) 496 norb=norb+(2*(l+1)+1) 497 if(norb.ge.io) goto 40 49815 continue 49918 continue 50020 continue 501 write(6,*) 'LOFIO: ERROR: ORBITAL INDEX IO=',IO 502 write(6,*) 'LOFIO: ERROR: NOT FOUND' 503 call die 504 50530 lofio=l 506 return 507 50840 lofio=l+1 509 return 510 511 elseif(io.lt.0) then 512 513 nkb=0 514 do 50 l=0,lmxkbsave(is) 515 do 45 izeta=1,nkblsave(l,is) 516 nkb=nkb-(2*l+1) 517 if(nkb.le.io) goto 60 51845 continue 51950 continue 520 52160 lofio=l 522 523c elseif (io.eq.0) then 524 else 525 526 lofio=0 527 528 endif 529 end function lofio 530! 531! 532 FUNCTION MOFIO (IS,IO) 533 integer mofio 534 integer, intent(in) :: is ! Species index 535 integer, intent(in) :: io ! Orbital index (within atom) 536 537C Returns magnetic quantum number of a given atomic basis 538C basis orbital or Kleynman-Bylander projector. 539 540C INTEGER IO : Orbital index (within atom) 541C IO > 0 => Basis orbitals 542C IO < 0 => Kleynman-Bylander projectors 543C************************OUTPUT***************************************** 544C INTEGER MOFIO : Quantum number M of orbital or KB projector 545 546 integer l, norb, izeta, ipol, nkb, lorb, lkb, nsm 547 548 call check_is('mofio',is) 549 if((io.gt.nomax(is)).or.(io.lt.-nkbmax(is))) then 550 write(6,*) 'MOFIO: THERE ARE NO DATA FOR IO=',IO 551 write(6,*) 'MOFIO: IOMIN= ',-nkbmax(is), 552 . ' IOMAX= ',nomax(is) 553 CALL DIE 554 endif 555 556 if (io.gt.0) then 557 558 norb=0 559 do 10 l=0,lmxosave(is) 560 do 8 nsm=1,nsemicsave(l,is)+1 561 do 5 izeta=1,nzetasave(l,nsm,is) 562 norb=norb+(2*l+1) 563 if(norb.ge.io) goto 30 564 5 continue 565 8 continue 56610 continue 567 568 do 20 l=0,lmxosave(is) 569 do 18 nsm=1,nsemicsave(l,is)+1 570 do 15 ipol=1, npolorbsave(l,nsm,is) 571 norb=norb+(2*(l+1)+1) 572 if(norb.ge.io) goto 40 57315 continue 57418 continue 57520 continue 576 write(6,*) 'MOFIO: ERROR: ORBITAL INDEX IO=',IO 577 write(6,*) 'MOFIO: ERROR: NOT FOUND' 578 call die() 579 58030 lorb=l 581 mofio=io-norb+lorb 582 return 583 58440 lorb=l+1 585 mofio=io-norb+lorb 586 return 587 588 589 elseif(io.lt.0) then 590 591 592 nkb=0 593 do 50 l=0,lmxkbsave(is) 594 do 45 izeta=1,nkblsave(l,is) 595 nkb=nkb-(2*l+1) 596 if(nkb.le.io) goto 60 59745 continue 59850 continue 599 60060 lkb=l 601 mofio=-io+nkb+lkb 602c elseif (io.eq.0) then 603 else 604 605 mofio=0 606 607 endif 608 609 end function mofio 610! 611 612! End of FIOs 613! 614 615 FUNCTION EPSKB (IS,IO) 616 real(dp) epskb 617 integer, intent(in) :: is ! Species index 618 integer, intent(in) :: io ! KB proyector index (within atom) 619 ! May be positive or negative 620 ! (only ABS(IO) is used). 621 622C Returns the energies epsKB_l of the Kleynman-Bylander projectors: 623C <Psi|V_KB|Psi'> = <Psi|V_local|Psi'> + 624C Sum_lm( epsKB_l * <Psi|Phi_lm> * <Phi_lm|Psi'> ) 625C where Phi_lm is returned by subroutine PHIATM. 626 627 628C REAL(DP) EPSKB : Kleynman-Bylander projector energy 629C Energy in Rydbergs. 630 631 integer ik, nkb, indx, l, ikb 632C 633C****************************************************************** 634C*****************Variables in the common blocks******************* 635C 636 call check_is('epskb',is) 637 638 ik=-abs(io) 639 if (ik.eq.0) then 640 write(6,*) 'EPSKB: FUNCTION CANNOT BE CALLED WITH' 641 . ,' ARGUMENT EQUAL TO ZERO' 642 CALL DIE 643 endif 644 645 if(ik.lt.-nkbmax(is)) then 646 write(6,*) 'EPSKB: THERE ARE NO DATA FOR IO=',IK 647 write(6,*) 'EPSKB: IOMIN= ',-nkbmax(is) 648 CALL DIE() 649 endif 650 651 nkb=0 652 indx=0 653 do 10 l=0,lmxkbsave(is) 654 do 5 ikb=1,nkblsave(l,is) 655 indx=indx+1 656 nkb=nkb-(2*l+1) 657 if(nkb.le.ik) goto 20 6585 continue 65910 continue 660 66120 indx=-indx 662 663 epskb=table(2,indx,is) 664 665 end function epskb 666! 667! 668 function rcut(is,io) 669 real(dp) rcut 670 integer, intent(in) :: is ! Species index 671 integer, intent(in) :: io ! Orbital index (within atom) 672 ! io> => basis orbitals 673 ! io<0 => KB projectors 674 ! io=0 : Local screened pseudopotential 675 676C Returns cutoff radius of Kleynman-Bylander projectors and 677C atomic basis orbitals. 678C Distances in Bohr 679 680 integer l, norb, lorb, nzetorb, izeta, ipol, nkb,lkb,nsm 681 integer indx, nsmorb 682C 683 call check_is('rcut',is) 684 685 if ((io.gt.nomax(is)).or.(io.lt.-nkbmax(is))) then 686 write(6,*) 'RCUT: THERE ARE NO DATA FOR IO=',IO 687 write(6,*) 'RCUT: IOMIN= ',-nkbmax(is), 688 . ' IOMAX= ',nomax(is) 689 CALL DIE 690 endif 691 692 693 if (io.gt.0) then 694 695 norb=0 696 indx=0 697 do 10 l=0,lmxosave(is) 698 do 8 nsm=1,nsemicsave(l,is)+1 699 do 5 izeta=1,nzetasave(l,nsm,is) 700 norb=norb+(2*l+1) 701 indx=indx+1 702 if(norb.ge.io) goto 30 703 5 continue 704 8 continue 70510 continue 706 707 indx=0 708 do 20 l=0,lmxosave(is) 709 do 18 nsm=1,nsemicsave(l,is)+1 710 do 15 ipol=1, npolorbsave(l,nsm,is) 711 norb=norb+(2*(l+1)+1) 712 indx=indx+1 713 if(norb.ge.io) goto 40 71415 continue 71518 continue 71620 continue 717 write(6,*) 'RCUT: ERROR: ORBITAL INDEX IO=',IO 718 write(6,*) 'RCUT: ERROR: NOT FOUND' 719 call die() 720 72130 lorb=l 722 nzetorb=izeta 723 nsmorb=nsm 724 rcut=rcotb(nzetorb,lorb,nsmorb,is) 725c rcut=table(1,indx,is)*(ntbmax-1) 726 return 727 72840 lorb=l 729 nzetorb=ipol 730 nsmorb=nsm 731 rcut=rcpoltb(nzetorb,lorb,nsmorb,is) 732c rcut=tabpol(1,indx,is)*(ntbmax-1) 733 return 734 735 736 elseif(io.lt.0) then 737 738 739 nkb=0 740 do 50 l=0,lmxkbsave(is) 741 do 45 izeta=1,nkblsave(l,is) 742 nkb=nkb-(2*l+1) 743 if(nkb.le.io) goto 60 74445 continue 74550 continue 746 74760 lkb=l 748 indx=-(lkb+1) 749 rcut=rctb(izeta,lkb,is) 750c rcut=table(1,indx,is)*(ntbmax-1) 751 return 752 753c elseif (io.eq.0) then 754 else 755 756 rcut=table(2,0,is) 757 758 endif 759 760 end function rcut 761! 762! 763 end module old_atmfuncs 764