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! 9 module dftu_specs 10! 11! Javier Junquera, January 2016, based on previous redftuproj 12! 13! Processes the information in an fdf file 14! to generate the projectors for the DFT+U simulations, 15! and populates the "projectors specifications" data structures. 16! 17! Here is a guide to the behavior of the main routine "read_dftu_specs": 18! 19! * Read the generation method of the projectors: 20! The DFT+U projectors are the localized functions used 21! to calculate the local populations used in a Hubbard-like term 22! that modifies the LDA Hamiltonian and energy. 23! It is important to recall that DFT+U projectors should be 24! quite localized functions. 25! Otherwise the calculated populations loose their atomic character 26! and physical meaning. Even more importantly, 27! the interaction range can increase so much that jeopardizes 28! the efficiency of the calculation. 29! 30! Two methods are currently implemented (accepted values are 1 and 2): 31! - If method_gen_dftu_proj = 1 32! Projectors are slightly-excited numerical atomic orbitals 33! similar to those used as an automatic basis set by SIESTA. 34! The radii of these orbitals are controlled using 35! the parameter DFTU.EnergyShift and/or the data 36! in block DFTU.proj (quite similar to the data block PAO.Basis used 37! to specify the basis set, see below). 38! - If method_gen_dftu_proj = 2 39! Projectors are exact solutions of the pseudoatomic 40! problem (and, in principle, are not strictly localized) which are 41! cutoff using a Fermi function $1/\{1+\exp[(r-r_c)\omega]\}$. 42! The values of $r_c$ and $\omega$ are controlled using 43! the parameter DFTU.CutoffNorm and/or the data 44! block DFTU.proj. 45! The default value is method_gen_dftu_proj = 2 46! 47! * Read the energy shift to generate the DFT+U projectors 48! Energy increased used to define the localization radious 49! of the DFTU projectors (similar to the parameter PAO.EnergyShift). 50! 51! * Allocate storage for the data structures 52! in particular the projector pointer that will be used later 53! in 54! * Determine any "global" basis specification parameters: 55! DFTU.proj - This is the most complex block, very flexible but in 56! need of spelling-out the specific restrictions. 57! It follows the same spirit as the PAO.Basis block. 58! Line by line, the specific info is: 59! 60! 1st: Species_label number_of_l_shells [basis_type] [ionic_charge] 61! 62! For each l_shell: 63! [n= n] l [E vcte rinn] 64! where 'n= n' is *mandatory* if the species has any semicore states, 65! E (soft confinement potential) section is optional. 66! we assume that only one projector per (n,l) quantum numbers is allowed 67! (i.e., nzeta = 1) 68! 69! U and J (in eV) 70! rc and omega (if method_gen_dftu_proj = 2), 71! where rc and omega are, respectively, 72! the equivalent of the Fermi energy and width 73! of the Fermi functions used to cut the 74! long wave functions 75! 76! or 77! 78! U and J (in eV) 79! rc (if method_gen_dftu_proj = 1), 80! where rc is the cutoff radius of the (short) wave 81! function used to generate the projector. 82! 83! The cutoff radii in Bohrs. These lines are mandatory. 84! 85! A line containing contraction (scale) factors is optional, but if it 86! appears, the values *must be* real numbers. 87! -------------------------------------------------------------------- 88! 89! After processing DFTU.proj block, whatever PAO information 90! which is not already available is determined in routine 'autobasis', 91! using the following defaults: 92! 93! rc(1:nzeta) is set to 0.0 94! lambda(1:nzeta) is set to 1.0 (this is a change from old practice) 95! 96! ---------------------------------- 97! 98! - If method_gen_dftu_proj = 1 99! The Schrodinger equation for the isolated atom is solved, with 100! the cutoff radius specified by the user. 101! The Schrodinger equation is solved using the 102! electrostatic potential generated by the "scaled" valence charge density. 103! 104! If the cutoff radius is set to zero, then it is determined by 105! the DFTU.EnergyShift parameter. If this is the case, the Schrodinger 106! equation has to be solved a first time, and this is done with the 107! electrostatic potential generated by the valence charge 108! density, readed from the pseudo-file. 109! 110! - If method_gen_dftu_proj = 2 111! The Schrodinger equation for the isolated atom is solved, with a 112! very large, arbitrary, cutoff radius. 113! Here, it is fixed by the rmax parameter to 60.0 Bohrs. 114! 115! The potential energy included in the Schrodinger equation is computed 116! from the rescaled charge density, as it was done in the generation of 117! the basis set. In other words, if we angularly integrate the 118! rescaled charge density between 0 and infinity it amounts 119! to (zval + charge), where zval is the nominal charge of the ion 120! (Zval (Si) = 4.0, Zval(Ba) = 10, and so on), 121! and charge is the charge included in the PAO.Basis block. 122! 123! Since the cutoff radius is very large, no soft-confinement is considered. 124! 125! Once the eigenfunctions are found for each angular momentum shell, 126! we search the point where the radial part drops below a threshold, 127! determined by the parameter min_func_val, fixed to 1.0d-6 128! The wave functions for larger distances will not be considered. 129! 130! Then, a Fermi-Dirac distribution is defined as 131! 1/[1+exp(r-rc)/width], 132! where width is defined in the DFTU.proj block and stored in the variable 133! dftu%width 134! and rc is explicitly given in the DFTU.proj block or computed 135! using the DFTU.CutoffNorm label, that defines the norm of the 136! original orbital contained in a sphere of the radius given by 137! this parameter. 138! 139! Finally, the long wave function is multiplied by the Fermi function 140! 141! To determine the cutoff of the DFTU+U projector, 142! we select the point where the previous prodcut is smaller than 143! a small tolerance, set up to 1.d-4. 144! 145! ======================================================================= 146! 147 use precision 148 149 use sys, only : die ! Termination routine 150 use basis_specs, only : label2species ! Function that reads the 151 ! label of a species and 152 ! converts it to the 153 ! corresponding index 154 ! according to the 155 ! Chemical_Species_Label block 156 use basis_types, only : basis_def_t ! Derived type where all the 157 ! information relative to the 158 ! definition of the basis set 159 ! is defined 160 use basis_types, only : shell_t ! Derived type of PAO shells 161 use basis_types, only : dftushell_t ! Derived type where all the 162 ! information relative to the 163 ! atomic orbitals where the U 164 ! correction will be applied 165 ! is defined 166 use basis_types, only : basis_parameters ! Derived type where all the 167 ! information about the 168 ! - basis set 169 ! - Kleinman-Bylander proj. 170 ! - DFT+U proj. 171 ! ... 172 ! for all the species 173 ! are defined 174 use basis_types, only: initialize ! Subroutine to initialize 175 ! the values of some derived 176 ! types 177 use basis_types, only: print_dftushell ! Subroutine to print 178 ! the values of the projectors 179 ! for DFT+U calculations 180 use basis_types, only : nsp ! Number of different 181 ! chemical species 182 use basis_types, only : charge ! Ionic charge to generate the 183 ! the basis set 184 use atmparams, only : nrmax ! Maximum number of points 185 ! in the logarithmic grid 186 use atmparams, only : lmaxd ! Maximum angular momentum 187 ! for both orbitals and 188 ! projectors. 189 use atmparams, only : NTBMAX ! Maximum number of points 190 ! in the tables defining 191 ! orbitals, projectors and 192 ! local neutral-atom pseudo 193 use pseudopotential, only: pseudopotential_t ! Derived type where all 194 ! the information about 195 ! the pseudopotential 196 ! is stored 197 use atom, only : schro_eq ! Subroutine to solve the 198 ! radial part of the 199 ! Schrodinger equation 200 use atom, only : rc_vs_e ! Subroutine to determine 201 ! the cutoff radius from the 202 ! energy shift 203 use atom, only : build_vsoft ! Subroutine to construct 204 ! the soft-confinement potent. 205 use atom_options,only: write_ion_plot_files ! Subroutine to plot the 206 ! basis functions and other 207 ! atomic functions 208 use atm_types, only : species_info ! Derived type with all the info 209 ! about the radial functions 210 ! (PAOs, KB projectors, 211 ! DFT+U proj, 212 ! VNA potentials, etc) 213 ! for a given atomic specie 214 use atm_types, only : species ! Actual array where the 215 ! previous information is 216 ! stored 217 use atm_types, only : nspecies ! Total number of different 218 ! atomic species 219 use units, only : pi, eV ! Conversions 220 use alloc, only : re_alloc ! Allocation routines 221 use radial ! Derived type for the radial 222 ! functions 223 use interpolation, only: spline ! set spline interpolation 224 use interpolation, only: polint ! polynomial interpolation 225 226 227 implicit none 228 229 integer :: method_gen_dftu_proj ! Method used to generate the 230 ! DFT+U projectors 231 ! Default value: exact solution 232 ! of the pseudoatomic problem 233 ! cutted with a Fermi function 234 real(dp) :: energy_shift_dftu ! Energy increase used to define 235 ! the localization radious of the DFT+U 236 ! projectors (similar to the parameter 237 ! PAO.EnergyShift) 238 ! Default value: 0.05 Ry 239 real(dp) :: dnrm_rc ! Parameter used to define the cutoff 240 ! radius that enters the 241 ! Fermi distribution to cut the 242 ! DFT+U projectors. 243 ! Only used if method_gen_dftu_proj = 2 244 ! It is the norm of the original 245 ! pseudoatomic orbital contained in 246 ! a sphere of radius dnrm_rc. 247 ! Default value: 0.90 248 real(dp) :: width_fermi_dftu ! Parameter used to define the width of 249 ! Fermi distribution to cut the 250 ! DFT+U projectors. 251 ! Only used if method_gen_dftu_proj = 2 252 ! Default value: 0.05 253 real(dp) :: dtol_dftupop ! Parameter that defines the 254 ! convergence criterium for the 255 ! DFT+U local population 256 real(dp) :: dDmax_threshold ! Parameter that defines the 257 ! criterium required to start or update 258 ! the calculation of the populations of 259 ! the DFT+U projections 260 logical :: dftu_init ! Flag that determines whether the 261 ! local populations are calculated 262 ! on the first iteration 263 logical :: dftu_shift ! Flag that determines whether the 264 ! parameter is interpreted 265 ! as a local potential shift 266 real(dp), pointer :: projector(:,:,:) ! Radial part of the DFT+U projector 267 integer, save, public, pointer :: nprojsdftu(:) 268 ! Total number of DFT+U projectors 269 ! (including the different angular 270 ! dependencies): i.e. a radial projector 271 ! with d-character counts as 5 different 272 ! DFT+U projectors 273 real(dp), parameter :: rmax = 60.0_dp 274 ! Arbitrary long localization radius 275 real(dp), parameter :: min_func_val = 1.e-6_dp 276 ! Minimum value of the 277 ! wave function times r, below which 278 ! the long wave function is no longer 279 ! considered 280 logical, save :: switch_dftu = .false. 281 ! Switch that determines whether 282 ! and DFT+U simulation is required or not 283 284 ! Routines 285 public :: read_dftu_specs 286 public :: dftu_proj_gen 287 public :: populate_species_info_dftu 288 289 ! Variables 290 public :: switch_dftu 291 public :: dftu_shift 292 public :: dftu_init 293 public :: dtol_dftupop 294 public :: dDmax_threshold 295 296 private 297 298 CONTAINS 299 300! subroutine read_dftu_specs : Subroutine that reads all the 301! info in the fdf file related with the 302! DFT+U projectors and 303! allocate some space for 304! the projector pointer 305! subroutine dftu_proj_gen : Subroutine that solves the 306! Schrodinger equation for the 307! isolated atom and generates the 308! DFT+U projectors 309! subroutine fermicutoff : Subroutine that computes the Fermi 310! function used to cut the long 311! atomic wave functions and produce 312! the DFT+U projectors. 313! only used if 314! method_gen_dftu_proj = 2 315! subroutine populate_species_info_dftu: Subroutine that populates the 316! data structures related with the DFT+U 317! projectors in the species 318! derived types. 319! Called from the atm_transfer subroutine 320 321!--- 322 subroutine read_dftu_specs() 323 324 use fdf 325 use m_cite, only: add_citation 326 use parallel, only: IONode 327 328 type(basis_def_t), pointer :: basp 329 type(dftushell_t), pointer :: dftu 330 type(dftushell_t), pointer :: lsdftu 331 type(shell_t), pointer :: shell 332 333 334 type(block_fdf) :: bfdf 335 type(parsed_line), pointer :: pline 336 337 integer :: isp ! Dummy parameter to account for the 338 ! species label 339 integer :: ish, jsh, i ! Dummy parameters to account for the 340 ! loop on shells 341 integer :: indexp ! Dummy parameters to account for the 342 ! reading of lines in DFTU.proj block 343 integer :: l ! Angular quantum number 344 integer :: maxnumberproj ! Maximum number of projectors 345 ! considered in a given species 346 347 logical :: bool 348 349! Default Soft-confinement parameters set by the user 350 logical, save :: lsoft 351 real(dp), save :: softRc, softPt 352 353!------------------------------------------------------------------------ 354! Read the generation method for the DFT+U projectors 355 method_gen_dftu_proj = 356 & fdf_get('LDAU.ProjectorGenerationMethod',2) 357 method_gen_dftu_proj = 358 & fdf_get('DFTU.ProjectorGenerationMethod',method_gen_dftu_proj) 359 360! Read the energy-shift to define the cut-off radius of the DFT+U projectors 361 energy_shift_dftu = 362 & fdf_get('LDAU.EnergyShift',0.05_dp,'Ry') 363 energy_shift_dftu = 364 & fdf_get('DFTU.EnergyShift',energy_shift_dftu,'Ry') 365 366! Read the parameter used to define the cutoff radius used in the Fermi 367! distribution 368! (only used if method_gen = 2) 369 dnrm_rc = fdf_get('LDAU.CutoffNorm',0.90_dp) 370 dnrm_rc = fdf_get('DFTU.CutoffNorm',dnrm_rc) 371 372! Read information about defaults for soft confinement 373 lsoft = fdf_get('PAO.SoftDefault', .false. ) 374 softRc = fdf_get('PAO.SoftInnerRadius', 0.9d0 ) 375 softPt = fdf_get('PAO.SoftPotential', 40.0d0 ) 376! Sanity checks on values 377 softRc = max(softRc,0.00d0) 378 softRc = min(softRc,0.99d0) 379 softPt = abs(softPt) 380 381! Read the parameter that defines the criterium required to start or update 382! the calculation of the populations of 383! the DFT+U projections 384 dDmax_threshold = fdf_get('LDAU.ThresholdTol', 1.e-2_dp) 385 dDmax_threshold = fdf_get('DFTU.ThresholdTol', dDmax_threshold) 386 387! Read the parameter that defines the convergence criterium for the 388! DFT+U local population 389 dtol_dftupop = fdf_get('LDAU.PopTol',1.e-3_dp) 390 dtol_dftupop = fdf_get('DFTU.PopTol',dtol_dftupop) 391 392! Read the flag that determines whether the U parameter is interpreted 393! as a local potential shift 394 dftu_shift = fdf_get('LDAU.PotentialShift', .false.) 395 dftu_shift = fdf_get('DFTU.PotentialShift', dftu_shift) 396 397 ! Read the flag that determines whether the local populations are 398 ! calculated on the first iteration 399 dftu_init = fdf_get('LDAU.FirstIteration', dftu_shift ) 400 dftu_init = fdf_get('DFTU.FirstIteration', dftu_init ) 401 ! When the local potential shift is applied 402 ! the initial iteration is forced to calculate 403 ! the DFT+U terms 404 if ( dftu_shift ) dftu_init = .true. 405 406! Allocate and initialize the array with the number of projectors per 407! atomic specie 408 nullify( nprojsdftu ) 409 call re_alloc( nprojsdftu, 1, nsp, 'nprojsdftu', 410 . 'read_dftu_specs' ) 411 nprojsdftu(:) = 0 412 413! Read the DFTU.proj block 414 if ( .not. fdf_block('LDAU.proj', bfdf) ) then 415 if (.not. fdf_block('DFTU.proj', bfdf)) RETURN 416 end if 417 418 ! Add citation 419 if ( IONode ) then 420 call add_citation("10.1103/PhysRevB.57.1505") 421 end if 422 423 do while(fdf_bline(bfdf, pline)) !! over species 424 if (.not. fdf_bmatch(pline,'ni')) 425 . call die('Wrong format in DFTU.proj') 426 isp = label2species(fdf_bnames(pline,1)) 427 if (isp .eq. 0) then 428 write(6,'(a,1x,a)') 429 . 'WRONG species symbol in DFTU.proj:', 430 . trim(fdf_bnames(pline,1)) 431 call die() 432 endif 433 434 basp => basis_parameters(isp) 435 436! Read on how many orbitals of a given atomic species 437! we are going to apply the U correction 438 basp%ndftushells = fdf_bintegers(pline,1) 439 440! Allocate space in the derived type basis_parameters 441! to host the information on the atomic orbitals where 442! the U will be applied 443 allocate(basp%dftushell(basp%ndftushells)) 444 445! Loop on all the different orbitals where the U will be applied 446 shells: do ish = 1, basp%ndftushells 447 dftu => basp%dftushell(ish) 448 call initialize(dftu) 449 450 if (.not. fdf_bline(bfdf, pline)) 451 . call die('Not enough information on the AO in DFTU.proj') 452 453! Read the principal and angular quantum numbers of the atomic orbital 454! where the U will be applied 455! Also we check what is the maximum value of the angular quantum 456! number between all of them that are read 457 458! In the DFTU.proj block, the information about the projectors 459! can be given as: 460! n=3 2 # n, l 461! i.e. with a string "n=" and then two integers... 462 if (fdf_bmatch(pline,'nii')) then 463 dftu%n = fdf_bintegers(pline,1) 464 dftu%l = fdf_bintegers(pline,2) 465 basp%lmxdftupj = max(basp%lmxdftupj,dftu%l) 466 467! or deleting the string "n=" 468! 3 2 # n, l 469! i.e. only with two integers 470 elseif (fdf_bmatch(pline,'ii')) then 471 dftu%n = fdf_bintegers(pline,1) 472 dftu%l = fdf_bintegers(pline,2) 473 basp%lmxdftupj = max(basp%lmxdftupj,dftu%l) 474 475! or only with one integer. In this case, this is the 476! angular quantum number. 477! If the semicore states is included in the valence, then 478! this is not valid and the principal quantum number has to be given 479! explicitly 480! 2 # l 481 elseif (fdf_bmatch(pline,'i')) then 482 if (basp%semic) 483 . call die('Please specify n if there are semicore states') 484 dftu%l = fdf_bintegers(pline,1) 485 dftu%n = basp%ground_state%n(dftu%l) 486 basp%lmxdftupj = max(basp%lmxdftupj,dftu%l) 487 else 488 call die('Bad format of (n), l line in DFTU.proj') 489 endif 490 491! Check for consistency in the sequence of principal and 492! angular quantum numbers 493 do jsh = 1, ish-1 494 if( dftu%l .eq. basp%dftushell(jsh)%l .and. 495 . dftu%n .eq. basp%dftushell(jsh)%n ) 496 . call die( 497 . 'DFTU projs. with the same l need different values of n') 498 enddo 499 500 ! Check that the principal and angular quantum numbers 501 ! are already an PAO 502 bool = .false. 503 do i = 0 , basp%lmxo 504 do jsh = 1 , basp%lshell(i)%nn 505 shell => basp%lshell(i)%shell(jsh) 506 if ( shell%nzeta == 0 ) cycle 507 bool = bool .or. 508 & (shell%n == dftu%n .and. shell%l == dftu%l) 509 end do 510 end do 511 if ( .not. bool ) then 512 call die( 513 &'DFTU projs require quantum numbers to exist. Check n and L') 514 end if 515 516! Check whether soft-confinement will be used 517 if (fdf_bsearch(pline,'E',indexp)) then 518 if (fdf_bmatch(pline,'vv',after=indexp)) then 519 dftu%vcte = fdf_bvalues(pline,ind=1,after=indexp) 520 dftu%rinn = fdf_bvalues(pline,ind=2,after=indexp) 521 else 522 call die('Need vcte and rinn after E in DFTU.proj') 523 endif 524 elseif (lsoft) then 525 dftu%vcte = softPt 526 dftu%rinn = -softRc 527 else 528 dftu%vcte = 0.0_dp 529 dftu%rinn = 0.0_dp 530 endif 531 532! Read the U and J parameters for this atomic orbital 533 if (.not. fdf_bline(bfdf, pline)) 534 . call die('No information for the U and J parameters...') 535 536 if ( fdf_bnvalues(pline) .ne. 2) 537 . call die('Insert values for the U and J parameters') 538 539 if (fdf_bmatch(pline,'vv')) then 540 dftu%u = fdf_bvalues(pline,1)*eV 541 dftu%j = fdf_bvalues(pline,2)*eV 542 endif 543 544! Read the cutoff radii (rc) to generate the projectors 545! and the contraction functions (lambda) 546 if ( .not. fdf_bline(bfdf, pline) ) 547 . call die('No information for the rc for projectors...') 548 549 if( method_gen_dftu_proj .eq. 1 ) then 550 if ( fdf_bnvalues(pline) .ne. 1 ) 551 . call die('Insert one value for the rc') 552 dftu%rc = fdf_bvalues(pline,1) 553 554 elseif( method_gen_dftu_proj .eq. 2 ) then 555 if ( fdf_bnvalues(pline) .ne. 2 ) 556 . call die('Insert one value for the rc and width') 557 dftu%rc = fdf_bvalues(pline,1) 558 dftu%dnrm_rc = dnrm_rc 559 if ( fdf_bvalues(pline,2) .lt. 1.d-4 ) then 560! Default value of parameter used to define the width of 561! the Fermi distribution to produce the DFT+U projectors 562! (only used if method_gen = 2) 563 dftu%width = 0.05_dp 564 else 565 dftu%width = fdf_bvalues(pline,2) 566 endif 567 568 endif 569 570! Optional: read the value for the contraction factor (lambda) 571 if ( .not. fdf_bline(bfdf,pline) ) then 572 if (ish.ne.basp%ndftushells) 573 . call die('Not enough shells') 574! Dafault values for the scale factors 575 else 576 if (.not. fdf_bmatch(pline,'r')) then 577 ! New shell or species 578 ! Default values for the scale factors 579 if ( .not. fdf_bbackspace(bfdf) ) 580 . call die('read_dftu_specs: ERROR in DFTU.proj block') 581! cycle shells 582 else 583 if ( fdf_bnreals(pline) .ne. 1 ) 584 . call die('One optional value of lambda') 585 dftu%lambda = fdf_breals(pline,1) 586 endif 587 endif 588 589 enddo shells ! end of loop over shells for species isp 590 591 592! Count the total number of projectors 593 nprojsdftu(isp) = 0 594 do ish = 1, basp%ndftushells 595 dftu => basp%dftushell(ish) 596 l = dftu%l 597 nprojsdftu(isp) = nprojsdftu(isp) + (2*l + 1) 598 enddo 599 basp%ndftuprojs_lm = nprojsdftu(isp) 600!! For debugging 601! write(6,'(a,i5)')'read_dftu_specs: lmxkb = ', 602! . basp%lmxkb 603! write(6,'(a,i5)')'read_dftu_specs: lmxdftupj = ', 604! . basp%lmxdftupj 605! write(6,'(a,i5)')'read_dftu_specs: nprojsdftu(isp) = ', 606! . nprojsdftu(isp) 607!! End debugging 608 609 enddo ! end loop over species 610 611! Allocate and initialize the array where the radial part of the 612! projectors in the logarithmic grid will be stored 613 maxnumberproj = 0 614 do isp = 1, nsp 615 basp => basis_parameters(isp) 616 maxnumberproj = max( maxnumberproj, basp%ndftushells ) 617 enddo 618 nullify( projector ) 619 call re_alloc( projector, 620 . 1, nsp, 621 . 1, maxnumberproj, 622 . 1, nrmax, 623 . 'projector', 'dftu_proj_gen' ) 624 projector = 0.0_dp 625 626 627!! For debugging 628! call die('Testing read_dftu_specs') 629!! End debugging 630 631 632 end subroutine read_dftu_specs 633 634! --------------------------------------------------------------------- 635! --------------------------------------------------------------------- 636 637 subroutine dftu_proj_gen( isp ) 638! --------------------------------------------------------------------- 639! Generation of DFT+U projectors 640! DFT+U projectors are, basically, pseudo-atomic-orbitals 641! with artificially small radii. 642! Written by D. Sanchez-Portal, Aug. 2008 after module basis_gen 643! Rewritten by Javier Junquera to merge with the top of the trunk 644! (Siesta 4.0), Feb. 2016 645! --------------------------------------------------------------------- 646 use basis_specs, only : restricted_grid 647 use basis_specs, only : rmax_radial_grid 648 649 use siestaXC, only : atomXC 650 651 integer, intent(in) :: isp ! Species index 652 653! Internal variables 654 integer :: n ! Principal quantum number of the projector 655 integer :: l ! Angular quantum number of the projector 656 integer :: lpseudo ! Angular quantum number of the pseudopotential 657 integer :: iproj ! Counter for the loop on projectors 658 integer :: ir ! Counter for the loop on points in the log grid 659 integer :: ndown ! Counter for the loop on l for the pseudos 660 integer :: ndftupj ! Number of DFT+U projectors that will be computed 661 ! for a given specie (here we consider only 662 ! different radial parts) 663 integer :: nodd ! Check whether we have and odd number of points 664 ! in the logarithmic grid 665 integer :: nnodes ! Number of nodes in the radial part of the 666 ! eigenfunctions of the Schrodinger equation 667 integer :: nprin ! Principal quantum number within the pseudoatom 668 real(dp) :: U ! Value of the U parameter 669 real(dp) :: J ! Value of the J parameter 670 real(dp) :: r2 ! Square of the distance to the nuclei 671 real(dp) :: rco ! Cutoff radius 672 real(dp) :: rc ! Cutoff radius (auxiliary variable to fit in an 673 ! odd number of points in the log grid) 674 real(dp) :: phi ! Wave function times r at a given point in 675 ! the log grid 676 real(dp) :: lambda ! Contraction factor 677 real(dp) :: el ! Energy of the eigenvalue after adding the 678 ! energy shift 679 real(dp) :: dnorm ! Norm of the projector 680 real(dp) :: rinn ! Inner radius where the soft-confinement potent. 681 ! starts off 682 real(dp) :: vcte ! Prefactor of the soft-confinement potent. 683 real(dp) :: ionic_charge! Ionic charge to generate the basis set. 684! Variables used only in the call to atomxc 685 real(dp) :: ex ! Total exchange energy 686 real(dp) :: ec ! Total correlation energy 687 real(dp) :: dx ! IntegralOf( rho * (eps_x - v_x) ) 688 real(dp) :: dc ! IntegralOf( rho * (eps_c - v_c) ) 689 690 real(dp) :: eigen(0:lmaxd) ! Eigenvalues of the Schrodinger equation 691 real(dp) :: rphi(nrmax,0:lmaxd) ! Eigenvectors of the Schrodinger equation 692 real(dp) :: vsoft(nrmax) ! Soft-confinement potential 693 real(dp) :: fermi_func(nrmax) ! Fermi function used to cut the 694 ! long pseudowave functions and 695 ! produce the DFT+U projectors 696 697! 698! Derived types where some information on the different shells are stored 699! 700 701 type(basis_def_t), pointer :: basp ! Parameters that define the 702 ! basis set, KB projectors, 703 ! DFT+U projectors, pseudopot 704 ! etc for a given species 705 type(dftushell_t), pointer :: shell ! Information about 706 ! DFT+U projectors 707 type(pseudopotential_t), pointer :: vps ! Pseudopotential information 708 709! 710! Variables related with the radial logarithmic grid 711! 712 integer :: nr ! Number of points required to 713 ! store the pseudopotential and 714 ! the wave functions in the 715 ! logarithmic grid 716 ! (directly read from the 717 ! pseudopotential file) 718 integer :: nrval ! Actual number of points in the 719 ! logarithmic grid 720 integer :: nrc ! Number of points required to 721 ! store the pseudowave functions 722 ! in the logarithmic grid 723 ! after being strictly confined. 724 integer :: nrwf ! Actual number of points in the 725 ! logarithmic grid to solve the 726 ! Schrodinger equation of the isolated 727 ! atom when no cutoff radius is 728 ! specified 729 ! In these cases, an arbitrary long 730 ! localization radius of 60.0 Bohrs 731 ! is assumed 732 integer :: nrwf_new ! 733 real(dp) :: a ! Step parameter of log. grid 734 ! (directly read from the 735 ! pseudopotential file) 736 real(dp) :: b ! Scale parameter of log. grid 737 ! (directly read from the 738 ! pseudopotential file) 739 real(dp) :: rofi(nrmax) ! Radial points of the 740 ! logarithmic grid 741 ! rofi(r)=b*[exp(a*(i-1)) - 1] 742 ! (directly read from the 743 ! pseudopotential file) 744 real(dp) :: drdi(nrmax) ! Derivative of the radial 745 ! distance respect the mesh index 746 ! Computed after the radial mesh 747 ! is read 748 real(dp) :: s(nrmax) ! Metric array 749 ! Computed after the radial mesh 750 ! is read 751 real(dp) :: rpb, ea ! Local variables used in the 752 ! calculation of the log. grid 753 754! 755! Variable used to store the semilocal component of the pseudopotential 756! 757! 758 character*4 :: nicore ! Flag that determines whether 759 ! non-linear core corrections 760 ! are included 761 character*3 :: irel ! Flag that determines whether 762 ! the atomic calculation is 763 ! relativistic or not 764 real(dp) :: vpseudo(nrmax,0:lmaxd) ! Semilocal components of the 765 ! pseudopotentials 766 ! (directly read from the 767 ! pseudopotential file) 768 real(dp) :: zval ! Valence charge of the atom 769 ! (directly read from the 770 ! pseudopotential file) 771 ! This value is the nominal one 772 773 774! 775! Variable used to store the semilocal component of the pseudopotential 776! 777 real(dp) :: chgvps ! Valence charge of the pseudoion 778 ! for which the pseudo was 779 ! generated in the ATM code 780 ! (it might not correspond with 781 ! the nominal valence charge 782 ! of the atom if the pseudo 783 ! has been generated for an ionic 784 ! configuration, for instance 785 ! when semicore has been 786 ! explicitly included in the 787 ! valence). 788 ! For instance, for Ba with 789 ! the semicore in valence, 790 ! (atomic reference configuration 791 ! 5s2 5p6 5d0 4f0), 792 ! chgvps = 8 (two in the 5s 793 ! and six in the 5p) 794 ! zval = 10 (two in the 5s, 795 ! six in the 5p, 796 ! and two in the 6s. 797 ! These two last electrons were 798 ! not included in the 799 ! reference atomic configuration) 800 real(dp) :: rho(nrmax) ! Valence charge density 801 ! As read from the pseudo file, 802 ! it is angularly integrated 803 ! (i.e. multiplied by 4*pi*r^2). 804 real(dp) :: rho_PAO(nrmax) ! Valence charge density 805 ! it is angularly integrated 806 ! (i.e. multiplied by 4*pi*r^2). 807 real(dp) :: ve(nrmax) ! Electrostatic potential 808 ! generated by the valence charge 809 ! density, readed from the 810 ! pseudo file 811 real(dp) :: vePAO(nrmax) ! Electrostatic potential 812 ! generated by the "scaled" 813 ! valence charge density 814 real(dp) :: vePAOsoft(nrmax) ! vePAO + the soft-confinement pot. 815 real(dp) :: vxc(nrmax) ! Exchange and correlation potentil 816 real(dp) :: chcore(nrmax) ! Core charge density 817 ! As read from the pseudo file, 818 ! it is angularly integrated 819 ! (i.e. multiplied by 4*pi*r^2). 820 real(dp) :: auxrho(nrmax) ! Sum of the valence charge and 821 ! core charge (if NLCC included) 822 ! densities to compute the 823 ! atomic exchange and correl. 824 ! potential. 825 ! auxrho is NOT angularly integr. 826 ! (not multiplied by 4*pi*r^2) 827 integer :: irelt ! Flag that determines whether the 828 ! atomic calculation to 829 ! generate the pseudopotential 830 ! was relativistic (irelt = 1) 831 ! or no relativistic (irelt = 0) 832 real(dp), parameter :: eps = 1.0e-4_dp ! Epsilon value used to 833 ! determine the cutoff of 834 ! the DFT+U projector 835 ! if method = 2 836 837 838! Associate the pointer so it points to the variable where all the 839! parameters defining the basis sets of the given species are stored 840 basp => basis_parameters(isp) 841 842! Determine if something has to be done regarding the 843! generation of the DFT+U projectors. 844! If DFT+U is not required (number of DFT+U projectors equal to zero), 845! then do nothing and return. 846 847! Compute how many DFT+U projector we are going to compute 848! for this species 849 ndftupj = basp%ndftushells 850 851 852! Determine whether the calculation of DFT+U projectors is required or not 853! for this atomic species 854 if( .not. ndftupj > 0 ) return 855 856! This switch will be used afterwards in setup_hamiltonian 857! to determine whether the call to the Hubbard subroutine 858! is required or not 859 switch_dftu = .true. 860 861! Associate the pointer so it points to the variable where all the 862! parameters defining the basis sets of the given species are stored 863 vps => basp%pseudopotential 864 865! 866! Read all the required information from the pseudopotentials that 867! will be required to solve the Schrodinger equation for the isolated atoms 868! 869 nr = vps%nr 870 b = vps%b 871 a = vps%a 872 zval = vps%zval 873 nicore = vps%nicore 874 irel = vps%irel 875 ionic_charge = charge(isp) 876 877 nrval = nr + 1 878 if (rmax_radial_grid /= 0.0_dp) then 879 nrval = nint(log(rmax_radial_grid/b+1.0d0)/a)+1 880 write(6,"(a,f10.5,i5)") 881 . 'Maximum radius (at nrval) set to ', 882 . rmax_radial_grid, nrval 883 endif 884 885 if (restricted_grid) then 886 nodd = mod(nrval,2) 887 nrval = nrval -1 + nodd ! Will be less than or equal to vp%nrval 888 endif 889 890 if ( nrval .gt. nrmax ) then 891 write(6,'(a,i4)') 892 . 'dftu_proj_gen: ERROR: Nrmax must be increased to at least', 893 . nrval 894 call die 895 endif 896 897! Read the radial logarithmic mesh 898 rofi(1:nrval) = vps%r(1:nrval) 899 900! Calculate drdi and s 901! drdi is the derivative of the radial distance respect to the mesh index 902! i.e. rofi(ir)= b*[ exp( a*(i-1) ) - 1 ] and therefore 903! drdi=dr/di =a*b*exp(a*(i-1))= a*[rofi(ir)+b] 904 905 rpb = b 906 ea = dexp(a) 907 do ir = 1, nrval 908 drdi(ir) = a * rpb 909 s(ir) = dsqrt( a * rpb ) 910 rpb = rpb * ea 911 enddo 912 913!! For debugging 914!! Differences with respect Daniel's grid implementation 915!! can appear in the number of points nrval. 916!! In this latest version, the option restricted_grid is activated 917!! by default. 918!! This was not yet implemented in the version where Daniel started 919!! Even more, the value of s printed by Daniel corresponds 920!! with the redefinition given below for the integration of the 921!! Schrodinger equation (s = drdi^2). 922!! The one defined here is required for the solution of the Poisson equation 923! do ir = 1, nrval 924! write(6,'(i5,3f20.12)') ir, rofi(ir), drdi(ir), s(ir) 925! enddo 926! call die() 927!! End debugging 928 929! 930! Read the ionic pseudopotentials (Only 'down' used) 931! These are required to solve the Schrodinger equation for the isolated 932! atoms. 933! Here we read all the semilocal components of the pseudopotential 934! independently of whether the DFT+U projector for a particular 935! angular momentum is required or not. 936! 937 do 20 ndown = 1, basp%lmxdftupj+1 938 939 lpseudo = vps%ldown(ndown) 940 941 if ( lpseudo .ne. ndown-1 ) then 942 write(6,'(a)') 943 . 'dftu_proj_gen: Unexpected angular momentum for pseudopotential' 944 write(6,'(a)') 945 . 'dftu_proj_gen: Pseudopot. should be ordered by increasing l' 946 endif 947 948 vpseudo(1:nrval,lpseudo) = vps%vdown(ndown,1:nrval) 949 950 do ir = 2, nrval 951 vpseudo(ir,lpseudo) = vpseudo(ir,lpseudo)/rofi(ir) 952 enddo 953 vpseudo(1,lpseudo) = vpseudo(2,lpseudo) ! AG 954 955 20 enddo 956 957!! For debugging 958!! Up to this point, these are the same pseudos as read in 959!! Daniel's version of DFT+U 960!! The only difference might be at the number of points in 961!! the log grid 962! do lpseudo = 0, basp%lmxdftupj 963! write(6,'(/a,i5)') 964! . ' dftu_proj_gen: Reading pseudopotential for l = ', 965! . lpseudo 966! 967! do ir = 1, nrval 968! write(6,'(a,i5,2f20.12)') 969! . ' ir, rofi, vpseudo = ', ir, rofi(ir), vpseudo(ir,lpseudo) 970! enddo 971! enddo 972!! End debugging 973 974! Read the valence charge density from the pseudo file 975! and scale it if the ionic charge of the reference configuration 976! is not the same as the nominal valence charge of the atom 977 chgvps = vps%gen_zval 978 do ir = 1, nrval 979 rho(ir) = chgvps * vps%chval(ir)/zval 980 enddo 981 982! Find the Hartree potential created by a radial electron density 983! using the Numerov's method to integrate the radial Poisson equation. 984! The input charge density at this point has to be angularly integrated. 985 call vhrtre( rho, ve, rofi, drdi, s, nrval, a ) 986 987 988! Set 'charge': 989! 1. If 'charge' is not set in the fdf file (in the PAO.Basis block, 990! an charge can be included to generate the pseudopotential) 991! then set it to zval-chgvps. 992! 2. If 'charge' is equal to zval-chgvps, set it to that. 993! 994 if( ( abs(ionic_charge) .eq. huge(1.0_dp) ) .or. 995 . ( abs( ionic_charge-(zval-chgvps) ) .lt. 1.0d-3) ) then 996 ionic_charge = zval - chgvps 997 endif 998 999! For DFT+U projector calculations 1000! We use the "scaled" charge density of an ion of total charge "charge" 1001! As seen above, this ion could be the one involved in ps generation, 1002! or another specified at the time of basis generation. 1003! Example: Ba: ps ionic charge: +2 1004! basis gen specified charge: +0.7 1005! if we integrate the charge density in rho_PAO, the integral 1006! would amount to (zval - ionic_charge) = 10.0 - 0.7 = 9.3 1007 1008 do ir = 2,nrval 1009 rho_PAO(ir) = (zval-ionic_charge) * rho(ir) / chgvps 1010 enddo 1011 rho_PAO(1) = rho_PAO(2) 1012 1013!! For debugging: check the normalization condition of the rescaled charge 1014!! density 1015! dnorm = 0.0_dp 1016! do ir = 2,nrval 1017! dnorm = dnorm + rho_PAO(ir) * drdi(ir) 1018! enddo 1019! write(6,'(a,4f12.5)')'Total charge in rho_PAO = ', 1020! . dnorm, zval, ionic_charge, (zval-ionic_charge) 1021!! End debugging 1022 1023 call vhrtre( rho_PAO, vePAO, rofi, drdi, s, nrval, a ) 1024 1025!! For debugging 1026! write(6,'(a,3f12.5)')' zval, chgvps, ionic_charge = ', 1027! . zval, chgvps, ionic_charge 1028! do ir = 1, nrval 1029! write(6,'(a,i5,4f20.12)') 1030! . ' ir, rofi, rho, ve, vePAO = ', 1031! . ir, rofi(ir), rho(ir), ve(ir), vePAO(ir) 1032! enddo 1033!! call die() 1034!! End debugging 1035 1036! Read the core charge density from the pseudo file 1037 chcore(1:nrval) = vps%chcore(1:nrval) 1038 1039! Compute the exchange and correlation potential in the atom 1040! Note that atomxc expects true rho(r), not 4 * pi * r^2 * rho(r) 1041! We use auxrho for that 1042! 1043 do ir = 2, nrval 1044 r2 = rofi(ir)**2 1045 r2 = 4.0_dp * pi * r2 1046 dc = rho(ir) / r2 1047 if( nicore .ne. 'nc ') dc = dc + chcore(ir) / r2 1048 auxrho(ir) = dc 1049 enddo 1050 r2 = rofi(2) / (rofi(3)-rofi(2)) 1051 auxrho(1) = auxrho(2) - ( auxrho(3) - auxrho(2) ) * r2 1052 1053! Determine whether the atomic calculation to generate the pseudopotential 1054! is relativistic or not 1055 if (irel.eq.'rel') irelt=1 1056 if (irel.ne.'rel') irelt=0 1057 1058! Compute the exchange and correlation potential 1059 call atomxc( irelt, nrval, nrmax, rofi, 1060 & 1, auxrho, ex, ec, dx, dc, vxc ) 1061 1062!! For debugging 1063! write(6,'(a,i5)') 'irelt = ', irelt 1064! write(6,'(a,i5)') 'nrval = ', nrval 1065! write(6,'(a,i5)') 'nrmax = ', nrmax 1066! do ir = 1, nrval 1067! write(6,'(a,i5,3f20.12)') 1068! . ' ir, rofi, auxrho, vxc = ', 1069! . ir, rofi(ir), auxrho(ir), vxc(ir) 1070! enddo 1071! call die() 1072!! End debugging 1073 1074! Add the exchange and correlation potential to the Hartree potential 1075 ve(1:nrval) = ve(1:nrval) + vxc(1:nrval) 1076 1077!! For debugging 1078! write(6,'(a,f20.12)')' chg = ', chgvps 1079! write(6,'(a,f20.12)')' a = ', a 1080! write(6,'(a,f20.12)')' b = ', b 1081! do ir = 1, nrval 1082! write(6,'(a,i5,3f20.12)') 1083! . ' ir, rofi, vxc+vhr, vpseudo = ', 1084! . ir, rofi(ir), ve(ir), vpseudo(ir,0) 1085! enddo 1086! call die() 1087!! End debugging 1088 1089 do ir = 2,nrval 1090 r2 = rofi(ir)**2 1091 r2 = 4.0d0*pi*r2 1092 dc = rho_PAO(ir)/r2 1093 if ( nicore .ne. 'nc ' ) dc = dc + chcore(ir)/r2 1094 auxrho(ir) = dc 1095 enddo 1096 r2 = rofi(2)/(rofi(3)-rofi(2)) 1097 auxrho(1) = auxrho(2) -(auxrho(3)-auxrho(2))*r2 1098 1099 call atomxc( irelt, nrval, nrmax, rofi, 1100 & 1, auxrho, ex, ec, dx, dc, vxc ) 1101 1102 vePAO(1:nrval) = vePAO(1:nrval) + vxc(1:nrval) 1103 1104!! For debugging 1105! do ir = 1, nrval 1106! write(6,'(a,i5,4f20.12)') 1107! . ' ir, rofi, rho, ve, vePAO = ', 1108! . ir, rofi(ir), rho(ir), ve(ir), vePAO(ir) 1109! enddo 1110!! call die() 1111!! End debugging 1112 1113! 1114! Redefine the array s for the Schrodinger equation integration 1115! 1116 s(2:nrval) = drdi(2:nrval) * drdi(2:nrval) 1117 s(1) = s(2) 1118 1119! Loop over all the projectors that will be generated 1120 loop_projectors: do iproj = 1, ndftupj 1121 shell => basp%dftushell(iproj) 1122 n = shell%n 1123 l = shell%l 1124 U = shell%u 1125 J = shell%j 1126 rco = shell%rc 1127 rinn = shell%rinn 1128 vcte = shell%vcte 1129 1130! If the compression factor is negative or zero, 1131! the orbitals are left untouched 1132 if( shell%lambda .le. 0.0d0 ) shell%lambda=1.0d0 1133 lambda = shell%lambda 1134 1135! Check whether the cutoff radius for the DFT+U projector 1136! is explicitly determined in the input file or automatically 1137! controlled by 1138! - the EnergyShift parameter (method_gen_dftu_proj = 1) 1139! - the cutoff of the Fermi distribution (method_gen_dftu_proj = 2) 1140! For the first generation method, and if we rely on the automatic 1141! determination, we have to compute the rc from the EnergyShift. 1142! This is done from a cut-and-paste from the corresponding lines 1143! in the generation of the PAOs for the basis sets. 1144! For the second generation method, 1145! the rc is determined by the Fermi function, 1146! and it is done in fermicutoff distribution 1147 if ( rco .lt. 1.0d-5 ) then 1148 1149! Cutoff controled by the energy shift parameter: 1150 if( method_gen_dftu_proj .eq. 1) then 1151! Some required variables to solve the Schrodinger 1152! equation are defined below 1153 1154! Determine the number of nodes in the radial part of 1155! the eigenfunction 1156! THIS HAS TO BE UPDATED WITH THE SUBROUTINES 1157! OF THE NEW PSEUDOS: 1158! FROM THE KNOWLEDGE OF n AND l, IT SHOULD BE POSSIBLE 1159! TO DETERMINE THE NUMBER OF NODES 1160 nnodes = 1 1161! Determine the principal quantum number within the pseudoatom 1162 nprin = l + 1 1163 1164 nrwf = nrval 1165 if (restricted_grid) nrwf = nrwf + 1 - mod(nrwf,2) 1166 1167!! For debugging 1168! write(6,'(/a,i2)') 1169! . 'DFTUprojs with principal quantum number n = ', n 1170! write(6,'(a,i2)') 1171! . 'DFTUprojs with angular momentum l = ', l 1172! write(6,'(a,f12.5)') 1173! . 'DFTUprojs with U = ',U 1174! write(6,'(a,f12.5)') 1175! . 'DFTUprojs with J = ',J 1176! write(6,'(a,f12.5)') 1177! . 'DFTUprojs with lambda = ',shell%lambda 1178! write(6,'(a,f12.5)') 1179! . 'DFTUprojs with rc = ',shell%rc 1180! write(6,'(a,i5)') 1181! . 'DFTUprojs with nnodes = ',nnodes 1182! write(6,'(a,i5)') 1183! . 'DFTUprojs with nprin = ',nprin 1184! write(6,'(a,i5)') 1185! . 'DFTUprojs with nrval = ',nrval 1186! write(6,'(a,i5)') 1187! . 'DFTUprojs with nrwf = ',nrwf 1188! write(6,'(a,i5)') 1189! . 'DFTUprojs with nrmax = ',nrmax 1190! write(6,'(a,f12.5)') 1191! . 'DFTUprojs with zval = ',zval 1192! write(6,'(a,f12.5)') 1193! . 'DFTUprojs with a = ',a 1194! write(6,'(a,f12.5)') 1195! . 'DFTUprojs with b = ',b 1196!! End debugging 1197 1198! Initialize the eigenfunctions 1199 rphi(:,l) = 0.0_dp 1200! Solve the Schrodinger for the long cutoff 1201 call schro_eq( zval, rofi, vpseudo(1,l), ve, s, drdi, 1202 . nrwf, l, a, b, nnodes, nprin, 1203 . eigen(l), rphi(1,l) ) 1204 1205! 1206! Compute the cutoff radius of the DFT+U projectors 1207! as given by energy_shift_dftu 1208! 1209 if( eigen(l) .gt. 0.0_dp ) then 1210 write(6,'(/a,i2,a)') 1211 . 'dftu_proj_gen: ERROR Orbital with angular momentum L=', 1212 . l, ' not bound in the atom' 1213 write(6,'(a)') 1214 . 'dftu_proj_gen: an rc radius must be explicitly given' 1215 call die() 1216 endif 1217 1218 if( abs(energy_shift_dftu) .gt. 1.0d-5 ) then 1219 el = eigen(l) + energy_shift_dftu 1220 call rc_vs_e( a, b, rofi, vpseudo(1,l), ve, nrval, l, 1221 . el, nnodes, rco ) 1222 else 1223 rco = rofi(nrval-2) 1224 endif 1225 1226! Store the new variable for the cutoff radii 1227! automatically determined 1228 shell%rc = rco 1229 1230 write(6,'(/,a,/,a,f10.6,a)') 1231 . 'dftu_proj_gen: PAO cut-off radius determined from an', 1232 . 'dftu_proj_gen: energy shift =',energy_shift_dftu,' Ry' 1233 write(6,'(a,f10.6,a)') 1234 . 'dftu_proj_gen: rco =',rco,' Bohr' 1235 1236 endif 1237 1238 endif ! End if automatic determination of the rc 1239 1240! At this point, independently of the method, 1241! we should now the cutoff radius of the DFT+U projector. 1242! Now, we compute it 1243 rco = shell%rc 1244 1245! Store the radial point of the logarithmic grid where the 1246! DFT+U projector vanishes 1247 nrc = nint(log(rco/b+1.0_dp)/a)+1 1248 shell%nrc = nrc 1249 1250! Determine the number of nodes 1251 nnodes = 1 1252! Determine the principal quantum number within the pseudoatom 1253 nprin = l + 1 1254 1255 if( method_gen_dftu_proj .eq. 1) then 1256! Build the soft confinement potential 1257 vsoft = 0.0_dp 1258! Scale the orbitals with the contraction factor 1259 rc = rco / lambda 1260 call build_vsoft( isp, l, 1, rinn, vcte, 1261 . 0.0_dp, 0.0_dp, 0.0_dp, 1262 . a, b, rc, rofi, nrval, 1263 . vsoft, plot=write_ion_plot_files ) 1264!! For debugging 1265! write(6,'(/a,i5)') '# l = ' , l 1266! write(6,'(a,f12.5)')'# Eigenvalue = ', eigen(l) 1267! write(6,'(a)') '# Soft-confinement ' 1268! write(6,'(a,f12.5)')'# Inner radius = ' , rinn 1269! write(6,'(a,f12.5)')'# Prefactor = ', vcte 1270! write(6,'(a,f12.5)')'# Cutoff radius = ', rco 1271! do ir = 1, nrval 1272! write(6,'(2f20.12)')rofi(ir), vsoft(ir) 1273! enddo 1274!! End debugging 1275 1276 do ir = 1, nrval 1277 vePAOsoft(ir) = vePAO(ir) + vsoft(ir) 1278 enddo 1279 1280! 1281! If rc is negative, treat it as a fractional value 1282 1283 if (rco .lt. 0.0_dp) then 1284 call die("rc < 0 for first-zeta orbital") 1285 endif 1286 1287! Find the required number of points in the logarithmic grid 1288! to solve the Scrodingcer equation 1289 nrc = nint(log(rc/b+1.0_dp)/a)+1 1290 1291! Note that rco is redefined here, to make it fall on an odd-numbered 1292! grid point. 1293! 1294 if (restricted_grid) then 1295 nodd = mod(nrc,2) 1296 if( nodd .eq. 0 ) then 1297 nrc = nrc + 1 1298 endif 1299 endif 1300 1301 rc = b*(exp(a*(nrc-1))-1.0d0) 1302 1303! Solve the Schrodinger equation for the required cutoff 1304! and with the Hartree potential from the scaled charge density 1305 1306! Initialize the eigenfunctions 1307 rphi(:,l) = 0.0_dp 1308 call schro_eq( zval, rofi, vpseudo(1,l), vePAOsoft, s, drdi, 1309 . nrc, l, a, b, nnodes, nprin, 1310 . eigen(l), rphi(1,l) ) 1311 1312! Normalize the eigenfunctions 1313! and divide them by r^(l+1) 1314! In the previous subroutine, we compute r * phi, 1315! where phi is the radial part of the wave functions. 1316! In Siesta, we store in the tables phi/r^l. 1317! Therefore, we need to divide the previous solution by 1318! r^(l+1) 1319 1320 projector(isp,iproj,:)=rphi(:,l)/(rofi(:)**(l+1)) 1321 projector(isp,iproj,1)=projector(isp,iproj,2) 1322 1323 dnorm = 0.0_dp 1324 do ir = 2, nrc 1325 dnorm = dnorm + drdi(ir) * 1326 . (projector(isp,iproj,ir)*rofi(ir)**(l+1))**2 1327 enddo 1328 dnorm = sqrt(dnorm) 1329 projector(isp,iproj,:) = projector(isp,iproj,:)/dnorm 1330 projector(isp,iproj,1) = projector(isp,iproj,2) 1331 1332 shell%nrc = nrc 1333 shell%rc = rc 1334 1335 else if( method_gen_dftu_proj .eq. 2) then 1336! An arbitrary long localization radius for these orbitals 1337! is set up with the parameter rmax (=60.0 Bohr by default). 1338! This was suggested in the original implementation by Daniel 1339! and is kept here for backwards compatibility 1340 nrwf = nint(log(rmax/b+1.0d0)/a)+1 1341 nrwf = min(nrwf,nrval) 1342 if (restricted_grid) nrwf = nrwf + 1 - mod(nrwf,2) 1343 1344! For debugging 1345 write(6,'(/a,i2)') 1346 . 'DFTUprojs with principal quantum number n = ', n 1347 write(6,'(a,i2)') 1348 . 'DFTUprojs with angular momentum l = ', l 1349 write(6,'(a,f12.5)') 1350 . 'DFTUprojs with U = ',U 1351 write(6,'(a,f12.5)') 1352 . 'DFTUprojs with J = ',J 1353 write(6,'(a,f12.5)') 1354 . 'DFTUprojs with lambda = ',shell%lambda 1355 write(6,'(a,f12.5)') 1356 . 'DFTUprojs with rc = ',shell%rc 1357 write(6,'(a,i5)') 1358 . 'DFTUprojs with nnodes = ',nnodes 1359 write(6,'(a,i5)') 1360 . 'DFTUprojs with nprin = ',nprin 1361 write(6,'(a,i5)') 1362 . 'DFTUprojs with nrval = ',nrval 1363 write(6,'(a,i5)') 1364 . 'DFTUprojs with nrwf = ',nrwf 1365 write(6,'(a,i5)') 1366 . 'DFTUprojs with nrmax = ',nrmax 1367 write(6,'(a,f12.5)') 1368 . 'DFTUprojs with zval = ',zval 1369 write(6,'(a,f12.5)') 1370 . 'DFTUprojs with a = ',a 1371 write(6,'(a,f12.5)') 1372 . 'DFTUprojs with b = ',b 1373! End debugging 1374 1375! Initialize the eigenfunctions 1376 rphi(:,l) = 0.0_dp 1377 1378! Solve the Schrodinger for the long cutoff 1379 call schro_eq( zval, rofi, vpseudo(1,l), vePAO, s, drdi, 1380 . nrwf, l, a, b, nnodes, nprin, 1381 . eigen(l), rphi(1,l) ) 1382 1383! We consider only the solutions to the Schrodinger 1384! equation up to the point where its value is smaller than 1385! a given tolerance, setup by the min_func_val parameter 1386 nrwf_new = nrwf 1387 do ir = nrwf, 2, -1 1388 if( abs(rphi(ir,l) ) .gt. min_func_val ) then 1389 nrwf_new = ir + 1 1390 write(6,'(a,f20.12,a)') 1391 . 'dftu_proj_gen: updating the rc to', 1392 . rofi(nrwf_new), ' Bohr' 1393 exit 1394 endif 1395 enddo 1396 1397! Divide the eigenfunctions by r^(l+1) and normalize them. 1398! In the previous subroutine, we compute r * phi, 1399! where phi is the radial part of the wave functions. 1400! In Siesta, we store in the tables phi/r^l. 1401! Therefore, we need to divide the previous solution by 1402! r^(l+1) 1403 do ir = 2, nrwf_new 1404 rphi(ir,l)=rphi(ir,l)/(rofi(ir)**(l+1)) 1405 enddo 1406 rphi(1,l)=rphi(2,l) 1407! Nullify the rest of the solution 1408 rphi(nrwf_new+1:nrmax,l) = 0.0_dp 1409 1410 dnorm = 0.0_dp 1411 do ir = 1, nrwf_new 1412 phi = rphi(ir,l) 1413! To compute the norm, we need to integrate 1414! r^2 \times wave_function^2. 1415! Since we have stored wave_function/r^l, we need to 1416! multiply it by r^(l+2) 1417 dnorm = dnorm + drdi(ir) * (phi * rofi(ir)**(l+1))**2 1418 enddo 1419 dnorm = dsqrt(dnorm) 1420 do ir = 2, nrwf_new 1421 rphi(ir,l) = rphi(ir,l)/dnorm 1422 enddo 1423 1424! Now, define the Fermi distribution that will be used 1425! to cut the long eigenfunction 1426! The width of the Fermi distribution is defined by 1427! the shell%width parameter 1428! while the equivalent of the Fermi energy is determined by 1429! the shell%dnrm parameter 1430 call fermicutoff( nrmax, nrwf_new, rofi, drdi, 1431 . rphi(:,l), shell, fermi_func ) 1432 1433!! For debugging 1434! write(6,'(/a,i5)') '# l = ' , l 1435! write(6,'(a,f12.5)')'# Eigenvalue = ', eigen(l) 1436! write(6,'(a,f12.5)')'# rc = ', shell%rc 1437! write(6,'(a,f12.5)')'# width = ', shell%width 1438! write(6,'(a,f12.5)')'# Norm = ', dnorm 1439! write(6,'(a)') '# Eigenfunction ' 1440! do ir = 1, nrwf_new 1441! write(6,'(3f20.12)')rofi(ir), rphi(ir,l), fermi_func(ir) 1442! enddo 1443!! End debugging 1444 1445! Here we multiply the long wave function times the 1446! Fermi-Dirac distribution to cut it. 1447 projector(isp,iproj,:) = 0.0_dp 1448 do ir = 1, nrwf_new 1449 projector(isp,iproj,ir) = fermi_func(ir) * rphi(ir,l) 1450 enddo 1451 1452! Normalize the projector 1453 dnorm = 0.0_dp 1454 do ir = 1, nrwf_new 1455! Here we have stored projector/r^l, where projector is the 1456! radial part of the DFT+U projector. 1457! To compute the norm in spherical coordinates, 1458! we have to integrate \int r^{2} projector^{2} dr, 1459! and this implies that we have to mutiply rphi**2 times r^(l+1)**2 1460 dnorm = dnorm + drdi(ir) * 1461 . (projector(isp,iproj,ir)*rofi(ir)**(l+1))**2 1462 enddo 1463 dnorm = dsqrt(dnorm) 1464 projector(isp,iproj,:) = projector(isp,iproj,:) / dnorm 1465 projector(isp,iproj,1) = projector(isp,iproj,2) 1466 1467! 1468! Set up the cutoff for the DFT+U projector 1469! 1470 do ir = nrwf_new, 1, -1 1471 if(dabs(projector(isp,iproj,ir)) .gt. eps) exit 1472 enddo 1473 shell%nrc = ir+1 1474 shell%rc = rofi(ir+1) 1475 1476! Normalize after the cut 1477 dnorm = 0.0_dp 1478 do ir = 1, shell%nrc 1479! The projector that has been stored in the array projector 1480! is written in the same format as the atomic orbitals in the 1481! inners of Siesta, i. e., in the format of R/r^l, 1482! where R is the radial part of the projector. 1483! To check if it is normalized, 1484! we have to integrate \int r^{2} R^{2} dr, 1485! and this implies just to take projector**2 1486! and multiply by r^(2l+2) = r^(2*(l+1)) 1487 dnorm = dnorm + drdi(ir)*(projector(isp,iproj,ir)**2)* 1488 . rofi(ir)**(2*(l+1)) 1489 enddo 1490 dnorm = dsqrt(dnorm) 1491 projector(isp,iproj,:) = projector(isp,iproj,:) / dnorm 1492 projector(isp,iproj,1) = projector(isp,iproj,2) 1493 1494 endif 1495 1496!! For debugging 1497! dnorm = 0.0_dp 1498! do ir = 1, shell%nrc 1499!! The projector that has been stored in the array projector 1500!! is written in the same format as the atomic orbitals in the 1501!! inners of Siesta, i. e., in the format of R/r^l, 1502!! where R is the radial part of the projector. 1503!! To check if it is normalized, 1504!! we have to integrate \int r^{2} R^{2} dr, 1505!! and this implies just to take projector**2 1506!! and multiply by r^(2l+2) = r^(2*(l+1)) 1507! dnorm = dnorm + drdi(ir)*(projector(isp,iproj,ir)**2)* 1508! . rofi(ir)**(2*(l+1)) 1509! enddo 1510! 1511! write(6,'(/a,i5)') '# l = ', l 1512! write(6,'(a,f12.5)')'# Eigenvalue = ', eigen(l) 1513! write(6,'(a,f12.5)')'# Cutoff = ', shell%rc 1514! write(6,'(a,i7)') '# Number of points = ', shell%nrc 1515! write(6,'(a,f12.5)')'# Width = ', shell%width 1516! write(6,'(a)') '# Projector ' 1517! write(6,'(a,f12.5)')'# Norm = ', dnorm 1518! do ir = 1, shell%nrc 1519! write(6,'(2f20.12)')rofi(ir), projector(isp,iproj,ir) 1520! enddo 1521!! End debugging 1522 1523 enddo loop_projectors ! End the loop on projectors 1524 1525!! For debugging 1526! call die("Testing dftu_proj_gen") 1527!! End debugging 1528 1529 1530 end subroutine dftu_proj_gen 1531! --------------------------------------------------------------------- 1532 1533 subroutine fermicutoff( nrmax, nrval, rofi, drdi, rphi, 1534 . dftushell, fermi_func ) 1535! 1536! This subroutine defines the fermi function used to cut the long 1537! atomic wave functions and produce the DFT+U projectors 1538! Only used if method_gen_dftu_proj = 2 1539! 1540 1541 integer, intent(in) :: nrmax ! Maximum number of points 1542 ! of the log grid 1543 ! (required to define 1544 ! the dimensions of 1545 ! some arrays). 1546 integer, intent(in) :: nrval ! Number of points of the 1547 ! logarithmic grid where 1548 ! the long wave function 1549 ! is computed 1550 real(dp), intent(in) :: rofi(nrmax) ! Points of the log grid 1551 real(dp), intent(in) :: drdi(nrmax) ! Distance between 1552 ! consecutive points of 1553 ! the log grid 1554 real(dp), intent(in) :: rphi(nrmax) ! Long wave functions 1555 ! (eigenfunctions) 1556 ! of the Schrodinger 1557 ! equation 1558 type(dftushell_t),intent(inout) :: dftushell 1559 real(dp), intent(out) :: fermi_func(nrmax) ! Fermi function 1560 1561! Internal vars 1562 integer :: ir ! Counter for the loops on real space 1563 ! grids 1564 integer :: l ! Angular momentum of the shell 1565 real(dp) :: rc ! "Fermi energy" of the Fermi function 1566 real(dp) :: width ! Width of the Fermi function 1567 real(dp) :: a ! Auxiliary function to compute the 1568 ! Fermi function 1569 real(dp) :: dnorm ! Norm of the original pseudoatomic 1570 ! wave function 1571 real(dp), parameter :: gexp = 60.0_dp 1572 real(dp), parameter :: eps = 1.0e-4_dp ! A small value (epsilon) 1573 ! for comparison 1574 1575! Initialize the angular momentum quantum number. 1576 l = dftushell%l 1577 1578! If no cutoff distance is explicitly given in the input file 1579! (DFTU.proj block) then compute the cutoff distance for the Fermi function 1580! For this, we have to check at which radial distance 1581! the norm of the original pseudo atomic orbital equals 1582! the value introduced in DFTU.CutoffNorm 1583 1584 if ( dftushell%rc .lt. eps ) then 1585 1586 dnorm = 0.0_dp 1587 do ir = 1, nrmax 1588! In rphi we have stored phi/r^l, where phi is the radial part of the 1589! wave function. 1590! To compute the norm in spherical coordinates, 1591! we have to integrate \int r^{2} R^{2} dr, 1592! and this implies that we have to mutiply rphi**2 times r^(l+1)**2 1593 dnorm = dnorm + drdi(ir) * (rphi(ir)*rofi(ir)**(l+1))**2 1594 if( dnorm .gt. dnrm_rc ) exit 1595 enddo 1596 dftushell%rc = rofi(ir) 1597 endif 1598 1599! Initialize Fermi function 1600 fermi_func = 0.0_dp 1601 1602! Determine the parameters of the Fermi distribution 1603 rc = dftushell%rc 1604 width = dftushell%width 1605 1606 do ir = 1, nrval 1607 a = ( rofi(ir) - rc ) / width 1608 if( a .lt. -gexp ) then 1609 fermi_func(ir) = 1.0_dp 1610 else if( a .gt. gexp ) then 1611 fermi_func(ir) = 0.0_dp 1612 else 1613 fermi_func(ir) = 1.0_dp / ( 1.0_dp+dexp(a) ) 1614 endif 1615 enddo 1616 1617!! For debugging 1618! write(6,'(a,f12.5)')'# Fermi function computed with rc = ', rc 1619! write(6,'(a,f12.5)')'# and width = ', width 1620! do ir = 1, nrval 1621! write(6,'(2f20.12)') 1622! . rofi(ir), fermi_func(ir) 1623! enddo 1624!! call die() 1625!! End debugging 1626! 1627 end subroutine fermicutoff 1628 1629! ---------------------------------------------------------------------- 1630 subroutine populate_species_info_dftu 1631! 1632! In this subroutine, we populate the variables in the species_info 1633! derived type related with the DFT+U projectors. 1634! It is called from atm_transfer. 1635! 1636 use alloc, only : de_alloc 1637 type(species_info), pointer :: spp 1638 type(basis_def_t), pointer :: basp 1639 type(dftushell_t), pointer :: dftushell 1640 type(rad_func), pointer :: pp 1641 type(pseudopotential_t), pointer :: vps 1642 1643! Internal variables 1644 integer :: is ! Counter for the loop on atomic species 1645 integer :: iproj ! Counter for the loop on projectors 1646 integer :: ir ! Counter for the loop on real space points 1647 integer :: l ! Quantum angular momentum of a given DFT+U proj. 1648 integer :: im ! Counter for the loop magnetic quantum number 1649 integer :: imcount ! 1650 integer :: nr ! Point in the log. grid closest to the linear grid 1651 integer :: nn ! Total number of points in the log grid considered 1652 ! for the interpolation 1653 integer :: nmin ! nr - npoint (see below for the meaning of npoint) 1654 integer :: nmax ! nr + npoint (see below for the meaning of npoint) 1655 real(dp) :: rc ! Cutoff radius of the different DFT+U proj. 1656 integer :: nrc ! Point in the log. grid where the DFT+U proj. vanish 1657 real(dp) :: delta ! Interval between consecutive points in the grid 1658 ! where the DFT+U projectors are stored 1659 real(dp) :: rpoint ! Coordinate of the real space points 1660 real(dp) :: projint ! Interpolated value of the DFT+U projector at rpoint 1661 real(dp) :: dy ! Function derivative at point rpoint 1662 real(dp) :: a ! Parameters of the logarithmic grid 1663 real(dp) :: b ! Parameters of the logarithmic grid 1664 real(dp) :: yp1 ! First derivative at the first point of the grid 1665 real(dp) :: ypn ! First derivative at the last point of the grid 1666 real(dp) :: rofi(nrmax) ! Radial points of the 1667 ! logarithmic grid 1668 ! rofi(r)=b*[exp(a*(i-1)) - 1] 1669 ! (directly read from the 1670 ! pseudopotential file) 1671 real(dp) :: projinputint(nrmax) ! Radial part of the projector that 1672 ! enters the interpolation routines 1673 1674 integer, parameter :: npoint = 4 ! Number of points used by polint 1675 ! for the interpolation 1676 1677! Loop on different atomic species 1678 loop_species: do is = 1, nspecies 1679 spp => species(is) 1680 basp => basis_parameters(is) 1681 vps => basp%pseudopotential 1682 1683! Read the parameters for the logarithmic grid 1684 a = vps%a 1685 b = vps%b 1686 nr = vps%nr 1687! Read the radial logarithmic mesh 1688 rofi(1:nr) = vps%r(1:nr) 1689 1690! Store the total number of DFT+U projectors 1691! counting the "m copies" 1692! (including the (2l + 1) factor for each l). 1693 spp%nprojsdftu = nprojsdftu(is) 1694 1695! Number of DFT+U projectors 1696! not counting the "m copies" 1697 spp%n_pjdftunl = basp%ndftushells 1698 1699! Store the maximum angular momentum of the DFT+U projectors 1700! for each atomic specie 1701 spp%lmax_dftu_projs = basp%lmxdftupj 1702 1703! Loop on all the projectors for a given specie 1704! This loop is done only on the different radial shapes, 1705! without considering the (2l + 1) possible angular dependencies 1706 imcount = 0 1707 loop_projectors: do iproj = 1, spp%n_pjdftunl 1708 dftushell => basp%dftushell(iproj) 1709 1710 spp%pjdftunl_n(iproj) = 1 1711 spp%pjdftunl_l(iproj) = dftushell%l 1712 spp%pjdftunl_U(iproj) = dftushell%U 1713 spp%pjdftunl_J(iproj) = dftushell%J 1714 1715 l = spp%pjdftunl_l(iproj) 1716 1717 do im = -l, l 1718 imcount = imcount + 1 1719 spp%pjdftu_n(imcount) = dftushell%n 1720 spp%pjdftu_l(imcount) = dftushell%l 1721 spp%pjdftu_m(imcount) = im 1722 spp%pjdftu_index(imcount) = iproj 1723 enddo 1724 enddo loop_projectors ! End loop on projectors for a given specie 1725 if( imcount .ne. spp%nprojsdftu ) call die('DFT+U indexing...') 1726 1727! Allocate the derived types pjdftu, of radial kind, 1728! where the radial components of the DFT+U projectors will be stored 1729! There will be as many radial functions of this kind 1730! as different DFT+U projectors, without including the m copies. 1731 allocate ( spp%pjdftu(spp%n_pjdftunl) ) 1732 1733 do iproj = 1, spp%n_pjdftunl 1734 dftushell => basp%dftushell(iproj) 1735 pp => spp%pjdftu(iproj) 1736 call rad_alloc(pp,NTBMAX) 1737 rc = dftushell%rc 1738 nrc = dftushell%nrc 1739 delta = rc/(dble(ntbmax-1)+1.0d-20) 1740 pp%cutoff = rc 1741 pp%delta = delta 1742 1743 projinputint(:) = projector(is,iproj,:) 1744 1745! Interpolate the projectors from the logarithmic grid to the 1746! linear grid 1747 do ir = 1, ntbmax-1 1748 rpoint = delta * (ir-1) 1749 nr = nint(log(rpoint/b+1.0d0)/a)+1 1750 nmin = max( 1, nr-npoint ) 1751 nmax = min( nrc, nr+npoint ) 1752 nn = nmax - nmin + 1 1753 call polint( rofi(nmin), projinputint(nmin), 1754 . nn, rpoint, projint, dy ) 1755 pp%f(ir) = projint 1756 enddo 1757 1758! Compute the second derivative of the projectors 1759 call rad_setup_d2(pp,yp1=0.0_dp,ypn=huge(1.0_dp)) 1760 1761 enddo 1762 1763 enddo loop_species ! End loop on atomic species 1764 1765!! For debugging 1766! do is = 1, nspecies 1767! write(6,'(/a,i5)') 1768! . '# populate_species_info_dftu: specie number = ', 1769! . is 1770! 1771! basp => basis_parameters(is) 1772! spp => species(is) 1773! 1774! write(6,'(a,i5)') 1775! . '#populate_species_info_dftu: specie, spp%lmax_dftu_projs = ', 1776! . spp%lmax_dftu_projs 1777! write(6,'(a,i5)') 1778! . '#populate_species_info_dftu: specie, spp%n_pjdftunl = ', 1779! . spp%n_pjdftunl 1780! write(6,'(a)') 1781! . '#populate_species_info_dftu: Loop over different projectors' 1782! write(6,'(a)') 1783! . '#populate_species_info_dftu: not considering m copies ' 1784! write(6,'(a)') 1785! . '#populate_species_info_dftu: iproj, pjdftu_n, pjdftunl_l' 1786! 1787! do iproj = 1, spp%n_pjdftunl 1788! write(6,'(a,3i5)') 1789! . '#populate_species_info_dftu:', 1790! . iproj, spp%pjdftunl_n(iproj), spp%pjdftunl_l(iproj) 1791! pp => spp%pjdftu(iproj) 1792! write(6,'(a,f20.12)') 1793! . '#populate_species_info_dftu: cutoff = ', pp%cutoff 1794! write(6,'(a,f20.12)') 1795! . '#populate_species_info_dftu: delta = ', pp%delta 1796! do ir = 1, ntbmax-1 1797! rpoint = pp%delta * (ir-1) 1798! write(6,'(3f20.12)') rpoint, pp%f(ir), pp%d2(ir) 1799! enddo 1800! write(6,*) 1801! enddo 1802! 1803! write(6,'(a,i5)') 1804! . '#populate_species_info_dftu: specie, spp%nprojsdftu = ', 1805! . spp%nprojsdftu 1806! write(6,'(a)') 1807! . '#populate_species_info_dftu: Loop over different projectors' 1808! write(6,'(a)') 1809! . '#populate_species_info_dftu: considering m copies ' 1810! write(6,'(a)') 1811! . '#populate_species_info_dftu: iproj, pjdftu_n, l , m, index' 1812! do iproj = 1, spp%nprojsdftu 1813! dftushell => basp%dftushell(spp%pjdftu_index(iproj)) 1814! write(6,'(5i5,f12.5)') 1815! . iproj, spp%pjdftu_n(iproj), spp%pjdftu_l(iproj), 1816! . spp%pjdftu_m(iproj), spp%pjdftu_index(iproj),dftushell%rc 1817! enddo 1818! enddo 1819! call die('End testing populate_species_info_dftu') 1820!! End debugging 1821 1822 call de_alloc( projector, 'projector', 'dftu_proj_gen') 1823 call de_alloc( nprojsdftu, 'nprojsdftu', 'read_dftu_specs') 1824 1825 end subroutine populate_species_info_dftu 1826 1827 end module dftu_specs 1828