1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 module basis_specs 9! 10! Alberto Garcia, August 2000, based on 'classic' redbasis. 11! 12! Processes the basis information in an fdf file and populates 13! the "basis specifications" data structures. 14! 15! Here is a guide to the behavior of the main routine "read_basis_specs": 16! 17! * Find out the number of species 18! * Allocate storage for the data structures 19! * Determine any "global" basis specification parameters: 20! - basis_size (sz, dz, dzp, etc) 21! - basis_type (split, nodes, etc) ("USER" is no longer valid) 22! * Process the mandatory "Chemical-species-label" block to read the 23! species labels (not the same as chemical symbols) and atomic numbers. 24! (note that negative atomic numbers have a special meaning) 25! * Assign default values to 'basis_size', 'basis_type'. 26! * Assign default values to the 'mass' based on the atomic number. 27! * Set up the information about the ground state of the atom. 28! * Read information about the valence charge density in the 29! pseudopotential file and determine whether there are semicore 30! states. (Note that in this case the user is explicitly asked 31! (see below) to input the 'n' quantum numbers in the PAO.Basis block.) 32! * Read the optional fdf blocks: 33! AtomicMass (routine remass) 34! PAO.BasisSize - Overrides the default 'basis_size' on a per-species basis. 35! PAO.Basis - This is the most complex block, very flexible but in 36! need of spelling-out the specific restrictions. Line by 37! line, these are: 38! 39! 1st: Species_label number_of_l_shells [basis_type] [ionic_charge] 40! 41! For each l_shell: 42! [n= n] l nzeta [P [ nzeta_pol]] [E vcte rinn] [Q qcoe [qyuk [qwid]]] [F cutoff] 43! where 'n= n' is *mandatory* if the species has any semicore states, 44! and the 'P' (polarization), E (soft confinement potential), 45! 'Q' (charge confinement), and 'F' (filteret) sections are optional and can appear 46! in any order after nzeta. 47! 48! rc_1 rc_2 .... rc_nzeta 49! 50! are the cutoff radii in bohrs. This line is mandatory 51! If the number of rc's for a given shell is less than the number of 52! 'zetas', the program will assign the last rc value to the remaining 53! zetas, rather than stopping with an error. This is particularly 54! useful for Bessel suites of orbitals. 55! A line containing contraction (scale) factors is optional, but if it 56! appears, the values *must be* real numbers, and there must be at 57! least nzeta of them. 58! -------------------------------------------------------------------- 59! 60! After processing PAO.Basis (if it exists), whatever PAO information 61! which is not already available is determined in routine 'autobasis', 62! using the following defaults: 63! 64! (There is a first check to make sure that species with semicore 65! states, or with Bessel floating orbitals, have been included in 66! the PAO.Basis block) 67! 68! The maximum angular momentum for the basis (lmxo) (excluding any 69! polarization orbitals) is that of the ground state, as returned by 70! routine 'lmxofz' from Z. 71! 72! Each l-shell contains just one shell of PAOS, with n set 73! to the appropriate ground state value. 74! 75! Nzeta is determined by the first two characters of the 'basis_size' 76! string: 1 for 'sz' and 2 for 'dz'. 77! 78! There are no 'per l-shell' polarization orbitals, except if the 79! third character of 'basis_size' is 'p' (as in 'dzp'), in which 80! case polarization orbitals are defined so they have the minimum 81! angular momentum l such that there are no occupied orbitals 82! with the same l in the valence shell of the ground-state 83! atomic configuration. They polarize the corresponding l-1 shell. 84! 85! The Soft-Confinement parameters 'rinn' and 'vcte' are set to 0.0 86! The Charge-Confinement parameters 'qcoe', 'qyuk' and 'qwid' 87! are set to 0.0, 0.0 and 0.01 88! 89! rc(1:nzeta) is set to 0.0 90! lambda(1:nzeta) is set to 1.0 (this is a change from old practice) 91! 92! ---------------------------------- 93! 94! Next come the blocks associated to the KB projectors: 95! 96! Block PS.lmax (routine relmxkb) assigns a per-species maximum l for 97! the KB projectors. Internally, it sets the 'lmxb_requested' 98! component of the data structure, which by default is -1. 99! 100! Block PS.KBprojectors (optional) sets the number of KB projectors 101! per angular momentum. The same routine which reads this block 102! (readkb) also sets any appropriate defaults: 103! 104! 105! If the species does not appear in the PS.KBprojectors block, then 106! Setting of lmxkb: 107! If lmxkb is set in PS.lmax, 108! set it to that 109! else 110! Use PAO information: (routine set_default_lmxkb) 111! (Set it to a meaningless value for floating orbitals) 112! Set it to lmxo+1, or to lpol+1, where lpol is the angular 113! momentum of the highest-l polarization orbital. 114! endif 115! 116! There is a hard limit for lmxkb: if the pseudopotential file 117! contains semilocal pseudopotentials up to lmax_pseudo, then 118! lmxkb <= lmax_pseudo. 119! 120! The number of KB projectors per l is set to the number of 121! n-shells in the corresponding PAO shell with the same l. For l 122! greater than lmxo, it is set to 1. The reference energies are 123! in all cases set to huge(1.d0) to mark the default. 124! 125! Archaeological note: The first implementation of the basis-set generation 126! module had only non-polarization orbitals. Polarization orbitals were added 127! later as "second-class" companions. This shows in details like "lmxo" (the 128! maximum l of the basis set) not taking into account polarization orbitals. 129! Polarization orbitals were tagged at the end, without maintaining l-shell 130! ordering. 131! Even later, support for "semicore" orbitals was added. A new ('nsm') index was 132! used to distinguish the different orbitals in a l-shell. Polarization orbitals 133! were not brought into this classification. 134! The code is strained when semicore and polarization orbitals coexist for the same l, 135! as in Ti, whose electronic structure is []3s2 3p6 3d2 4s2 (4p0*) with 4p as the 136! polarization orbital. 137! Here 'nsemic' (the number of semicore states) is kept at zero for l=1 138! (the polarization orbital is not counted), with the side-effect that only one KB 139! projector is generated for l=1. 140! Special checks have been implemented to cover these cases. 141! 142! Future work should probably remove the separate treatment of polarization orbitals. 143! (Note that if *all* orbitals are specified in a PAO.Basis block, in effect turning 144! perturbative polarization orbitals into 'normal' orbitals, this problem is not present.) 145! 146! ======================================================================= 147! 148 use precision 149 use basis_types, only: basis_def_t, shell_t, lshell_t, kbshell_t 150 use basis_types, only: nsp, basis_parameters, ground_state_t 151 use basis_types, only: destroy, copy_shell, initialize 152 use pseudopotential, only: pseudo_read, pseudo_reparametrize 153 use pseudopotential, only: pseudo_init_constant 154 use periodic_table, only: qvlofz, lmxofz, cnfig, atmass 155 use chemical 156 use sys 157 use fdf 158 159 Implicit None 160 161 type(basis_def_t), pointer :: basp => null() 162 type(shell_t), pointer :: s => null() 163 type(lshell_t), pointer :: ls => null() 164 type(kbshell_t), pointer :: k => null() 165 166 character(len=1), parameter :: 167 $ sym(0:4) = (/ 's','p','d','f','g' /) 168 169! Default Soft-confinement parameters set by the user 170! 171 logical, save :: lsoft 172 real(dp), save :: softRc, softPt 173 174! Default norm-percentage for the automatic definition of 175! multiple-zeta orbitals with the 'SPLIT' option 176 177 real(dp), parameter :: splnorm_default=0.15_dp 178 real(dp), parameter :: splnormH_default=-1.0_dp 179 real(dp), save :: global_splnorm, global_splnorm_H 180 integer isp ! just an index dummy variable for the whole module 181 182! 183 logical, save, public :: restricted_grid 184 logical, parameter :: restricted_grid_default = .true. 185 186 real(dp), save, public :: rmax_radial_grid 187 188 public :: read_basis_specs 189 public :: label2species 190 191 private 192 193 CONTAINS 194 195!--- 196 subroutine read_basis_specs() 197 198 character(len=15), parameter :: basis_size_default='standard' 199 character(len=10), parameter :: basistype_default='split' 200 201 character(len=15) :: basis_size 202 character(len=10) :: basistype_generic 203 204 type(block_fdf) :: bfdf 205 type(parsed_line), pointer :: pline 206 207 type(ground_state_t), pointer :: gs 208 209 integer nns, noccs, i, ns_read, l 210 logical synthetic_atoms, found, reparametrize_pseudos 211 real(dp) :: new_a, new_b 212 213!------------------------------------------------------------------------ 214 reparametrize_pseudos = 215 $ fdf_boolean('ReparametrizePseudos',.false.) 216! 217! r(i) = b * ( exp(a*(i-1)) - 1 ) 218! spacing near zero: delta = a*b 219! If the desired spacing at r=R is delta*(1+beta), 220! then: 221! a = beta*delta/R 222! b = R/beta 223! N at R = 1 + R*ln(1+beta)/(beta*delta) 224! 225 if (reparametrize_pseudos) then 226 ! These default values will provide a grid spacing 227 ! of 1.0e-5 near r=0, and 0.01 near r=10 (bohr units) 228 ! with N at Rmax=100 on the order of 10000 points. 229 new_a = fdf_double("NewAParameter",0.001_dp) 230 new_b = fdf_double("NewBParameter",0.01_dp) 231 endif 232 233! 234! Whether or not to restrict integration grids to an odd number 235! of points (and hence to shift rcs accordingly) 236! 237 restricted_grid = 238 $ fdf_boolean("Restricted.Radial.Grid", 239 $ restricted_grid_default) 240! 241! If non-zero, the value will be the maximum value of the 242! radial coordinate 243! 244 if (reparametrize_pseudos) then 245 rmax_radial_grid = fdf_double('Rmax.Radial.Grid',50.0_dp) 246 else 247 rmax_radial_grid = fdf_double('Rmax.Radial.Grid',0.0_dp) 248 endif 249 250! 251 basis_size = fdf_string('PAO.BasisSize',basis_size_default) 252 call size_name(basis_size) 253 basistype_generic = fdf_string('PAO.BasisType',basistype_default) 254 call type_name(basistype_generic) 255 256C Read information about defaults for soft confinement 257 258 lsoft = fdf_boolean('PAO.SoftDefault',.false.) 259 softRc = fdf_double('PAO.SoftInnerRadius',0.9d0) 260 softPt = fdf_double('PAO.SoftPotential',40.0d0) 261 262C Sanity checks on values 263 264 softRc = max(softRc,0.00d0) 265 softRc = min(softRc,0.99d0) 266 softPt = abs(softPt) 267! 268! Read defaults for split_norm parameter 269 270 global_splnorm = fdf_double('PAO.SplitNorm',splnorm_default) 271 global_splnorm_H = fdf_double('PAO.SplitNormH',splnormH_default) 272 if (global_splnorm_H < 0.0_dp) global_splnorm_H = global_splnorm 273 274!------------------------------------------------------------------ 275! 276! Use standard routine in chemical module to process the 277! chemical species 278! 279 nsp = size(basis_parameters) 280 281 synthetic_atoms = .false. 282 283 do isp=1,nsp 284 basp=>basis_parameters(isp) 285 286 basp%label = species_label(isp) 287 basp%z = atomic_number(isp) 288 basp%floating = is_floating(isp) 289 basp%bessel = is_bessel(isp) 290 basp%synthetic = is_synthetic(isp) 291 292 basp%basis_size = basis_size 293 basp%basis_type = basistype_generic 294 if (basp%floating) then 295 basp%mass = 1.d40 ! big but not too big, as it is used 296 ! later in computations 297 else if (basp%synthetic) then 298 basp%mass = -1.0_dp ! Signal -- Set later 299 else 300 basp%mass = atmass(abs(int(basp%z))) 301 endif 302 if (basp%bessel) then 303 ! Initialize a constant pseudo 304 call pseudo_init_constant(basp%pseudopotential) 305 else if (basp%synthetic) then 306 synthetic_atoms = .true. 307 ! Will set gs later 308 call pseudo_read(basp%label,basp%pseudopotential) 309 else 310 call ground_state(abs(int(basp%z)),basp%ground_state) 311 call pseudo_read(basp%label,basp%pseudopotential) 312 endif 313 if (reparametrize_pseudos.and. .not. basp%bessel) 314 . call pseudo_reparametrize(p=basp%pseudopotential, 315 . a=new_a, b=new_b,label=basp%label) 316 enddo 317 318 if (synthetic_atoms) then 319 320 found = fdf_block('SyntheticAtoms',bfdf) 321 if (.not. found ) 322 . call die("Block SyntheticAtoms does not exist.") 323 ns_read = 0 324 do while(fdf_bline(bfdf, pline)) 325 326 ns_read = ns_read + 1 327 if (.not. fdf_bmatch(pline,'i')) 328 . call die("Wrong format in SyntheticAtoms") 329 isp = fdf_bintegers(pline,1) 330 if (isp .gt. nsp .or. isp .lt. 1) 331 . call die("Wrong specnum in SyntheticAtoms") 332 basp => basis_parameters(isp) 333 gs => basp%ground_state 334 if (.not. fdf_bline(bfdf, pline)) call die("No n info") 335 nns = fdf_bnintegers(pline) 336 if (nns .lt. 4) 337 . call die("Please give all valence n's " // 338 . "in SyntheticAtoms block") 339 gs%n = 0 340 do i = 1, nns 341 gs%n(i-1) = fdf_bintegers(pline,i) 342 enddo 343 if (.not. fdf_bline(bfdf, pline)) 344 . call die("No occupation info") 345 noccs = fdf_bnvalues(pline) 346 if (noccs .lt. nns) call die("Need more occupations") 347 gs%occupation(:) = 0.0_dp 348 do i = 1, noccs 349 gs%occupation(i-1) = fdf_bvalues(pline,i) 350 enddo 351 ! Determine the last occupied state in the atom 352 do i = nns, 1, -1 353 if (gs%occupation(i-1) /=0) then 354 gs%lmax_valence = i-1 355 exit 356 endif 357 enddo 358 gs%occupied(0:3) = (gs%occupation .gt. 0.0_dp) 359 gs%occupied(4) = .false. 360 gs%z_valence = sum(gs%occupation(0:noccs-1)) 361 write(6,'(a,i2)',advance='no') 362 . 'Ground state valence configuration (synthetic): ' 363 do l=0,3 364 if (gs%occupied(l)) 365 . write(6,'(2x,i1,a1,f8.5)',advance='no') 366 . gs%n(l),sym(l), gs%occupation(l) 367 enddo 368 write(6,'(a)') '' 369 370 enddo 371 write(6,"(a,i2)") "Number of synthetic species: ", ns_read 372 373 endif 374! 375! Defer this here in case there are synthetic atoms 376 do isp=1,nsp 377 basp=>basis_parameters(isp) 378 if (basp%synthetic) then 379 gs => basp%ground_state 380 if (gs%z_valence .lt. 0.001) 381 . call die("Synthetic species not detailed") 382 endif 383 call semicore_check(isp) 384 enddo 385 386 call remass() 387 call resizes() 388 call repaobasis() 389 call autobasis() 390 call relmxkb() 391 call readkb() 392 393! do isp=1,nsp 394! call print_basis_def(basis_parameters(isp)) 395! enddo 396 397 end subroutine read_basis_specs 398!----------------------------------------------------------------------- 399 400 function label2species(str) result(index) 401 integer index 402 character(len=*), intent(in) :: str 403 404 integer i 405 406 index = 0 407 do i=1,nsp 408 if(leqi(basis_parameters(i)%label,str)) index = i 409 enddo 410 end function label2species 411!----------------------------------------------------------------------- 412 413 subroutine ground_state(z,gs) 414 integer, intent(in) :: z 415 type(ground_state_t), intent(out) :: gs 416! 417! Determines ground state valence configuration from Z 418! 419 integer l, latm 420 421 gs%z_valence = 0.d0 422 do l=0,3 423 gs%occupation(l)=0.0d0 424 enddo 425 426 call lmxofz(z,gs%lmax_valence,latm) 427 call qvlofz(z,gs%occupation(:)) 428 do l=0,gs%lmax_valence 429 gs%z_valence = gs%z_valence + gs%occupation(l) 430 enddo 431 call cnfig(z,gs%n(0:3)) 432 433 write(6,'(a,i2)',advance='no') 434 . 'Ground state valence configuration: ' 435 gs%occupied(4) = .false. !! always 436 do l=0,3 437 gs%occupied(l) = (gs%occupation(l).gt.0.1d0) 438 if (gs%occupied(l)) 439 . write(6,'(2x,i1,a1,i2.2)',advance='no') 440 . gs%n(l),sym(l),nint(gs%occupation(l)) 441 enddo 442 write(6,'(a)') '' 443 444 end subroutine ground_state 445 446!--------------------------------------------------------------- 447 subroutine readkb() 448 integer lpol, isp, ish, i, l 449 character(len=20) unitstr 450 451 type(block_fdf) :: bfdf 452 type(parsed_line), pointer :: pline 453 454 integer :: lmax_pseudo 455 456 lpol = 0 457 458 if (fdf_block('PS.KBprojectors',bfdf) ) then 459 460! First pass to find out about lmxkb and set any defaults. 461 462 do while(fdf_bline(bfdf, pline)) !! over species 463 if (.not. fdf_bmatch(pline,'ni')) 464 . call die("Wrong format in PS.KBprojectors") 465 isp = label2species(fdf_bnames(pline,1)) 466 if (isp .eq. 0) then 467 write(6,'(a,1x,a)') 468 . "WRONG species symbol in PS.KBprojectors:", 469 . trim(fdf_bnames(pline,1)) 470 call die() 471 endif 472 basp => basis_parameters(isp) 473 basp%nkbshells = fdf_bintegers(pline,1) 474 do ish= 1, basp%nkbshells 475 if (.not. fdf_bline(bfdf, pline)) call die("No l nkbl") 476 if (.not. fdf_bmatch(pline,'ii')) 477 . call die("Wrong format l nkbl") 478 l = fdf_bintegers(pline,1) 479 if (l .gt. basp%lmxkb) basp%lmxkb = l 480 if (.not. fdf_bline(bfdf, pline)) then 481 if (ish .ne. basp%nkbshells) 482 . call die("Not enough kb shells for this species...") 483 ! There is no line with ref energies 484 else if (fdf_bmatch(pline,'ni')) then 485 ! We are seeing the next species' section 486 if (ish .ne. basp%nkbshells) 487 . call die("Not enough shells for this species...") 488 if (.not. fdf_bbackspace(bfdf)) 489 . call die('readkb: ERROR in PS.KBprojectors block') 490 else if (fdf_bmatch(pline,'ii')) then 491 ! We are seeing the next shell's section 492 if (ish .gt. basp%nkbshells) 493 . call die("Too many kb shells for this species...") 494 if (.not. fdf_bbackspace(bfdf)) 495 . call die('readkb: ERROR in PS.KBprojectors block') 496 endif 497 enddo ! end of loop over shells for species isp 498 499 enddo 500 endif 501 502 do isp=1,nsp 503! 504! Fix defaults and allocate kbshells 505! 506 basp=>basis_parameters(isp) 507 if (basp%lmxkb .eq. -1) then ! not set in KBprojectors 508 if (basp%lmxkb_requested.eq.-1) then ! not set in PS.lmax 509 basp%lmxkb = set_default_lmxkb(isp) ! Use PAO info 510 else 511 basp%lmxkb = basp%lmxkb_requested 512 endif 513 allocate(basp%kbshell(0:basp%lmxkb)) 514 do l=0,basp%lmxkb 515 call initialize(basp%kbshell(l)) 516 k=>basp%kbshell(l) 517 k%l = l 518 if (l.gt.basp%lmxo) then 519 k%nkbl = 1 520 else 521 ! Set equal to the number of PAO shells with this l 522 k%nkbl = basp%lshell(l)%nn 523 ! Should include polarization orbs (as in Ti case: 3p..4p*) 524 ! (See 'archaeological note' in the header of this file) 525 if (l>0) then 526 do i = 1, basp%lshell(l-1)%nn 527 if (basp%lshell(l-1)%shell(i)%polarized) then 528 k%nkbl = k%nkbl + 1 529 write(6,"(a,i1,a)") trim(basp%label) // 530 $ ': nkbl increased for l=',l, 531 $ ' due to the presence of a polarization orbital' 532 endif 533 enddo 534 endif 535 if (k%nkbl.eq.0) then 536 write(6,*) 'Warning: Empty PAO shell. l =', l 537 write(6,*) 'Will have a KB projector anyway...' 538 k%nkbl = 1 539 endif 540 endif 541 allocate(k%erefkb(1:k%nkbl)) 542 k%erefkb(1:k%nkbl) = huge(1.d0) 543 enddo 544 else ! Set in KBprojectors 545 if (basp%lmxkb_requested.ne.-1) then ! set in PS.lmax 546 if (basp%lmxkb.ne.basp%lmxkb_requested) then 547 call die("LmaxKB conflict between " // 548 $ "PS.Lmax and PS.KBprojectors blocks") 549 endif 550 endif 551 !! OK, we have a genuine lmxkb 552 allocate(basp%kbshell(0:basp%lmxkb)) 553 do l=0,basp%lmxkb 554 call initialize(basp%kbshell(l)) 555 enddo 556 endif 557 if (basp%z .le. 0) then 558 if (basp%lmxkb .ne. -1) 559 . call die("Floating orbs cannot have KB projectors...") 560 endif 561 enddo 562! 563! Now re-scan the block (if it exists) and fill in as instructed 564! 565 if (fdf_block('PS.KBprojectors',bfdf) ) then 566 567 do while(fdf_bline(bfdf, pline)) !! over species 568 if (.not. fdf_bmatch(pline,'ni')) 569 . call die("Wrong format in PS.KBprojectors") 570 isp = label2species(fdf_bnames(pline,1)) 571 if (isp .eq. 0) then 572 write(6,'(a,1x,a)') 573 . "WRONG species symbol in PS.KBprojectors:", 574 . trim(fdf_bnames(pline,1)) 575 call die() 576 endif 577 basp => basis_parameters(isp) 578 basp%nkbshells = fdf_bintegers(pline,1) 579 do ish=1, basp%nkbshells 580 if (.not. fdf_bline(bfdf,pline)) call die("No l nkbl") 581 if (.not. fdf_bmatch(pline,'ii')) 582 . call die("Wrong format l nkbl") 583 l = fdf_bintegers(pline,1) 584 k => basp%kbshell(l) 585 k%l = l 586 k%nkbl = fdf_bintegers(pline,2) 587 if (k%nkbl < 0) then 588 call die("nkbl < 0 in PS.KBprojectors") 589 endif 590 if (k%nkbl == 0) then 591 call message("WARNING","nkbl=0 in PS.KBprojectors") 592 endif 593 allocate(k%erefkb(k%nkbl)) 594 if (.not. fdf_bline(bfdf,pline)) then 595 if (ish .ne. basp%nkbshells) 596 . call die("Not enough shells for this species...") 597 ! There is no line with ref energies 598 ! Use default values 599 k%erefKB(1:k%nkbl) = huge(1.d0) 600 else if (fdf_bmatch(pline,'ni')) then 601 ! We are seeing the next species' section 602 if (ish .ne. basp%nkbshells) 603 . call die("Not enough kb shells for this species...") 604 ! Use default values for ref energies 605 k%erefKB(1:k%nkbl) = huge(1.d0) 606 if (.not. fdf_bbackspace(bfdf)) 607 . call die('readkb: ERROR in PS.KBprojectors block') 608 else if (fdf_bmatch(pline,'ii')) then 609 ! We are seeing the next shell's section 610 if (ish .gt. basp%nkbshells) 611 . call die("Too many kb shells for this species...") 612 ! Use default values for ref energies 613 k%erefKB(1:k%nkbl) = huge(1.d0) 614 if (.not. fdf_bbackspace(bfdf)) 615 . call die('readkb: ERROR in PS.KBprojectors block') 616 else 617 if (fdf_bnreals(pline) .ne. k%nkbl) 618 . call die("Wrong number of energies") 619 unitstr = 'Ry' 620 if (fdf_bnnames(pline) .eq. 1) 621 . unitstr = fdf_bnames(pline,1) 622 ! Insert ref energies in erefkb 623 do i= 1, k%nkbl 624 k%erefKB(i) = 625 . fdf_breals(pline,i)*fdf_convfac(unitstr,'Ry') 626 enddo 627 endif 628 enddo ! end of loop over shells for species isp 629 630 ! For those l's not specified in block, use default values 631 do l=0, basp%lmxkb 632 k => basp%kbshell(l) 633 if (k%l.eq.-1) then 634 k%l = l 635 k%nkbl = 1 636 allocate(k%erefkb(1)) 637 k%erefkb(1) = huge(1.d0) 638 endif 639 enddo 640 enddo !! Over species 641 endif 642 643 do isp=1,nsp 644! 645! Check that we have enough semilocal components... 646! 647 basp=>basis_parameters(isp) 648 lmax_pseudo = basp%pseudopotential%npotd - 1 649 if (basp%lmxkb > lmax_pseudo) then 650 write(6,'(a,i1,a)') 651 . trim(basp%label) // 652 . " pseudopotential only contains V_ls up to l=", 653 . lmax_pseudo, " -- lmxkb reset." 654 basp%lmxkb = lmax_pseudo 655 endif 656 enddo 657 658 end subroutine readkb 659!--------------------------------------------------------------- 660 661 subroutine repaobasis() 662 663 integer isp, ish, nn, i, ind, l, indexp, index_splnorm 664 integer nrcs_zetas 665 666 type(block_fdf) :: bfdf 667 type(parsed_line), pointer :: pline 668 669 if (.not. fdf_block('PAO.Basis',bfdf)) RETURN 670 671 do while(fdf_bline(bfdf, pline)) !! over species 672 if (.not. fdf_bmatch(pline,'ni')) 673 . call die("Wrong format in PAO.Basis") 674 isp = label2species(fdf_bnames(pline,1)) 675 if (isp .eq. 0) then 676 write(6,'(a,1x,a)') 677 . "WRONG species symbol in PAO.Basis:", 678 . trim(fdf_bnames(pline,1)) 679 call die() 680 endif 681 682 basp => basis_parameters(isp) 683 basp%label = fdf_bnames(pline,1) 684 basp%nshells_tmp = fdf_bintegers(pline,1) 685 basp%lmxo = 0 686 !! Check whether there are optional type and ionic charge 687 if (fdf_bnnames(pline) .eq. 2) 688 . basp%basis_type = fdf_bnames(pline,2) 689 if (fdf_bnvalues(pline) .eq. 2) 690 . basp%ionic_charge = fdf_bvalues(pline,2) 691 allocate(basp%tmp_shell(basp%nshells_tmp)) 692 693 shells: do ish= 1, basp%nshells_tmp 694 s => basp%tmp_shell(ish) 695 call initialize(s) 696 if (.not. fdf_bline(bfdf,pline)) call die("No l nzeta, etc") 697 698 if (fdf_bmatch(pline,'niii')) then 699 s%n = fdf_bintegers(pline,1) 700 s%l = fdf_bintegers(pline,2) 701 basp%lmxo = max(basp%lmxo,s%l) 702 s%nzeta = fdf_bintegers(pline,3) 703 elseif (fdf_bmatch(pline,'ii')) then 704 ! l, nzeta 705 706 if (basp%semic) 707 . call die("Please specify n if there are semicore states") 708 709 s%l = fdf_bintegers(pline,1) 710 s%n = basp%ground_state%n(s%l) 711 s%nzeta = fdf_bintegers(pline,2) 712 basp%lmxo = max(basp%lmxo,s%l) 713 else 714 call die("Bad format of (n), l, nzeta line in PAO.Basis") 715 endif 716! 717! If this is a filteret basis then the number of zetas input must be one 718! 719 if (basp%basis_type.eq.'filteret') then 720 s%nzeta = 1 721 endif 722! 723! Optional stuff: Polarization, Soft-confinement Potential and Filteret Cutoff 724! 725! Split norm 726! 727 if (fdf_bsearch(pline,'S',index_splnorm)) then 728 if (fdf_bmatch(pline,'v',after=index_splnorm)) then 729 s%split_norm = fdf_bvalues(pline, ind=1, 730 . after=index_splnorm) 731 if (s%split_norm .eq. 0.0_dp) 732 . write(6,"(a)") 733 . "WARNING: zero split_norm after S in PAO.Basis" 734 s%split_norm_specified = .TRUE. 735 else 736 call die("Specify split_norm after S in PAO.Basis") 737 endif 738 else 739 if (abs(basp%z) .eq. 1) then 740 s%split_norm = global_splnorm_H 741 else 742 s%split_norm = global_splnorm 743 endif 744 endif 745! 746! Polarization functions 747! 748 if (fdf_bsearch(pline,'P',indexp)) then 749 s%polarized = .TRUE. 750 if (fdf_bmatch(pline,'i',after=indexp)) then 751 s%nzeta_pol = fdf_bintegers(pline,ind=1,after=indexp) 752 else 753 s%nzeta_pol = 1 754 endif 755 endif 756! 757! Soft-confinement 758! 759 if (fdf_bsearch(pline,'E',indexp)) then 760 if (fdf_bmatch(pline,'vv',after=indexp)) then 761 s%vcte = fdf_bvalues(pline,ind=1,after=indexp) 762 s%rinn = fdf_bvalues(pline,ind=2,after=indexp) 763 else 764 call die("Need vcte and rinn after E in PAO.Basis") 765 endif 766 elseif (lsoft) then 767 s%vcte = softPt 768 s%rinn = -softRc 769 else 770 s%vcte = 0.0_dp 771 s%rinn = 0.0_dp 772 endif 773! 774! Charge confinement 775! 776 if (fdf_bsearch(pline,'Q',indexp)) then 777 if (fdf_bmatch(pline,'vvv',after=indexp)) then 778 s%qcoe = fdf_bvalues(pline,ind=1,after=indexp) 779 s%qyuk = fdf_bvalues(pline,ind=2,after=indexp) 780 s%qwid = fdf_bvalues(pline,ind=3,after=indexp) 781 elseif (fdf_bmatch(pline,'vv',after=indexp)) then 782 s%qcoe = fdf_bvalues(pline,ind=1,after=indexp) 783 s%qyuk = fdf_bvalues(pline,ind=2,after=indexp) 784 s%qwid = 0.01_dp 785 elseif (fdf_bmatch(pline,'v',after=indexp)) then 786 s%qcoe = fdf_bvalues(pline,ind=1,after=indexp) 787 s%qyuk = 0.0_dp 788 s%qwid = 0.01_dp 789 else 790 call die("Need one, two or three real numbers after Q in 791 . PAO.Basis") 792 endif 793 else 794 s%qcoe = 0.0_dp 795 s%qyuk = 0.0_dp 796 s%qwid = 0.01_dp 797 endif 798! 799! Filteret cutoff 800! 801 if (fdf_bsearch(pline,"F",indexp)) then 802 if (fdf_bmatch(pline,"v",after=indexp)) then 803 s%filtercut = fdf_bvalues(pline,ind=1,after=indexp) 804 else 805 call die("Need cut-off after F in PAO.Basis") 806 endif 807 else 808 s%filtercut = 0.0_dp 809 endif 810 811 allocate(s%rc(s%nzeta),s%lambda(s%nzeta)) 812 s%rc(:) = 0.d0 813 s%lambda(:) = 1.d0 814 if (.not. fdf_bline(bfdf,pline)) call die("No rc's") 815 816 ! Use the last rc entered for the successive zetas 817 ! if there are not enough values (useful for Bessel) 818 nrcs_zetas = fdf_bnvalues(pline) 819 if (nrcs_zetas < 1) then 820 call die("Need at least one rc per shell in PAO.Basis block") 821 endif 822 do i= 1, s%nzeta 823 if (i <= nrcs_zetas) then 824 s%rc(i) = fdf_bvalues(pline,i) 825 else 826 s%rc(i) = s%rc(nrcs_zetas) 827 endif 828 enddo 829 830 if (s%split_norm_specified) then 831 do i = 2,s%nzeta 832 if (s%rc(i) /= 0.0_dp) then 833 write(6,"(/,a,i1,a,f8.4,/)") 834 . "*Warning: Per-shell split_norm parameter " // 835 . "will not apply to zeta-", i, ". rc=", s%rc(i) 836 endif 837 enddo 838 endif 839 840 ! Optional scale factors. They MUST be reals, or else... 841 if (.not. fdf_bline(bfdf,pline)) then 842 if (ish .ne. basp%nshells_tmp) 843 . call die("Not enough shells") 844 ! Default values for scale factors 845 else 846 if (.not. fdf_bmatch(pline,'r')) then 847 ! New shell or species 848 ! Default values for the scale factors 849 if (.not. fdf_bbackspace(bfdf)) 850 . call die('repaobasis: ERROR in PAO.Basis block') 851 cycle shells 852 else 853 ! Read scale factors 854 ! Use the last scale factor entered for the successive zetas 855 ! if there are not enough values 856 nrcs_zetas = fdf_bnreals(pline) 857 if (nrcs_zetas < 1) then 858 call die("Need at least one scale factor in PAO.Basis") 859 endif 860 do i= 1, s%nzeta 861 if (i <= nrcs_zetas) then 862 s%lambda(i) = fdf_breals(pline,i) 863 else 864 s%lambda(i) = s%lambda(nrcs_zetas) 865 endif 866 enddo 867 endif 868 endif 869 870 enddo shells 871 ! Clean up for this species 872 enddo 873! 874! OK, now classify the states by l-shells 875! 876 do isp = 1, nsp 877 basp => basis_parameters(isp) 878 if (basp%lmxo .eq. -1) cycle !! Species not in block 879 !! 880 allocate (basp%lshell(0:basp%lmxo)) 881 loop_l: do l= 0, basp%lmxo 882 ls => basp%lshell(l) 883 call initialize(ls) 884 ls%l = l 885 ! Search for tmp_shells with given l 886 nn = 0 887 do ish= 1, basp%nshells_tmp 888 s => basp%tmp_shell(ish) 889 if (s%l .eq. l) nn=nn+1 890 enddo 891 ls%nn = nn 892 if (nn.eq.0) then 893 !! What else do we do here? 894 cycle loop_l 895 endif 896 !! 897 allocate(ls%shell(1:nn)) 898 ! Collect previously allocated shells 899 ind = 0 900 do ish=1, basp%nshells_tmp 901 s => basp%tmp_shell(ish) 902 if (s%l .eq. l) then 903 ind = ind+1 904 call copy_shell(source=s,target=ls%shell(ind)) 905 endif 906 enddo 907 if (nn.eq.1) then 908 ! If n was not specified, set it to ground state n 909 if (ls%shell(1)%n.eq.-1) 910 . ls%shell(1)%n=basp%ground_state%n(l) 911 endif 912 !! Do we have to sort by n value???? 913 !! 914 enddo loop_l 915 !! Destroy temporary shells in basp 916 !! Warning: This does seem to destroy information!! 917 call destroy(basp%tmp_shell) 918 enddo 919 920 end subroutine repaobasis 921!_______________________________________________________________________ 922 923 924 function set_default_lmxkb(is) result(lmxkb) 925 integer lmxkb 926 integer, intent(in) :: is 927 928 integer lpol, l, i 929 930 lmxkb = -1 931 basp=>basis_parameters(is) 932 if (basp%z .le. 0) return ! Leave it at -1 for floating orbs. 933 934 lmxkb = basp%lmxo + 1 935 936 ! But watch out for polarization orbitals... 937 ! We ASSUME that these do not count towards lshell%nn... 938 939 lpol = 0 940 do l = 0, basp%lmxo 941 ls=>basp%lshell(l) 942 do i = 1, ls%nn 943 s=>ls%shell(i) 944 if (s%polarized) lpol = s%l + 1 945 enddo 946 enddo 947 if (lpol .gt. basp%lmxo) lmxkb = lpol + 1 948 949 write(6,'(3a,i1,/,2a,/,a)') 'For ', trim(basp%label), 950 . ', standard SIESTA heuristics set lmxkb to ', 951 . lmxkb, 952 . ' (one more than the basis l,', 953 . ' including polarization orbitals).', 954 . 'Use PS.lmax or PS.KBprojectors blocks to override.' 955! 956! But there is an upper limit for sanity: f is the highest 957! 958 if (lmxkb.gt.3) then 959 write(6,'(3a,i1)') 'Warning: For ', trim(basp%label), 960 . ' lmxkb would have been set to ', lmxkb 961 write(6,'(a)') 962 . 'Setting it to maximum value of 3 (f projector)' 963 lmxkb = 3 964 endif 965 966 end function set_default_lmxkb 967 968!----------------------------------------------------------------------- 969 970 subroutine resizes() 971 972c Reading atomic basis sizes for different species. 973c 974c Reads fdf block. Not necessarily all species have to be given. The 975c ones not given at input will be assumed to have the basis sizes 976c given by the general input PAO.BasisSize, or its default value. 977 978 type(block_fdf) :: bfdf 979 type(parsed_line), pointer :: pline 980 981 integer isp 982 983 if (fdf_block('PAO.BasisSizes',bfdf)) then 984 do while(fdf_bline(bfdf,pline)) 985 if (.not. fdf_bmatch(pline,'nn')) 986 . call die("Wrong format in PAO.BasisSizes") 987 isp = label2species(fdf_bnames(pline,1)) 988 if (isp .eq. 0) then 989 write(6,'(a,1x,a)') 990 . "WRONG species symbol in PAO.BasisSizes:", 991 . trim(fdf_bnames(pline,1)) 992 call die() 993 else 994 basp => basis_parameters(isp) 995 basp%basis_size = fdf_bnames(pline,2) 996 call size_name(basp%basis_size) !!! DEPRECATED 997 write(6,'(4a)') 998 . 'resizes: Read basis size for species ', 999 . trim(fdf_bnames(pline,1)),' = ',basp%basis_size 1000 endif 1001 enddo 1002 endif 1003 1004 end subroutine resizes 1005 1006!----------------------------------------------------------------------- 1007 1008 subroutine relmxkb() 1009 1010c Reads the maximum angular momentum of the Kleinman-Bylander 1011c projectors for the different species. 1012c 1013c Reads fdf block. Not necessarily all species have to be given. 1014 1015 type(block_fdf) :: bfdf 1016 type(parsed_line), pointer :: pline 1017 1018 integer isp 1019 1020 if (fdf_block('PS.lmax',bfdf)) then 1021 do while(fdf_bline(bfdf,pline)) 1022 if (.not. fdf_bmatch(pline,'ni')) 1023 . call die("Wrong format in PS.lmax") 1024 isp = label2species(fdf_bnames(pline,1)) 1025 if (isp .eq. 0) then 1026 write(6,'(a,1x,a)') "WRONG species symbol in PS.lmax:", 1027 . trim(fdf_bnames(pline,1)) 1028 call die() 1029 else 1030 basp => basis_parameters(isp) 1031 basp%lmxkb_requested = fdf_bintegers(pline,1) 1032 write(6,"(a, i4, 2a)") 1033 . 'relmxkb: Read Max KB Ang. Momentum= ', 1034 . basp%lmxkb_requested, 1035 . ' for species ', trim(fdf_bnames(pline,1)) 1036 endif 1037 enddo 1038 endif 1039 1040 end subroutine relmxkb 1041!----------------------------------------------------------------------- 1042 1043 subroutine remass() 1044 1045 1046 type(block_fdf) :: bfdf 1047 type(parsed_line), pointer :: pline 1048 1049 integer isp 1050 1051c Read atomic masses of different species. 1052c 1053c Reads fdf block. Not necessarily all species have to be given. The 1054c ones not given at input will be assumed to have their natural mass 1055c (according to atmass subroutine). 1056 1057 if (fdf_block('AtomicMass',bfdf)) then 1058 do while(fdf_bline(bfdf,pline)) 1059 if (.not. fdf_bmatch(pline,'iv')) 1060 . call die("Wrong format in AtomicMass") 1061 isp = fdf_bintegers(pline,1) 1062 if (isp .gt. nsp .or. isp .lt. 1) 1063 . call die("Wrong specnum in AtomicMass") 1064 basp => basis_parameters(isp) 1065 basp%mass = fdf_bvalues(pline,2) 1066 write(6,"(a, i4, a, f12.5)") 1067 . 'remass: Read atomic mass for species ', isp, 1068 . ' as ', basp%mass 1069 enddo 1070 endif 1071 1072 end subroutine remass 1073!----------------------------------------------------------------------- 1074 1075 subroutine size_name(str) 1076 character(len=*), intent(inout) :: str 1077 1078 1079 if(leqi(str,'MINIMAL')) str='sz' 1080 if(leqi(str,'SZ')) str='sz' 1081 if(leqi(str,'SZP')) str='szp' 1082 if(leqi(str,'SZP1')) str='szp' 1083 if(leqi(str,'SZSP')) str='szp' 1084 if(leqi(str,'SZ1P')) str='szp' 1085! 1086 if(leqi(str,'DZ')) str='dz' 1087 if(leqi(str,'STANDARD')) str='dzp' 1088 if(leqi(str,'DZP')) str='dzp' 1089 if(leqi(str,'DZP1')) str='dzp' 1090 if(leqi(str,'DZ1P')) str='dzp' 1091 if(leqi(str,'DZSP')) str='dzp' 1092 if(leqi(str,'DZP2')) str='dzp2' 1093 if(leqi(str,'DZDP')) str='dzp2' 1094 if(leqi(str,'DZ2P')) str='dzp2' 1095! 1096 if(leqi(str,'TZ')) str='tz' 1097 if(leqi(str,'TZP')) str='tzp' 1098 if(leqi(str,'TZ1P')) str='tzp' 1099 if(leqi(str,'TZP1')) str='tzp' 1100 if(leqi(str,'TZSP')) str='tzp' 1101 if(leqi(str,'TZP2')) str='tzp2' 1102 if(leqi(str,'TZ2P')) str='tzp2' 1103 if(leqi(str,'TZDP')) str='tzp2' 1104 if(leqi(str,'TZP3')) str='tzp3' 1105 if(leqi(str,'TZ3P')) str='tzp3' 1106 if(leqi(str,'TZTP')) str='tzp3' 1107 1108 if ( (str.ne.'szp').and.(str.ne.'sz').and. 1109 . (str.ne.'dz') .and.(str.ne.'dzp') .and. 1110 . (str.ne.'tz') .and.(str.ne.'tzp') .and. 1111 . (str.ne.'dzp2') .and. 1112 . (str.ne.'tzp2') .and. (str.ne.'tzp3') ) then 1113 1114 write(6,'(/,2a,/,9(a,/))') 1115 . 'size_name: Incorrect basis-size option specified,', 1116 . ' active options are:', 1117 . ' SZ or MINIMAL', 1118 . ' SZP, SZSP, SZ1P, SZP1', 1119 . ' DZ ', 1120 . ' DZP, DZSP, DZP1, DZ1P or STANDARD', 1121 . ' DZDP, DZP2, DZ2P ', 1122 . ' TZ ', 1123 . ' TZP, TZSP, TZP1, TZ1P', 1124 . ' TZDP, TZP2, TZ2P', 1125 . ' TZTP, TZP3, TZ3P' 1126 1127 call die() 1128 endif 1129 1130 end subroutine size_name 1131!----------------------------------------------------------------------- 1132 1133 subroutine type_name(basistype) 1134 1135 character basistype*(*) 1136 1137 if(leqi(basistype,'NODES')) then 1138 basistype='nodes' 1139 elseif(leqi(basistype,'NONODES')) then 1140 basistype='nonodes' 1141 elseif(leqi(basistype,'SPLIT')) then 1142 basistype='split' 1143 elseif(leqi(basistype,'SPLITGAUSS')) then 1144 basistype='splitgauss' 1145 elseif (leqi(basistype,'FILTERET')) then 1146 basistype='filteret' 1147 else 1148 write(6,'(/,2a,(/,5(3x,a)),(/,2(3x,a)))') 1149 . 'type_name: Incorrect basis-type option specified,', 1150 . ' active options are:', 1151 . 'NODES','SPLIT','SPLITGAUSS','NONODES','FILTERET' 1152 call die 1153 endif 1154 1155 end subroutine type_name 1156!----------------------------------------------------------------------- 1157! 1158! Find out whether semicore states are implied by the valence 1159! charge density in the pseudopotential file. 1160! 1161 subroutine semicore_check(is) 1162 integer, intent(in) :: is 1163 1164 real(dp), parameter :: tiny = 1.d-5 1165 integer ndiff 1166 real(dp) zval, zval_vps, charge_loc 1167 1168 basp => basis_parameters(is) 1169 1170 basp%semic = .false. 1171 if (basp%bessel) return 1172 1173 zval_vps = basp%pseudopotential%zval 1174 zval = basp%ground_state%z_valence 1175 1176 if (abs(Zval-zval_vps).lt.tiny) return 1177 1178 ndiff = nint(abs(Zval-zval_vps)) 1179 if (abs(ndiff-abs(Zval-zval_vps)).gt.tiny) then 1180 write(6,'(2a)') 1181 . 'ERROR expected valence charge for species ', 1182 . basp%label 1183 write(6,'(a)') 1184 . 'ERROR and the value read from the vps file' 1185 write(6,'(a,f6.3,a,f6.3)') 1186 . 'ERROR differ: Zval(expected)= ', Zval, 1187 . ' Zval(vps)= ',zval_vps 1188 call die() 1189 endif 1190 1191 basp%semic = .true. 1192 charge_loc = Zval_vps-Zval 1193 write(6,'(a,i2,a)') 1194 . 'Semicore shell(s) with ', nint(charge_loc), 1195 . ' electrons included in the valence for', trim(basp%label) 1196 1197 end subroutine semicore_check 1198!---------------------------------------------------------------------- 1199 subroutine autobasis() 1200 1201! 1202! It sets the defaults if a species has not been included 1203! in the PAO.Basis block 1204! 1205 integer l, nzeta, nzeta_pol 1206 1207 loop: do isp=1, nsp 1208 basp=>basis_parameters(isp) 1209 if (basp%lmxo .ne. -1) cycle loop ! Species already set 1210 ! in PAO.Basis block 1211 if (basp%semic) then 1212 write(6,'(2a)') basp%label, 1213 . ' must be in PAO.Basis (it has semicore states)' 1214 call die() 1215 endif 1216 if (basp%bessel) then 1217 write(6,'(2a)') basp%label, 1218 . ' must be in PAO.Basis (it is a floating Bessel function)' 1219 call die() 1220 endif 1221 ! 1222 ! Set the default max l 1223 ! 1224 basp%lmxo = basp%ground_state%lmax_valence 1225 1226 allocate (basp%lshell(0:basp%lmxo)) 1227 1228 if (basp%basis_size(1:2) .eq. 'sz') nzeta = 1 1229 if (basp%basis_size(1:2) .eq. 'dz') nzeta = 2 1230 if (basp%basis_size(1:2) .eq. 'tz') nzeta = 3 1231 1232 loop_l: do l=0, basp%lmxo 1233 ls=>basp%lshell(l) 1234 call initialize(ls) 1235 ls%l = l 1236 ls%nn = 1 1237 allocate(ls%shell(1:1)) 1238 s => ls%shell(1) 1239 call initialize(s) 1240 s%l = l 1241 s%n = basp%ground_state%n(l) 1242 if (basp%ground_state%occupied(l)) then 1243 s%nzeta = nzeta 1244 else 1245 s%nzeta = 0 1246 endif 1247 s%polarized = .false. 1248 if (abs(basp%z).eq.1) then 1249 s%split_norm = global_splnorm_H 1250 else 1251 s%split_norm = global_splnorm 1252 endif 1253 s%nzeta_pol = 0 1254 1255 if (lsoft) then 1256 s%vcte = softPt 1257 s%rinn = -softRc 1258 else 1259 s%rinn = 0.d0 1260 s%vcte = 0.d0 1261 endif 1262 1263 ! Default filteret cutoff for shell 1264 s%filtercut = 0.0d0 1265 1266 if (s%nzeta .ne.0) then 1267 allocate(s%rc(1:s%nzeta)) 1268 allocate(s%lambda(1:s%nzeta)) 1269 s%rc(1:s%nzeta) = 0.0d0 1270 s%lambda(1:s%nzeta) = 1.0d0 1271 endif 1272 enddo loop_l 1273 1274 if (basp%basis_size(3:3) .eq. 'p') then 1275 1276 ! Polarization orbitals are defined so they have the minimum 1277 ! angular momentum l such that there are not occupied orbitals 1278 ! with the same l in the valence shell of the ground-state 1279 ! atomic configuration. They polarize the corresponding l-1 1280 ! shell. 1281 1282 ! note that we go up to l=4 to make the loop simpler. 1283 ! (that is the reason why 'occupied' is dimensioned to 0:4) 1284 1285 select case (basp%basis_size(4:4)) 1286 case (' ', '1') 1287 nzeta_pol = 1 1288 case ('2') 1289 nzeta_pol = 2 1290 case ('3') 1291 nzeta_pol = 3 1292 end select 1293 1294 loop_angmom: do l=1,4 1295 if (.not. basp%ground_state%occupied(l)) then 1296 ls=>basp%lshell(l-1) 1297 s => ls%shell(1) 1298 1299 ! Check whether shell to be polarized is occupied in the gs. 1300 ! (i.e., whether PAOs are going to be generated for it) 1301 ! If it is not, mark it for PAO generation anyway. 1302 ! This will happen for confs of the type s0 p0 dn 1303 1304 ! Default filter cutoff for shell 1305 s%filtercut = 0.0d0 1306 1307 if (s%nzeta == 0) then 1308 write(6,"(a,i2,a)") "Marking shell with l=", 1309 . l-1, " for PAO generation and polarization." 1310 s%nzeta = nzeta 1311 allocate(s%rc(1:s%nzeta)) 1312 allocate(s%lambda(1:s%nzeta)) 1313 s%rc(1:s%nzeta) = 0.0d0 1314 s%lambda(1:s%nzeta) = 1.0d0 1315 endif 1316 s%polarized = .true. 1317 s%nzeta_pol = nzeta_pol 1318 1319 exit loop_angmom ! Polarize only one shell! 1320 endif 1321 enddo loop_angmom 1322 1323 endif 1324 1325 enddo loop 1326 1327 end subroutine autobasis 1328 1329 End module basis_specs 1330