1!! Copyright (C) 2002-2012 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module ps_oct_m 22 use atomic_oct_m 23 use global_oct_m 24 use io_oct_m 25 use lalg_adv_oct_m 26 use loct_math_oct_m 27 use parser_oct_m 28 use logrid_oct_m 29 use messages_oct_m 30 use namespace_oct_m 31 use profiling_oct_m 32 use ps_cpi_oct_m 33 use ps_fhi_oct_m 34 use ps_hgh_oct_m 35 use ps_xml_oct_m 36 use ps_in_grid_oct_m 37#ifdef HAVE_PSPIO 38 use fpspio_m 39#endif 40 use ps_psf_oct_m 41 use pseudo_oct_m 42 use splines_oct_m 43 use spline_filter_oct_m 44 implicit none 45 46 private 47 public :: & 48 ps_t, & 49 ps_init, & 50 ps_pspio_init, & 51 ps_separate, & 52 ps_filter, & 53 ps_getradius, & 54 ps_derivatives, & 55 ps_debug, & 56 ps_niwfs, & 57 ps_bound_niwfs, & 58 ps_end, & 59 ps_has_density, & 60 ps_has_nlcc, & 61 ps_density_volume 62 63 integer, parameter, public :: & 64 PS_FILTER_NONE = 0, & 65 PS_FILTER_TS = 2, & 66 PS_FILTER_BSB = 3 67 68 integer, public, parameter :: & 69 PROJ_NONE = 0, & 70 PROJ_HGH = 1, & 71 PROJ_KB = 2, & 72 PROJ_RKB = 3 73 74 integer, parameter, public :: INVALID_L = 333 75 76 character(len=4), parameter :: ps_name(PSEUDO_FORMAT_UPF1:PSEUDO_FORMAT_PSP8) = & 77 (/"upf1", "upf2", "qso ", "psml", "psf ", "cpi ", "fhi ", "hgh ", "psp8"/) 78 79 type ps_t 80 ! Components are public by default 81 integer :: projector_type 82 character(len=10), private :: label 83 84 integer, private :: ispin !< Consider spin (ispin = 2) or not (ispin = 1) 85 FLOAT, private :: z 86 FLOAT :: z_val 87 type(valconf_t) :: conf 88 type(logrid_t), private :: g 89 type(spline_t), pointer :: ur(:, :) !< (1:conf%p, 1:ispin) atomic wavefunctions, as a function of r 90 type(spline_t), pointer, private :: ur_sq(:, :) !< (1:conf%p, 1:ispin) atomic wavefunctions, as a function of r^2 91 logical, allocatable :: bound(:, :) !< (1:conf%p, 1:ispin) is the state bound or not 92 93 ! Kleinman-Bylander projectors stuff 94 integer :: lmax !< maximum value of l to take 95 integer :: llocal !< which component to take as local 96 97 type(spline_t) :: vl !< local part 98 99 FLOAT :: projectors_sphere_threshold !< The projectors are localized in real 100 !! space, and so they are contained in a 101 !! sphere whose radius is computed by 102 !! making sure that the projector 103 !! functions absolute value is below this 104 !! threshold, for points outside the 105 !! sphere. 106 FLOAT :: rc_max !< The radius of the spheres that contain the projector functions. 107 108 integer :: kbc !< Number of KB components (1 or 2 for TM ps, 3 for HGH) 109 FLOAT, pointer :: h(:,:,:), k(:, :, :) 110 type(spline_t), pointer :: kb(:, :) !< Kleinman-Bylander projectors 111 type(spline_t), pointer :: dkb(:, :) !< derivatives of KB projectors 112 113 logical :: nlcc !< .true. if the pseudo has non-linear core corrections. 114 type(spline_t) :: core !< normalization \int dr 4 pi r^2 rho(r) = N 115 type(spline_t) :: core_der !< derivative of the core correction 116 117 118 !LONG-RANGE PART OF THE LOCAL POTENTIAL 119 120 logical, private :: has_long_range 121 122 type(spline_t), private :: vlr !< the long-range part of the local potential 123 type(spline_t) :: vlr_sq !< the long-range part of the 124 !< local potential in terms of r^2, to avoid the sqrt 125 type(spline_t) :: nlr !< the charge density associated with the long-range part 126 127 FLOAT :: sigma_erf !< the a constant in erf(r/(sqrt(2)*sigma))/r 128 129 logical, private :: has_density !< does the species have a density? 130 type(spline_t), pointer :: density(:) !< the atomic density for each spin 131 type(spline_t), pointer :: density_der(:) !< the radial derivative for the atomic density for each spin 132 133 logical, private :: is_separated 134 logical :: local 135 logical :: hamann 136 integer, private :: file_format 137 integer, private :: pseudo_type 138 integer :: exchange_functional 139 integer :: correlation_functional 140 end type ps_t 141 142 FLOAT, parameter :: eps = CNST(1.0e-8) 143 144contains 145 146 147 ! --------------------------------------------------------- 148 subroutine ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename) 149 type(ps_t), intent(out) :: ps 150 type(namespace_t), intent(in) :: namespace 151 character(len=10), intent(in) :: label 152 integer, intent(in) :: user_lmax 153 integer, intent(in) :: user_llocal 154 integer, intent(in) :: ispin 155 FLOAT, intent(in) :: z 156 character(len=*), intent(in) :: filename 157 158 integer :: l, ii, ll, is, ierr 159 type(ps_psf_t) :: ps_psf !< SIESTA pseudopotential 160 type(ps_cpi_t) :: ps_cpi !< Fritz-Haber pseudopotential 161 type(ps_fhi_t) :: ps_fhi !< Fritz-Haber pseudopotential (from abinit) 162 type(ps_hgh_t) :: ps_hgh !< In case Hartwigsen-Goedecker-Hutter ps are used. 163 type(ps_xml_t) :: ps_xml !< For xml based pseudopotentials 164 logical, save :: xml_warned = .false. 165 FLOAT, allocatable :: eigen(:, :) !< eigenvalues 166 167 PUSH_SUB(ps_init) 168 169 ps%exchange_functional = PSEUDO_EXCHANGE_UNKNOWN 170 ps%correlation_functional = PSEUDO_CORRELATION_UNKNOWN 171 172 ! Fix the threshold to calculate the radius of the projector-function localization spheres: 173 174 call messages_obsolete_variable(namespace, 'SpecieProjectorSphereThreshold', 'SpeciesProjectorSphereThreshold') 175 176 !%Variable SpeciesProjectorSphereThreshold 177 !%Type float 178 !%Default 0.001 179 !%Section System::Species 180 !%Description 181 !% The pseudopotentials may be composed of a local part, and a linear combination of nonlocal 182 !% operators. These nonlocal projectors have "projector" form, <math> \left| v \right> \left< v \right| </math> 183 !% (or, more generally speaking, <math> \left| u \right> \left< v \right| </math>). 184 !% These projectors are localized in real space -- that is, the function <math>v</math> 185 !% has a finite support around the nucleus. This region where the projectors are localized should 186 !% be small or else the computation time required to operate with them will be very large. 187 !% 188 !% In practice, this localization is fixed by requiring the definition of the projectors to be 189 !% contained in a sphere of a certain radius. This radius is computed by making sure that the 190 !% absolute value of the projector functions, at points outside the localization sphere, is 191 !% below a certain threshold. This threshold is set by <tt>SpeciesProjectorSphereThreshold</tt>. 192 !%End 193 call parse_variable(namespace, 'SpeciesProjectorSphereThreshold', CNST(0.001), ps%projectors_sphere_threshold) 194 if(ps%projectors_sphere_threshold <= M_ZERO) call messages_input_error(namespace, 'SpeciesProjectorSphereThreshold') 195 196 ps%file_format = pseudo_detect_format(filename) 197 198 if(ps%file_format == PSEUDO_FORMAT_FILE_NOT_FOUND) then 199 call messages_write("Cannot open pseudopotential file '"//trim(filename)//"'.") 200 call messages_fatal(namespace=namespace) 201 end if 202 203 if(ps%file_format == PSEUDO_FORMAT_UNKNOWN) then 204 call messages_write("Cannot determine the pseudopotential type for species '"//trim(label)//"' from", & 205 new_line = .true.) 206 call messages_write("file '"//trim(filename)//"'.") 207 call messages_fatal(namespace=namespace) 208 end if 209 210 ps%label = label 211 ps%ispin = ispin 212 ps%hamann = .false. 213 ps%projector_type = PROJ_KB 214 215 select case(ps%file_format) 216 case(PSEUDO_FORMAT_PSF, PSEUDO_FORMAT_HGH) 217 ps%has_density = .true. 218 case default 219 ps%has_density = .false. 220 end select 221 222 select case(ps%file_format) 223 case(PSEUDO_FORMAT_PSF) 224 ps%pseudo_type = PSEUDO_TYPE_SEMILOCAL 225 226 call ps_psf_init(ps_psf, ispin, filename, namespace) 227 228 call valconf_copy(ps%conf, ps_psf%conf) 229 ps%z = z 230 ps%conf%z = nint(z) ! atomic number 231 ps%kbc = 1 ! only one projector per angular momentum 232 233 ps%lmax = ps_psf%ps_grid%no_l_channels - 1 234 235 if(user_lmax /= INVALID_L) then 236 ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components. 237 if(user_lmax /= ps%lmax) then 238 message(1) = "lmax in Species block for " // trim(label) // & 239 " is larger than number available in pseudopotential." 240 call messages_fatal(1, namespace=namespace) 241 end if 242 end if 243 244 ps%conf%p = ps_psf%ps_grid%no_l_channels 245 if(ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0. 246 247 ! the local part of the pseudo 248 if(user_llocal == INVALID_L) then 249 ps%llocal = 0 250 else 251 ps%llocal = user_llocal 252 end if 253 254 call ps_psf_process(ps_psf, namespace, ps%lmax, ps%llocal) 255 call logrid_copy(ps_psf%ps_grid%g, ps%g) 256 257 case(PSEUDO_FORMAT_CPI, PSEUDO_FORMAT_FHI) 258 ps%pseudo_type = PSEUDO_TYPE_SEMILOCAL 259 260 call valconf_null(ps%conf) 261 262 if(ps%file_format == PSEUDO_FORMAT_CPI) then 263 call ps_cpi_init(ps_cpi, trim(filename), namespace) 264 ps%conf%p = ps_cpi%ps_grid%no_l_channels 265 else 266 call ps_fhi_init(ps_fhi, trim(filename), namespace) 267 ps%conf%p = ps_fhi%ps_grid%no_l_channels 268 end if 269 270 ps%conf%z = nint(z) 271 ps%conf%symbol = label(1:2) 272 ps%conf%type = 1 273 do l = 1, ps%conf%p 274 ps%conf%l(l) = l - 1 275 end do 276 277 ps%z = z 278 ps%kbc = 1 ! only one projector per angular momentum 279 280 ps%lmax = ps%conf%p - 1 281 282 if(user_lmax /= INVALID_L) then 283 ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components. 284 if(user_lmax /= ps%lmax) then 285 message(1) = "lmax in Species block for " // trim(label) // & 286 " is larger than number available in pseudopotential." 287 call messages_fatal(1, namespace=namespace) 288 end if 289 end if 290 291 if(ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0. 292 293 ! the local part of the pseudo 294 if(user_llocal == INVALID_L) then 295 ps%llocal = 0 296 else 297 ps%llocal = user_llocal 298 end if 299 300 if(ps%file_format == PSEUDO_FORMAT_CPI) then 301 call ps_cpi_process(ps_cpi, ps%llocal, namespace) 302 call logrid_copy(ps_cpi%ps_grid%g, ps%g) 303 else 304 call ps_fhi_process(ps_fhi, ps%lmax, ps%llocal, namespace) 305 call logrid_copy(ps_fhi%ps_grid%g, ps%g) 306 end if 307 308 case(PSEUDO_FORMAT_HGH) 309 ps%pseudo_type = PSEUDO_TYPE_KLEINMAN_BYLANDER 310 ps%projector_type = PROJ_HGH 311 312 call hgh_init(ps_hgh, trim(filename), namespace) 313 call valconf_copy(ps%conf, ps_hgh%conf) 314 315 ps%z = z 316 ps%kbc = 3 317 ps%llocal = -1 318 ps%lmax = ps_hgh%l_max 319 320 call hgh_process(ps_hgh, namespace) 321 call logrid_copy(ps_hgh%g, ps%g) 322 323 case(PSEUDO_FORMAT_QSO, PSEUDO_FORMAT_UPF1, PSEUDO_FORMAT_UPF2, PSEUDO_FORMAT_PSML, PSEUDO_FORMAT_PSP8) 324 325 if(.not. xml_warned) then 326 call messages_experimental('XML (QSO, UPF, and PSML, PSP8) pseudopotential support') 327 xml_warned = .true. 328 end if 329 330 call ps_xml_init(ps_xml, namespace, trim(filename), ps%file_format, ierr) 331 332 ps%pseudo_type = pseudo_type(ps_xml%pseudo) 333 ps%exchange_functional = pseudo_exchange(ps_xml%pseudo) 334 ps%correlation_functional = pseudo_correlation(ps_xml%pseudo) 335 336 call valconf_null(ps%conf) 337 338 ps%z = z 339 ps%conf%z = nint(z) 340 341 if(ps_xml%kleinman_bylander) then 342 ps%conf%p = ps_xml%nwavefunctions 343 else 344 ps%conf%p = ps_xml%lmax + 1 345 end if 346 347 do ll = 0, ps_xml%lmax 348 ps%conf%l(ll + 1) = ll 349 end do 350 351 ps%kbc = ps_xml%nchannels 352 ps%lmax = ps_xml%lmax 353 354 if(ps_xml%kleinman_bylander) then 355 ps%llocal = ps_xml%llocal 356 else 357 ! we have several options 358 ps%llocal = 0 ! the default 359 if(ps_xml%llocal >= 0) ps%llocal = ps_xml%llocal ! the one given in the pseudopotential file 360 if(user_llocal /= INVALID_L) ps%llocal = user_llocal ! user supplied local component 361 ASSERT(ps%llocal >= 0) 362 ASSERT(ps%llocal <= ps%lmax) 363 end if 364 365 nullify(ps%g%drdi, ps%g%s) 366 367 ps%g%nrval = ps_xml%grid_size 368 369 SAFE_ALLOCATE(ps%g%rofi(1:ps%g%nrval)) 370 SAFE_ALLOCATE(ps%g%r2ofi(1:ps%g%nrval)) 371 372 do ii = 1, ps%g%nrval 373 ps%g%rofi(ii) = ps_xml%grid(ii) 374 ps%g%r2ofi(ii) = ps_xml%grid(ii)**2 375 end do 376 377 end select 378 379 ps%local = (ps%lmax == 0 .and. ps%llocal == 0 ) .or. (ps%lmax == -1 .and. ps%llocal == -1) 380 381 ! We allocate all the stuff 382 SAFE_ALLOCATE(ps%kb (0:ps%lmax, 1:ps%kbc)) 383 SAFE_ALLOCATE(ps%dkb (0:ps%lmax, 1:ps%kbc)) 384 SAFE_ALLOCATE(ps%ur (1:ps%conf%p, 1:ps%ispin)) 385 SAFE_ALLOCATE(ps%ur_sq(1:ps%conf%p, 1:ps%ispin)) 386 SAFE_ALLOCATE(ps%bound(1:ps%conf%p, 1:ps%ispin)) 387 SAFE_ALLOCATE(ps%h (0:ps%lmax, 1:ps%kbc, 1:ps%kbc)) 388 SAFE_ALLOCATE(ps%density(1:ps%ispin)) 389 SAFE_ALLOCATE(ps%density_der(1:ps%ispin)) 390 391 nullify(ps%k) 392 393 call spline_init(ps%kb) 394 call spline_init(ps%dkb) 395 call spline_init(ps%vl) 396 call spline_init(ps%core) 397 call spline_init(ps%core_der) 398 call spline_init(ps%density) 399 call spline_init(ps%density_der) 400 401 SAFE_ALLOCATE(eigen(1:ps%conf%p, 1:ps%ispin)) 402 eigen = M_ZERO 403 404 ! Now we load the necessary information. 405 select case(ps%file_format) 406 case(PSEUDO_FORMAT_PSF) 407 call ps_psf_get_eigen(ps_psf, eigen) 408 call ps_grid_load(ps, ps_psf%ps_grid) 409 call ps_psf_end(ps_psf) 410 case(PSEUDO_FORMAT_CPI) 411 call ps_grid_load(ps, ps_cpi%ps_grid) 412 call ps_cpi_end(ps_cpi) 413 case(PSEUDO_FORMAT_FHI) 414 call ps_grid_load(ps, ps_fhi%ps_grid) 415 call ps_fhi_end(ps_fhi) 416 case(PSEUDO_FORMAT_HGH) 417 call hgh_get_eigen(ps_hgh, eigen) 418 SAFE_ALLOCATE(ps%k (0:ps%lmax, 1:ps%kbc, 1:ps%kbc)) 419 call hgh_load(ps, ps_hgh) 420 call hgh_end(ps_hgh) 421 case(PSEUDO_FORMAT_QSO, PSEUDO_FORMAT_UPF1, PSEUDO_FORMAT_UPF2, PSEUDO_FORMAT_PSML, PSEUDO_FORMAT_PSP8) 422 call ps_xml_load(ps, ps_xml) 423 call ps_xml_end(ps_xml) 424 end select 425 426 if(ps_has_density(ps)) then 427 do is = 1, ps%ispin 428 call spline_der(ps%density(is), ps%density_der(is)) 429 end do 430 end if 431 432 if(ps_has_nlcc(ps)) then 433 call spline_der(ps%core, ps%core_der) 434 end if 435 436 call ps_check_bound(ps, eigen) 437 438 ps%has_long_range = .true. 439 ps%is_separated = .false. 440 441 call ps_info(ps, filename) 442 443 SAFE_DEALLOCATE_A(eigen) 444 445 POP_SUB(ps_init) 446 end subroutine ps_init 447 448 !------------------------------------------------------------------------ 449 450 subroutine ps_info(ps, filename) 451 type(ps_t), intent(in) :: ps 452 character(len=*), intent(in) :: filename 453 454 call messages_write(" Species '"//trim(ps%label)//"'", new_line = .true.) 455 call messages_write(" type : pseudopotential", new_line = .true.) 456 call messages_write(" file : '"//trim(filename)//"'") 457 call messages_info() 458 459 call messages_write(" file format :") 460 select case(ps%file_format) 461 case(PSEUDO_FORMAT_UPF1) 462 call messages_write(" UPF1") 463 case(PSEUDO_FORMAT_UPF2) 464 call messages_write(" UPF2") 465 case(PSEUDO_FORMAT_QSO) 466 call messages_write(" QSO") 467 case(PSEUDO_FORMAT_PSML) 468 call messages_write(" PSML") 469 case(PSEUDO_FORMAT_PSP8) 470 call messages_write(" PSP8") 471 case(PSEUDO_FORMAT_PSF) 472 call messages_write(" PSF") 473 case(PSEUDO_FORMAT_CPI) 474 call messages_write(" CPI") 475 case(PSEUDO_FORMAT_FHI) 476 call messages_write(" FHI") 477 case(PSEUDO_FORMAT_HGH) 478 call messages_write(" HGH") 479 end select 480 call messages_new_line() 481 482 call messages_write(" valence charge :") 483 call messages_write(ps%z_val, align_left = .true., fmt = '(f4.1)') 484 call messages_info() 485 486 call messages_write(" atomic number :") 487 call messages_write(nint(ps%z), fmt = '(i4)') 488 call messages_info() 489 490 call messages_write(" form on file :") 491 select case(ps%pseudo_type) 492 case(PSEUDO_TYPE_ULTRASOFT) 493 call messages_write(" ultrasoft") 494 case(PSEUDO_TYPE_SEMILOCAL) 495 call messages_write(" semilocal") 496 case(PSEUDO_TYPE_KLEINMAN_BYLANDER) 497 call messages_write(" kleinman-bylander") 498 case(PSEUDO_TYPE_PAW) 499 call messages_write(" paw") 500 end select 501 call messages_info() 502 503 if(ps%pseudo_type == PSEUDO_TYPE_SEMILOCAL) then 504 call messages_write(" orbital origin :") 505 select case(ps%file_format) 506 case(PSEUDO_FORMAT_PSF) 507 call messages_write(" calculated"); 508 case default 509 call messages_write(" from file"); 510 end select 511 call messages_info() 512 end if 513 514 call messages_write(" lmax :") 515 call messages_write(ps%lmax, fmt = '(i2)') 516 call messages_info() 517 518 call messages_write(" llocal :") 519 if(ps%llocal >= 0) then 520 call messages_write(ps%llocal, fmt = '(i2)') 521 else 522 call messages_write(ps%llocal, fmt = '(i3)') 523 end if 524 call messages_info() 525 526 call messages_write(" projectors per l :") 527 call messages_write(ps%kbc, fmt = '(i2)') 528 call messages_info() 529 530 call messages_write(" total projectors :") 531 if(ps%llocal < 0) then 532 call messages_write(ps%kbc*(ps%lmax + 1), fmt = '(i2)') 533 else 534 call messages_write(ps%kbc*ps%lmax, fmt = '(i2)') 535 end if 536 call messages_info() 537 538 if(ps%local) then 539 call messages_write(" application form : local") 540 else 541 call messages_write(" application form : kleinman-bylander") 542 end if 543 call messages_info() 544 545 call messages_write(" orbitals :") 546 call messages_write(ps_niwfs(ps), fmt='(i3)') 547 call messages_info() 548 call messages_write(" bound orbitals :") 549 call messages_write(ps_bound_niwfs(ps), fmt='(i3)') 550 call messages_info() 551 552 call messages_info() 553 554 end subroutine ps_info 555 556 557 ! --------------------------------------------------------- 558 !> separate the local potential into (soft) long-ranged and (hard) short-ranged parts 559 subroutine ps_separate(ps) 560 type(ps_t), intent(inout) :: ps 561 562 FLOAT, allocatable :: vsr(:), vlr(:), nlr(:) 563 FLOAT :: r, exp_arg 564 integer :: ii 565 566 PUSH_SUB(ps_separate) 567 568 ASSERT(ps%g%nrval > 0) 569 570 SAFE_ALLOCATE(vsr(1:ps%g%nrval)) 571 SAFE_ALLOCATE(vlr(1:ps%g%nrval)) 572 SAFE_ALLOCATE(nlr(1:ps%g%nrval)) 573 574 ps%sigma_erf = CNST(0.625) ! This is hard-coded to a reasonable value 575 576 vlr(1) = -ps%z_val*M_TWO/(sqrt(M_TWO*M_PI)*ps%sigma_erf) 577 578 do ii = 1, ps%g%nrval 579 r = ps%g%rofi(ii) 580 if (ii > 1) then 581 vlr(ii) = -ps%z_val*loct_erf(r/(ps%sigma_erf*sqrt(M_TWO)))/r 582 end if 583 vsr(ii) = spline_eval(ps%vl, r) - vlr(ii) 584 585 exp_arg = -M_HALF*r**2/ps%sigma_erf**2 586 if(exp_arg > CNST(-100)) then 587 nlr(ii) = -ps%z_val*M_ONE/(ps%sigma_erf*sqrt(M_TWO*M_PI))**3*exp(exp_arg) 588 else 589 nlr(ii) = M_ZERO 590 end if 591 end do 592 593 call spline_init(ps%vlr) 594 call spline_fit(ps%g%nrval, ps%g%rofi, vlr, ps%vlr) 595 596 call spline_init(ps%vlr_sq) 597 call spline_fit(ps%g%nrval, ps%g%r2ofi, vlr, ps%vlr_sq) 598 599 call spline_init(ps%nlr) 600 call spline_fit(ps%g%nrval, ps%g%rofi, nlr, ps%nlr) 601 602 !overwrite vl 603 call spline_end(ps%vl) 604 call spline_init(ps%vl) 605 call spline_fit(ps%g%nrval, ps%g%rofi, vsr, ps%vl) 606 607 SAFE_DEALLOCATE_A(vsr) 608 SAFE_DEALLOCATE_A(vlr) 609 SAFE_DEALLOCATE_A(nlr) 610 611 ps%is_separated = .true. 612 613 POP_SUB(ps_separate) 614 end subroutine ps_separate 615 616 617 ! --------------------------------------------------------- 618 subroutine ps_getradius(ps) 619 type(ps_t), intent(inout) :: ps 620 integer :: l, j 621 622 PUSH_SUB(ps_getradius) 623 624 ps%rc_max = CNST(0.0) 625 626 do l = 0, ps%lmax 627 if(l == ps%llocal) cycle 628 do j = 1, ps%kbc 629 ps%rc_max = max(ps%rc_max, spline_cutoff_radius(ps%kb(l, j), ps%projectors_sphere_threshold)) 630 end do 631 end do 632 633 POP_SUB(ps_getradius) 634 end subroutine ps_getradius 635 636 637 ! --------------------------------------------------------- 638 subroutine ps_derivatives(ps) 639 type(ps_t), intent(inout) :: ps 640 integer :: l, j 641 642 PUSH_SUB(ps_derivatives) 643 644 do l = 0, ps%lmax 645 do j = 1, ps%kbc 646 call spline_der(ps%kb(l, j), ps%dkb(l, j)) 647 end do 648 end do 649 650 POP_SUB(ps_derivatives) 651 end subroutine ps_derivatives 652 653 654 ! --------------------------------------------------------- 655 subroutine ps_filter(ps, filter, gmax) 656 type(ps_t), intent(inout) :: ps 657 integer, intent(in) :: filter 658 FLOAT, intent(in) :: gmax 659 660 integer :: l, k, ispin 661 type(profile_t), save:: prof 662 663 FLOAT :: alpha, beta_fs, rmax, rcut, gamma, beta_rs 664 665 PUSH_SUB(ps_filter) 666 call profiling_in(prof, "PS_FILTER") 667 668 select case(filter) 669 case(PS_FILTER_NONE) 670 671 case(PS_FILTER_TS) 672 alpha = CNST(1.1) 673 gamma = CNST(2.0) 674 675 rmax = spline_cutoff_radius(ps%vl, ps%projectors_sphere_threshold) 676 call spline_filter_mask(ps%vl, max(0, ps%llocal), rmax, gmax, alpha, gamma) 677 do l = 0, ps%lmax 678 if(l == ps%llocal) cycle 679 do k = 1, ps%kbc 680 call spline_filter_mask(ps%kb(l, k), l, ps%rc_max, gmax, alpha, gamma) 681 end do 682 end do 683 684 if(ps_has_nlcc(ps)) then 685 rmax = spline_cutoff_radius(ps%core, ps%projectors_sphere_threshold) 686 call spline_filter_mask(ps%core, 0, rmax, gmax, alpha, gamma) 687 end if 688 689 if(ps_has_density(ps)) then 690 do ispin = 1, ps%ispin 691 if(abs(spline_integral(ps%density(ispin))) > CNST(1.0e-12)) then 692 rmax = spline_cutoff_radius(ps%density(ispin), ps%projectors_sphere_threshold) 693 call spline_filter_mask(ps%density(ispin), 0, rmax, gmax, alpha, gamma) 694 call spline_force_pos(ps%density(ispin)) 695 end if 696 697 if(abs(spline_integral(ps%density_der(ispin))) > CNST(1.0e-12)) then 698 rmax = spline_cutoff_radius(ps%density_der(ispin), ps%projectors_sphere_threshold) 699 call spline_filter_mask(ps%density_der(ispin), 0, rmax, gmax, alpha, gamma) 700 end if 701 end do 702 end if 703 704 case(PS_FILTER_BSB) 705 alpha = CNST(0.7) ! The original was M_FOUR/CNST(7.0) 706 beta_fs = CNST(18.0) 707 rcut = CNST(2.5) 708 beta_rs = CNST(0.4) 709 710 call spline_filter_bessel(ps%vl, ps%llocal, gmax, alpha, beta_fs, rcut, beta_rs) 711 do l = 0, ps%lmax 712 if(l == ps%llocal) cycle 713 do k = 1, ps%kbc 714 call spline_filter_bessel(ps%kb(l, k), l, gmax, alpha, beta_fs, rcut, beta_rs) 715 end do 716 end do 717 718 if(ps_has_nlcc(ps)) then 719 call spline_filter_bessel(ps%core, 0, gmax, alpha, beta_fs, rcut, beta_rs) 720 end if 721 722 if(ps_has_density(ps)) then 723 do ispin = 1, ps%ispin 724 call spline_filter_bessel(ps%density(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs) 725 call spline_force_pos(ps%density(ispin)) 726 call spline_filter_bessel(ps%density_der(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs) 727 end do 728 end if 729 730 end select 731 732 call profiling_out(prof) 733 POP_SUB(ps_filter) 734 end subroutine ps_filter 735 736 ! --------------------------------------------------------- 737 subroutine ps_check_bound(ps, eigen) 738 type(ps_t), intent(inout) :: ps 739 FLOAT, intent(in) :: eigen(:,:) 740 741 integer :: i, is, ir 742 FLOAT :: ur1, ur2 743 744 PUSH_SUB(ps_check_bound) 745 746 ! Unbound states have positive eigenvalues 747 where(eigen > M_ZERO) 748 ps%bound = .false. 749 elsewhere 750 ps%bound = .true. 751 end where 752 753 ! We might not have information about the eigenvalues, so we need to check the wavefunctions 754 do i = 1, ps%conf%p 755 do is = 1, ps%ispin 756 if (.not. ps%bound(i, is)) cycle 757 758 do ir = ps%g%nrval, 3, -1 759 ! First we look for the outmost value that is not zero 760 if (abs(spline_eval(ps%ur(i, is), ps%g%rofi(ir))*ps%g%rofi(ir)) > M_ZERO) then 761 ! Usually bound states have exponentially decaying wavefunctions, 762 ! while unbound states have exponentially diverging 763 ! wavefunctions. Therefore we check if the wavefunctions 764 ! value is increasing with increasing radius. The fact 765 ! that we do not use the wavefunctions outmost value that 766 ! is not zero is on purpose, as some pseudopotential 767 ! generators do funny things with that point. 768 ur1 = spline_eval(ps%ur(i, is), ps%g%rofi(ir-2))*ps%g%rofi(ir-2) 769 ur2 = spline_eval(ps%ur(i, is), ps%g%rofi(ir-1))*ps%g%rofi(ir-1) 770 if ((ur1*ur2 > M_ZERO) .and. (abs(ur2) > abs(ur1))) ps%bound(i, is) = .false. 771 exit 772 end if 773 end do 774 end do 775 end do 776 777 POP_SUB(ps_check_bound) 778 end subroutine ps_check_bound 779 780 781 ! --------------------------------------------------------- 782 subroutine ps_debug(ps, dir, namespace) 783 type(ps_t), intent(in) :: ps 784 character(len=*), intent(in) :: dir 785 type(namespace_t), intent(in) :: namespace 786 787 ! We will plot also some Fourier transforms. 788 type(spline_t), allocatable :: fw(:, :) 789 FLOAT, parameter :: gmax = CNST(40.0) 790 791 integer :: iunit 792 integer :: j, k, l 793 794 PUSH_SUB(ps_debug) 795 796 ! A text file with some basic data. 797 iunit = io_open(trim(dir)//'/pseudo-info', namespace, action='write') 798 write(iunit,'(a,/)') ps%label 799 write(iunit,'(a,a,/)') 'Format : ', ps_name(ps%file_format) 800 write(iunit,'(a,f6.3)') 'z : ', ps%z 801 write(iunit,'(a,f6.3,/)') 'zval : ', ps%z_val 802 write(iunit,'(a,i4)') 'lmax : ', ps%lmax 803 write(iunit,'(a,i4)') 'lloc : ', ps%llocal 804 write(iunit,'(a,i4,/)') 'kbc : ', ps%kbc 805 write(iunit,'(a,f9.5,/)') 'rcmax : ', ps%rc_max 806 write(iunit,'(a,/)') 'h matrix:' 807 do l = 0, ps%lmax 808 do k = 1, ps%kbc 809 write(iunit,'(10f9.5)') (ps%h(l, k, j), j = 1, ps%kbc) 810 end do 811 end do 812 if(associated(ps%k)) then 813 write(iunit,'(/,a,/)') 'k matrix:' 814 do l = 0, ps%lmax 815 do k = 1, ps%kbc 816 write(iunit,'(10f9.5)') (ps%k(l, k, j), j = 1, ps%kbc) 817 end do 818 end do 819 end if 820 821 write(iunit,'(/,a)') 'orbitals:' 822 do j = 1, ps%conf%p 823 write(iunit,'(1x,a,i2,3x,a,i2,3x,a,f5.1,3x,a,l1)') 'n = ', ps%conf%n(j), 'l = ', ps%conf%l(j), & 824 'j = ', ps%conf%j(j), 'bound = ', all(ps%bound(j,:)) 825 end do 826 827 828 call io_close(iunit) 829 830 ! Local part of the pseudopotential 831 iunit = io_open(trim(dir)//'/local', namespace, action='write') 832 call spline_print(ps%vl, iunit) 833 call io_close(iunit) 834 835 ! Local part of the pseudopotential 836 iunit = io_open(trim(dir)//'/local_long_range', namespace, action='write') 837 call spline_print(ps%vlr, iunit) 838 call io_close(iunit) 839 840 ! Local part of the pseudopotential 841 iunit = io_open(trim(dir)//'/local_long_range_density', namespace, action='write') 842 call spline_print(ps%nlr, iunit) 843 call io_close(iunit) 844 845 ! Fourier transform of the local part 846 iunit = io_open(trim(dir)//'/local_ft', namespace, action='write') 847 SAFE_ALLOCATE(fw(1:1, 1:1)) 848 call spline_init(fw(1, 1)) 849 call spline_3dft(ps%vl, fw(1, 1), gmax = gmax) 850 call spline_print(fw(1, 1), iunit) 851 call spline_end(fw(1, 1)) 852 SAFE_DEALLOCATE_A(fw) 853 call io_close(iunit) 854 855 ! Kleinman-Bylander projectors 856 iunit = io_open(trim(dir)//'/nonlocal', namespace, action='write') 857 call spline_print(ps%kb, iunit) 858 call io_close(iunit) 859 860 iunit = io_open(trim(dir)//'/nonlocal_derivative', namespace, action='write') 861 call spline_print(ps%dkb, iunit) 862 call io_close(iunit) 863 864 iunit = io_open(trim(dir)//'/nonlocal_ft', namespace, action='write') 865 SAFE_ALLOCATE(fw(0:ps%lmax, 1:ps%kbc)) 866 call spline_init(fw) 867 do k = 0, ps%lmax 868 do j = 1, ps%kbc 869 call spline_3dft(ps%kb(k, j), fw(k, j), gmax = gmax) 870 end do 871 end do 872 call spline_print(fw, iunit) 873 call spline_end(fw) 874 SAFE_DEALLOCATE_A(fw) 875 call io_close(iunit) 876 877 ! Pseudo-wavefunctions 878 iunit = io_open(trim(dir)//'/wavefunctions', namespace, action='write') 879 call spline_print(ps%ur, iunit) 880 call io_close(iunit) 881 882 ! Density 883 if (ps%has_density) then 884 iunit = io_open(trim(dir)//'/density', namespace, action='write') 885 call spline_print(ps%density, iunit) 886 call io_close(iunit) 887 888 iunit = io_open(trim(dir)//'/density_derivative', namespace, action='write') 889 call spline_print(ps%density_der, iunit) 890 call io_close(iunit) 891 end if 892 893 ! Non-linear core-corrections 894 if(ps_has_nlcc(ps)) then 895 iunit = io_open(trim(dir)//'/nlcc', namespace, action='write') 896 call spline_print(ps%core, iunit) 897 call io_close(iunit) 898 end if 899 900 POP_SUB(ps_debug) 901 end subroutine ps_debug 902 903 904 ! --------------------------------------------------------- 905 subroutine ps_end(ps) 906 type(ps_t), intent(inout) :: ps 907 908 if(.not. associated(ps%kb)) return 909 910 PUSH_SUB(ps_end) 911 912 if(ps%is_separated) then 913 call spline_end(ps%vlr) 914 call spline_end(ps%vlr_sq) 915 call spline_end(ps%nlr) 916 end if 917 918 call spline_end(ps%kb) 919 call spline_end(ps%dkb) 920 call spline_end(ps%ur) 921 call spline_end(ps%ur_sq) 922 923 call spline_end(ps%vl) 924 call spline_end(ps%core) 925 call spline_end(ps%core_der) 926 927 if(associated(ps%density)) call spline_end(ps%density) 928 if(associated(ps%density_der)) call spline_end(ps%density_der) 929 930 call logrid_end(ps%g) 931 932 SAFE_DEALLOCATE_P(ps%kb) 933 SAFE_DEALLOCATE_P(ps%dkb) 934 SAFE_DEALLOCATE_P(ps%ur) 935 SAFE_DEALLOCATE_P(ps%ur_sq) 936 SAFE_DEALLOCATE_A(ps%bound) 937 SAFE_DEALLOCATE_P(ps%h) 938 SAFE_DEALLOCATE_P(ps%k) 939 SAFE_DEALLOCATE_P(ps%density) 940 SAFE_DEALLOCATE_P(ps%density_der) 941 942 POP_SUB(ps_end) 943 end subroutine ps_end 944 945 946 ! --------------------------------------------------------- 947 subroutine hgh_load(ps, ps_hgh) 948 type(ps_t), intent(inout) :: ps 949 type(ps_hgh_t), intent(inout) :: ps_hgh 950 951 integer :: l, ll 952 FLOAT :: x 953 954 PUSH_SUB(hgh_load) 955 956 ! Fixes some components of ps 957 ps%z_val = ps_hgh%z_val 958 ps%nlcc = .false. 959 if(ps%lmax>=0) then 960 ps%rc_max = CNST(1.1) * maxval(ps_hgh%kbr(0:ps%lmax)) ! Increase a little. 961 else 962 ps%rc_max = M_ZERO 963 end if 964 ps%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) 965 ps%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) 966 967 ! Fixes the occupations 968 if(ps%ispin == 2) then 969 do l = 1, ps%conf%p 970 ll = ps%conf%l(l) 971 x = ps%conf%occ(l, 1) 972 ps%conf%occ(l, 1) = min(x, TOFLOAT(2*ll+1)) 973 ps%conf%occ(l, 2) = x - ps%conf%occ(l, 1) 974 end do 975 end if 976 977 ! now we fit the splines 978 call get_splines() 979 980 POP_SUB(hgh_load) 981 982 contains 983 984 ! --------------------------------------------------------- 985 subroutine get_splines() 986 integer :: l, is, nrc, j, ip 987 FLOAT, allocatable :: hato(:), dens(:) 988 989 PUSH_SUB(hgh_load.get_splines) 990 991 SAFE_ALLOCATE(hato(1:ps_hgh%g%nrval)) 992 SAFE_ALLOCATE(dens(1:ps_hgh%g%nrval)) 993 994 ! Interpolate the KB-projection functions 995 do l = 0, ps_hgh%l_max 996 do j = 1, 3 997 hato = M_ZERO 998 nrc = nint(log(ps_hgh%kbr(l)/ps_hgh%g%b + M_ONE)/ps_hgh%g%a) + 1 999 hato(1:nrc) = ps_hgh%kb(1:nrc, l, j) 1000 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%kb(l, j)) 1001 end do 1002 end do 1003 1004 ! Now the part corresponding to the local pseudopotential 1005 ! where the asymptotic part is subtracted 1006 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, ps_hgh%vlocal, ps%vl) 1007 1008 ! Define the table for the pseudo-wavefunction components (using splines) 1009 ! with a correct normalization function 1010 do is = 1, ps%ispin 1011 dens = CNST(0.0) 1012 do l = 1, ps%conf%p 1013 hato(2:ps_hgh%g%nrval) = ps_hgh%rphi(2:ps_hgh%g%nrval, l)/ps_hgh%g%rofi(2:ps_hgh%g%nrval) 1014 hato(1) = hato(2) 1015 1016 do ip = 1, ps_hgh%g%nrval 1017 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(M_FOUR*M_PI) 1018 end do 1019 1020 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%ur(l, is)) 1021 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%r2ofi, hato, ps%ur_sq(l, is)) 1022 end do 1023 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, dens, ps%density(is)) 1024 end do 1025 1026 SAFE_DEALLOCATE_A(hato) 1027 SAFE_DEALLOCATE_A(dens) 1028 1029 POP_SUB(hgh_load.get_splines) 1030 end subroutine get_splines 1031 end subroutine hgh_load 1032 1033 1034 ! --------------------------------------------------------- 1035 subroutine ps_grid_load(ps, ps_grid) 1036 type(ps_t), intent(inout) :: ps 1037 type(ps_in_grid_t), intent(in) :: ps_grid 1038 1039 PUSH_SUB(ps_grid_load) 1040 1041 ! Fixes some components of ps, read in ps_grid 1042 ps%z_val = ps_grid%zval 1043 1044 ps%nlcc = ps_grid%core_corrections 1045 1046 ps%h(0:ps%lmax, 1, 1) = ps_grid%dkbcos(1:ps%lmax+1) 1047 1048 ! Increasing radius a little, just in case. 1049 ! I have hard-coded a larger increase of the cutoff for the filtering. 1050 ps%rc_max = maxval(ps_grid%kb_radius(1:ps%lmax+1)) * CNST(1.5) 1051 1052 ! now we fit the splines 1053 call get_splines(ps_grid%g) 1054 1055 ! Passes from Rydbergs to Hartrees. 1056 ps%h(0:ps%lmax,:,:) = ps%h(0:ps%lmax,:,:) / M_TWO 1057 1058 POP_SUB(ps_grid_load) 1059 1060 contains 1061 1062 subroutine get_splines(g) 1063 type(logrid_t), intent(in) :: g 1064 1065 FLOAT, allocatable :: hato(:), dens(:) 1066 integer :: is, l, ir, nrc, ip 1067 1068 PUSH_SUB(ps_grid_load.get_splines) 1069 1070 SAFE_ALLOCATE(hato(1:g%nrval)) 1071 SAFE_ALLOCATE(dens(1:g%nrval)) 1072 1073 ! the wavefunctions 1074 do is = 1, ps%ispin 1075 1076 dens = CNST(0.0) 1077 1078 do l = 1, ps_grid%no_l_channels 1079 hato(2:) = ps_grid%rphi(2:, l, 1+is)/g%rofi(2:) 1080 hato(1) = first_point_extrapolate(g%rofi, hato) 1081 1082 do ip = 1, g%nrval 1083 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(M_FOUR*M_PI) 1084 end do 1085 1086 call spline_fit(g%nrval, g%rofi, hato, ps%ur(l, is)) 1087 call spline_fit(g%nrval, g%r2ofi, hato, ps%ur_sq(l, is)) 1088 1089 end do 1090 1091 call spline_fit(g%nrval, g%rofi, dens, ps%density(is)) 1092 end do 1093 1094 1095 ! the Kleinman-Bylander projectors 1096 do l = 1, ps%lmax+1 1097 nrc = logrid_index(g, ps_grid%kb_radius(l)) + 1 1098 hato(1:nrc) = ps_grid%KB(1:nrc, l) 1099 hato(nrc+1:g%nrval) = M_ZERO 1100 1101 call spline_fit(g%nrval, g%rofi, hato, ps%kb(l-1, 1)) 1102 end do 1103 1104 ! Now the part corresponding to the local pseudopotential 1105 ! where the asymptotic part is subtracted 1106 hato(:) = ps_grid%vlocal(:)/M_TWO 1107 call spline_fit(g%nrval, g%rofi, hato, ps%vl) 1108 1109 if(ps_grid%core_corrections) then 1110 ! find cutoff radius 1111 hato(2:) = ps_grid%chcore(2:)/(M_FOUR*M_PI*g%rofi(2:)**2) 1112 1113 do ir = g%nrval-1, 2, -1 1114 if(hato(ir) > eps) then 1115 nrc = ir + 1 1116 exit 1117 end if 1118 end do 1119 1120 hato(nrc:g%nrval) = M_ZERO 1121 hato(1) = first_point_extrapolate(g%rofi, hato) 1122 1123 call spline_fit(g%nrval, g%rofi, hato, ps%core) 1124 end if 1125 1126 SAFE_DEALLOCATE_A(hato) 1127 SAFE_DEALLOCATE_A(dens) 1128 1129 POP_SUB(ps_grid_load.get_splines) 1130 end subroutine get_splines 1131 end subroutine ps_grid_load 1132 1133 ! --------------------------------------------------------- 1134 1135 subroutine ps_xml_load(ps, ps_xml) 1136 type(ps_t), intent(inout) :: ps 1137 type(ps_xml_t), intent(in) :: ps_xml 1138 1139 integer :: ll, ip, is, ic, jc, ir, nrc, ii 1140 FLOAT :: rr, kbcos, kbnorm, dnrm, avgv, volume_element 1141 FLOAT, allocatable :: vlocal(:), kbprojector(:), wavefunction(:), nlcc_density(:), dens(:) 1142 integer, allocatable :: cmap(:, :) 1143 FLOAT, allocatable :: matrix(:, :), eigenvalues(:) 1144 1145 PUSH_SUB(ps_xml_load) 1146 1147 ps%hamann = (ps_xml%kleinman_bylander .and. ps_xml%nchannels == 2 .and. ps_xml%llocal == -1) 1148 1149 ps%nlcc = ps_xml%nlcc 1150 1151 ps%z_val = ps_xml%valence_charge 1152 1153 ! the local potential 1154 SAFE_ALLOCATE(vlocal(1:ps%g%nrval)) 1155 1156 do ip = 1, ps%g%nrval 1157 rr = ps_xml%grid(ip) 1158 if(ip <= ps_xml%grid_size) then 1159 vlocal(ip) = ps_xml%potential(ip, ps%llocal) 1160 else 1161 vlocal(ip) = -ps_xml%valence_charge/rr 1162 end if 1163 end do 1164 1165 call spline_fit(ps%g%nrval, ps%g%rofi, vlocal, ps%vl) 1166 1167 SAFE_DEALLOCATE_A(vlocal) 1168 1169 SAFE_ALLOCATE(kbprojector(1:ps%g%nrval)) 1170 SAFE_ALLOCATE(wavefunction(1:ps%g%nrval)) 1171 1172 kbprojector = CNST(0.0) 1173 wavefunction = CNST(0.0) 1174 1175 ! the projectors and the orbitals 1176 if(ps_xml%kleinman_bylander) then 1177 1178 SAFE_ALLOCATE(cmap(0:ps_xml%lmax, 1:ps_xml%nchannels)) 1179 1180 ! the order of the channels is determined by spin orbit and the j value 1181 do ll = 0, ps_xml%lmax 1182 do ic = 1, ps_xml%nchannels 1183 cmap(ll, ic) = ic 1184 1185 if(ll == 0) cycle 1186 if(ll == ps_xml%llocal) cycle 1187 if(.not. pseudo_has_total_angular_momentum(ps_xml%pseudo)) cycle 1188 1189 ASSERT(ps_xml%nchannels == 2) 1190 if(pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll - 1) then 1191 ! this is Octopus convention 1192 cmap(ll, ic) = 2 1193 else 1194 ASSERT(pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll + 1) 1195 cmap(ll, ic) = 1 1196 end if 1197 1198 end do 1199 1200 ! check that all numbers are present for each l 1201 ASSERT(sum(cmap(ll, 1:ps_xml%nchannels)) == (ps_xml%nchannels + 1)*ps_xml%nchannels/2) 1202 end do 1203 1204 ASSERT(all(cmap >= 0 .and. cmap <= ps_xml%nchannels)) 1205 1206 SAFE_ALLOCATE(matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels)) 1207 SAFE_ALLOCATE(eigenvalues(1:ps_xml%nchannels)) 1208 1209 ps%h = CNST(0.0) 1210 1211 1212 if(pseudo_nprojectors(ps_xml%pseudo) > 0) then 1213 do ll = 0, ps_xml%lmax 1214 1215 if (is_diagonal(ps_xml%nchannels, ps_xml%dij(ll, :, :)) .or. & 1216 pseudo_has_total_angular_momentum(ps_xml%pseudo)) then 1217 matrix = CNST(0.0) 1218 do ic = 1, ps_xml%nchannels 1219 eigenvalues(ic) = ps_xml%dij(ll, ic, ic) 1220 matrix(ic, ic) = CNST(1.0) 1221 end do 1222 else 1223 ! diagonalize the coefficient matrix 1224 matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels) = ps_xml%dij(ll, 1:ps_xml%nchannels, 1:ps_xml%nchannels) 1225 call lalg_eigensolve(ps_xml%nchannels, matrix, eigenvalues) 1226 end if 1227 1228 do ic = 1, ps_xml%nchannels 1229 1230 do ip = 1, ps%g%nrval 1231 kbprojector(ip) = M_ZERO 1232 if(ip <= ps_xml%grid_size) then 1233 do jc = 1, ps_xml%nchannels 1234 kbprojector(ip) = kbprojector(ip) + matrix(jc, ic)*ps_xml%projector(ip, ll, jc) 1235 end do 1236 end if 1237 end do 1238 1239 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, cmap(ll, ic))) 1240 1241 ps%h(ll, cmap(ll, ic), cmap(ll, ic)) = eigenvalues(ic) 1242 1243 end do 1244 end do 1245 end if 1246 1247 SAFE_DEALLOCATE_A(matrix) 1248 SAFE_DEALLOCATE_A(eigenvalues) 1249 1250 ps%conf%p = ps_xml%nwavefunctions 1251 1252 do ii = 1, ps_xml%nwavefunctions 1253 1254 ps%conf%n(ii) = ps_xml%wf_n(ii) 1255 ps%conf%l(ii) = ps_xml%wf_l(ii) 1256 1257 if(ps%ispin == 2) then 1258 ps%conf%occ(ii, 1) = min(ps_xml%wf_occ(ii), CNST(2.0)*ps_xml%wf_l(ii) + CNST(1.0)) 1259 ps%conf%occ(ii, 2) = ps_xml%wf_occ(ii) - ps%conf%occ(ii, 1) 1260 else 1261 ps%conf%occ(ii, 1) = ps_xml%wf_occ(ii) 1262 end if 1263 1264 ps%conf%j(ii) = M_ZERO 1265 if(pseudo_has_total_angular_momentum(ps_xml%pseudo)) then 1266 ps%conf%j(ii) = M_HALF*pseudo_wavefunction_2j(ps_xml%pseudo, ii) 1267 end if 1268 1269 do ip = 1, ps%g%nrval 1270 if(ip <= ps_xml%grid_size) then 1271 wavefunction(ip) = ps_xml%wavefunction(ip, ii) 1272 else 1273 wavefunction(ip) = CNST(0.0) 1274 end if 1275 end do 1276 1277 do is = 1, ps%ispin 1278 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ii, is)) 1279 call spline_fit(ps%g%nrval, ps%g%r2ofi, wavefunction, ps%ur_sq(ii, is)) 1280 end do 1281 1282 end do 1283 1284 SAFE_DEALLOCATE_A(cmap) 1285 1286 else 1287 1288 do ll = 0, ps_xml%lmax 1289 ! we need to build the KB projectors 1290 ! the procedure was copied from ps_in_grid.F90 (r12967) 1291 dnrm = M_ZERO 1292 avgv = M_ZERO 1293 do ip = 1, ps_xml%grid_size 1294 rr = ps_xml%grid(ip) 1295 volume_element = rr**2*ps_xml%weights(ip) 1296 kbprojector(ip) = (ps_xml%potential(ip, ll) - ps_xml%potential(ip, ps%llocal))*ps_xml%wavefunction(ip, ll) 1297 dnrm = dnrm + kbprojector(ip)**2*volume_element 1298 avgv = avgv + kbprojector(ip)*ps_xml%wavefunction(ip, ll)*volume_element 1299 end do 1300 1301 kbcos = dnrm/(avgv + CNST(1.0e-20)) 1302 kbnorm = M_ONE/(sqrt(dnrm) + CNST(1.0e-20)) 1303 1304 if(ll /= ps%llocal) then 1305 ps%h(ll, 1, 1) = kbcos 1306 kbprojector = kbprojector*kbnorm 1307 else 1308 ps%h(ll, 1, 1) = CNST(0.0) 1309 end if 1310 1311 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, 1)) 1312 1313 ! wavefunctions, for the moment we pad them with zero 1314 do ip = 1, ps%g%nrval 1315 if(ip <= ps_xml%grid_size) then 1316 wavefunction(ip) = ps_xml%wavefunction(ip, ll) 1317 else 1318 wavefunction(ip) = CNST(0.0) 1319 end if 1320 end do 1321 1322 do is = 1, ps%ispin 1323 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ll + 1, is)) 1324 call spline_fit(ps%g%nrval, ps%g%r2ofi, wavefunction, ps%ur_sq(ll + 1, is)) 1325 end do 1326 end do 1327 1328 end if 1329 1330 ps%has_density = ps_xml%has_density 1331 1332 if(ps_has_density(ps)) then 1333 1334 SAFE_ALLOCATE(dens(1:ps%g%nrval)) 1335 1336 dens(1:ps_xml%grid_size) = ps_xml%density(1:ps_xml%grid_size)/ps%ispin 1337 dens(ps_xml%grid_size + 1:ps%g%nrval) = CNST(0.0) 1338 1339 do is = 1, ps%ispin 1340 call spline_fit(ps%g%nrval, ps%g%rofi, dens, ps%density(is)) 1341 end do 1342 1343 SAFE_DEALLOCATE_A(dens) 1344 end if 1345 1346 !Non-linear core-corrections 1347 if(ps_xml%nlcc) then 1348 1349 SAFE_ALLOCATE(nlcc_density(1:ps%g%nrval)) 1350 1351 nlcc_density(1:ps_xml%grid_size) = ps_xml%nlcc_density(1:ps_xml%grid_size) 1352 1353 ! find cutoff radius 1354 do ir = ps_xml%grid_size - 1, 1, -1 1355 if(nlcc_density(ir) > eps) then 1356 nrc = ir + 1 1357 exit 1358 end if 1359 end do 1360 1361 nlcc_density(nrc:ps%g%nrval) = M_ZERO 1362 1363 call spline_fit(ps%g%nrval, ps%g%rofi, nlcc_density, ps%core) 1364 1365 SAFE_DEALLOCATE_A(nlcc_density) 1366 end if 1367 1368 call ps_getradius(ps) 1369 1370 SAFE_DEALLOCATE_A(kbprojector) 1371 SAFE_DEALLOCATE_A(wavefunction) 1372 1373 POP_SUB(ps_xml_load) 1374 end subroutine ps_xml_load 1375 1376 ! --------------------------------------------------------- 1377 1378 logical function is_diagonal(dim, matrix) 1379 integer, intent(in) :: dim 1380 FLOAT, intent(in) :: matrix(:, :) 1381 1382 integer :: ii, jj 1383 1384 is_diagonal = .true. 1385 do ii = 1, dim 1386 do jj = 1, dim 1387 if(ii == jj) cycle 1388 if(abs(matrix(ii, jj)) > CNST(1e10)) is_diagonal = .false. 1389 end do 1390 end do 1391 1392 end function is_diagonal 1393 1394 ! --------------------------------------------------------- 1395 !> Returns the number of atomic orbitals taking into account then m quantum number multiplicity 1396 pure integer function ps_niwfs(ps) 1397 type(ps_t), intent(in) :: ps 1398 1399 integer :: i, l 1400 1401 ps_niwfs = 0 1402 do i = 1, ps%conf%p 1403 l = ps%conf%l(i) 1404 ps_niwfs = ps_niwfs + (2*l+1) 1405 end do 1406 1407 end function ps_niwfs 1408 1409 ! --------------------------------------------------------- 1410 !> Returns the number of bound atomic orbitals taking into account then m quantum number multiplicity 1411 pure integer function ps_bound_niwfs(ps) 1412 type(ps_t), intent(in) :: ps 1413 1414 integer :: i, l 1415 1416 ps_bound_niwfs = 0 1417 do i = 1, ps%conf%p 1418 l = ps%conf%l(i) 1419 if (any(.not. ps%bound(i,:))) cycle 1420 ps_bound_niwfs = ps_bound_niwfs + (2*l+1) 1421 end do 1422 1423 end function ps_bound_niwfs 1424 1425 !--------------------------------------- 1426 1427 pure logical function ps_has_density(ps) result(has_density) 1428 type(ps_t), intent(in) :: ps 1429 1430 has_density = ps%has_density 1431 1432 end function ps_has_density 1433 1434 !--------------------------------------- 1435 1436 pure logical function ps_has_nlcc(ps) result(has_nlcc) 1437 type(ps_t), intent(in) :: ps 1438 1439 has_nlcc = ps%nlcc 1440 1441 end function ps_has_nlcc 1442 1443 !--------------------------------------- 1444 FLOAT function ps_density_volume(ps, namespace) result(volume) 1445 type(ps_t), intent(in) :: ps 1446 type(namespace_t), intent(in) :: namespace 1447 1448 integer :: ip, ispin 1449 FLOAT :: rr 1450 FLOAT, allocatable ::vol(:) 1451 type(spline_t) :: volspl 1452 1453 PUSH_SUB(ps_density_volume) 1454 1455 if (.not. ps_has_density(ps)) then 1456 message(1) = "The pseudopotential does not contain an atomic density" 1457 call messages_fatal(1, namespace=namespace) 1458 end if 1459 1460 SAFE_ALLOCATE(vol(1:ps%g%nrval)) 1461 1462 do ip = 1, ps%g%nrval 1463 rr = ps%g%rofi(ip) 1464 vol(ip) = CNST(0.0) 1465 do ispin = 1, ps%ispin 1466 vol(ip) = vol(ip) + spline_eval(ps%density(ispin), rr)*CNST(4.0)*M_PI*rr**5 1467 end do 1468 end do 1469 1470 call spline_init(volspl) 1471 call spline_fit(ps%g%nrval, ps%g%rofi, vol, volspl) 1472 volume = spline_integral(volspl) 1473 call spline_end(volspl) 1474 1475 SAFE_DEALLOCATE_A(vol) 1476 1477 POP_SUB(ps_density_volume) 1478 end function ps_density_volume 1479 1480#include "ps_pspio_inc.F90" 1481 1482end module ps_oct_m 1483 1484!! Local Variables: 1485!! mode: f90 1486!! coding: utf-8 1487!! End: 1488