Octopus
ps.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2012 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module ps_oct_m
22 use atomic_oct_m
23 use debug_oct_m
24 use global_oct_m
25 use io_oct_m
26 use, intrinsic :: iso_fortran_env
30 use parser_oct_m
31 use logrid_oct_m
32 use math_oct_m
36 use ps_cpi_oct_m
37 use ps_fhi_oct_m
38 use ps_hgh_oct_m
39 use ps_xml_oct_m
41#ifdef HAVE_PSPIO
42 use fpspio_m
43#endif
44 use ps_psf_oct_m
45 use pseudo_oct_m
46 use sort_oct_m
49 use string_oct_m
50
51 implicit none
52
53 private
54 public :: &
55 ps_t, &
56 ps_init, &
59 ps_filter, &
62 ps_debug, &
63 ps_niwfs, &
65 ps_end, &
70
71 integer, parameter, public :: &
72 PS_FILTER_NONE = 0, &
73 ps_filter_ts = 2, &
75
76 integer, public, parameter :: &
77 PROJ_NONE = 0, &
78 proj_hgh = 1, &
79 proj_kb = 2, &
80 proj_rkb = 3
81
82 integer, public, parameter :: &
83 PROJ_J_INDEPENDENT = 0, & !< Non-relativisitic or scalar-relativistic pseudopotentials
84 proj_j_dependent = 1, &
86
87 integer, parameter, public :: INVALID_L = 333
88
89 character(len=4), parameter :: ps_name(PSEUDO_FORMAT_UPF1:PSEUDO_FORMAT_PSP8) = &
90 (/"upf1", "upf2", "qso ", "psml", "psf ", "cpi ", "fhi ", "hgh ", "psp8"/)
91
93 type ps_t
94 ! Components are public by default
95 integer :: projector_type
96 integer :: relativistic_treatment
98
99 character(len=10), private :: label
100
101 integer, private :: ispin
102
103 real(real64), private :: z
104 real(real64) :: z_val
105 type(valconf_t) :: conf
106
107 type(logrid_t), private :: g
108
109 type(spline_t), allocatable :: ur(:, :)
110 logical, allocatable :: bound(:, :)
111
112 ! Kleinman-Bylander projectors stuff
113 integer :: lmax
114 integer :: llocal
115
116 type(spline_t) :: vl
117 logical :: no_vl = .false.
118
119 real(real64) :: projectors_sphere_threshold
126 real(real64) :: rc_max
127
128 integer :: kbc
129 integer :: projectors_per_l(5)
130 real(real64), allocatable :: h(:,:,:), k(:, :, :)
131 type(spline_t), allocatable :: kb(:, :)
132 type(spline_t), allocatable :: dkb(:, :)
133
134 logical :: nlcc
135 type(spline_t) :: core
136 type(spline_t) :: core_der
137
138
139 !LONG-RANGE PART OF THE LOCAL POTENTIAL
140
141 logical, private :: has_long_range
142 logical, private :: is_separated
143
144 type(spline_t), private :: vlr
145 type(spline_t) :: vlr_sq
147 type(spline_t) :: nlr
148
149 real(real64) :: sigma_erf
150
151 logical, private :: has_density
152 type(spline_t), allocatable :: density(:)
153 type(spline_t), allocatable :: density_der(:)
154
155 logical :: local
156 integer :: file_format
157 integer, private :: pseudo_type
158 integer :: exchange_functional
159 integer :: correlation_functional
160 end type ps_t
161
162 real(real64), parameter :: eps = 1.0e-8_real64
163
164contains
165
167 ! ---------------------------------------------------------
168 subroutine ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
169 type(ps_t), intent(out) :: ps
170 type(namespace_t), intent(in) :: namespace
171 character(len=10), intent(in) :: label
172 integer, intent(in) :: user_lmax
173 integer, intent(in) :: user_llocal
174 integer, intent(in) :: ispin
175 real(real64), intent(in) :: z
176 character(len=*), intent(in) :: filename
178 integer :: l, ii, ll, is, ierr
179 type(ps_psf_t) :: ps_psf
180 type(ps_cpi_t) :: ps_cpi
181 type(ps_fhi_t) :: ps_fhi
182 type(ps_hgh_t) :: ps_hgh
183 type(ps_xml_t) :: ps_xml
184 real(real64), allocatable :: eigen(:, :)
185
186 push_sub(ps_init)
187
188 ps%exchange_functional = pseudo_exchange_unknown
189 ps%correlation_functional = pseudo_correlation_unknown
191 ! Fix the threshold to calculate the radius of the projector-function localization spheres:
192
193 call messages_obsolete_variable(namespace, 'SpecieProjectorSphereThreshold', 'SpeciesProjectorSphereThreshold')
195 !%Variable SpeciesProjectorSphereThreshold
196 !%Type float
197 !%Default 0.001
198 !%Section System::Species
199 !%Description
200 !% The pseudopotentials may be composed of a local part, and a linear combination of nonlocal
201 !% operators. These nonlocal projectors have "projector" form, <math> \left| v \right> \left< v \right| </math>
202 !% (or, more generally speaking, <math> \left| u \right> \left< v \right| </math>).
203 !% These projectors are localized in real space -- that is, the function <math>v</math>
204 !% has a finite support around the nucleus. This region where the projectors are localized should
205 !% be small or else the computation time required to operate with them will be very large.
206 !%
207 !% In practice, this localization is fixed by requiring the definition of the projectors to be
208 !% contained in a sphere of a certain radius. This radius is computed by making sure that the
209 !% absolute value of the projector functions, at points outside the localization sphere, is
210 !% below a certain threshold. This threshold is set by <tt>SpeciesProjectorSphereThreshold</tt>.
211 !%End
212 call parse_variable(namespace, 'SpeciesProjectorSphereThreshold', 0.001_real64, ps%projectors_sphere_threshold)
213 if (ps%projectors_sphere_threshold <= m_zero) call messages_input_error(namespace, 'SpeciesProjectorSphereThreshold')
215 ps%file_format = pseudo_detect_format(string_f_to_c(filename))
216
217 if (ps%file_format == pseudo_format_file_not_found) then
218 call messages_write("Cannot open pseudopotential file '"//trim(filename)//"'.")
219 call messages_fatal(namespace=namespace)
220 end if
222 if (ps%file_format == pseudo_format_unknown) then
223 call messages_write("Cannot determine the pseudopotential type for species '"//trim(label)//"' from", &
224 new_line = .true.)
225 call messages_write("file '"//trim(filename)//"'.")
226 call messages_fatal(namespace=namespace)
227 end if
228
229 ps%label = label
230 ps%ispin = ispin
231 ps%relativistic_treatment = proj_j_independent
232 ps%projector_type = proj_kb
233 ps%sigma_erf = 0.625_real64 ! This is hard-coded to a reasonable value
234
235 ps%projectors_per_l = 0
237 select case (ps%file_format)
239 ps%has_density = .true.
240 case default
241 ps%has_density = .false.
242 end select
243
244 select case (ps%file_format)
245 case (pseudo_format_psf)
246 ps%pseudo_type = pseudo_type_semilocal
248 call ps_psf_init(ps_psf, ispin, filename, namespace)
249
250 call valconf_copy(ps%conf, ps_psf%conf)
251 ps%z = z
252 ps%conf%z = nint(z) ! atomic number
253 ps%kbc = 1 ! only one projector per angular momentum
255 ps%lmax = ps_psf%ps_grid%no_l_channels - 1
256
257 if (user_lmax /= invalid_l) then
258 ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components.
259 if (user_lmax /= ps%lmax) then
260 message(1) = "lmax in Species block for " // trim(label) // &
261 " is larger than number available in pseudopotential."
262 call messages_fatal(1, namespace=namespace)
263 end if
264 end if
265
266 ps%conf%p = ps_psf%ps_grid%no_l_channels
267 if (ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0.
268
269 ! the local part of the pseudo
270 if (user_llocal == invalid_l) then
271 ps%llocal = 0
272 else
273 ps%llocal = user_llocal
274 end if
275
276 ps%projectors_per_l(1:ps%lmax+1) = 1
277 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
278
279 call ps_psf_process(ps_psf, namespace, ps%lmax, ps%llocal)
280 call logrid_copy(ps_psf%ps_grid%g, ps%g)
281
283 ps%pseudo_type = pseudo_type_semilocal
284
285 if (ps%file_format == pseudo_format_cpi) then
286 call ps_cpi_init(ps_cpi, trim(filename), namespace)
287 ps%conf%p = ps_cpi%ps_grid%no_l_channels
288 else
289 call ps_fhi_init(ps_fhi, trim(filename), namespace)
290 ps%conf%p = ps_fhi%ps_grid%no_l_channels
291 end if
292
293 ps%conf%z = nint(z)
294 ps%conf%symbol = label(1:2)
295 ps%conf%type = 1
296 do l = 1, ps%conf%p
297 ps%conf%l(l) = l - 1
298 end do
299
300 ps%z = z
301 ps%kbc = 1 ! only one projector per angular momentum
302
303 ps%lmax = ps%conf%p - 1
304
305 if (user_lmax /= invalid_l) then
306 ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components.
307 if (user_lmax /= ps%lmax) then
308 message(1) = "lmax in Species block for " // trim(label) // &
309 " is larger than number available in pseudopotential."
310 call messages_fatal(1, namespace=namespace)
311 end if
312 end if
313
314 if (ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0.
315
316 ! the local part of the pseudo
317 if (user_llocal == invalid_l) then
318 ps%llocal = 0
319 else
320 ps%llocal = user_llocal
321 end if
322
323 ps%projectors_per_l(1:ps%lmax+1) = 1
324 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
325
326 if (ps%file_format == pseudo_format_cpi) then
327 call ps_cpi_process(ps_cpi, ps%llocal, namespace)
328 call logrid_copy(ps_cpi%ps_grid%g, ps%g)
329 else
330 call ps_fhi_process(ps_fhi, ps%lmax, ps%llocal, namespace)
331 call logrid_copy(ps_fhi%ps_grid%g, ps%g)
332 end if
333
334 case (pseudo_format_hgh)
335 ps%pseudo_type = pseudo_type_kleinman_bylander
336 ps%projector_type = proj_hgh
337
338 call hgh_init(ps_hgh, trim(filename), namespace)
339 call valconf_copy(ps%conf, ps_hgh%conf)
340
341 ps%z = z
342 ps%z_val = ps_hgh%z_val
343 ps%kbc = 3
344 ps%llocal = -1
345 ps%lmax = ps_hgh%l_max
346 ps%conf%symbol = label(1:2)
347 ps%sigma_erf = ps_hgh%rlocal ! We use the correct value
348
349 ps%projectors_per_l(1:ps%lmax+1) = 1
350
351 ! Get the occupations from the valence charge of the atom
352 ps%conf%occ = m_zero
353 ! We impose here non-spin-polarized occupations, to preserve the behavior of the code
354 ! We might want to change this to get a better LCAO guess
355 call ps_guess_atomic_occupations(namespace, ps%z, ps%z_val, 1, ps%conf)
356 ! We need the information to solve the Schrodinder equation
357 call valconf_copy(ps_hgh%conf, ps%conf)
358
359 call hgh_process(ps_hgh, namespace)
360 call logrid_copy(ps_hgh%g, ps%g)
361
362 ! In case of spin-polarized calculations, we properly distribute the electrons
363 if (ispin == 2) then
365 end if
366
367 case (pseudo_format_qso, pseudo_format_upf1, pseudo_format_upf2, pseudo_format_psml, pseudo_format_psp8)
368
369 call ps_xml_init(ps_xml, namespace, trim(filename), ps%file_format, ierr)
370
371 ps%pseudo_type = pseudo_type(ps_xml%pseudo)
372 ps%exchange_functional = pseudo_exchange(ps_xml%pseudo)
373 ps%correlation_functional = pseudo_correlation(ps_xml%pseudo)
374
375 ps%z = z
376 ps%conf%z = nint(z)
377
378 if (ps_xml%kleinman_bylander) then
379 ps%conf%p = ps_xml%nwavefunctions
380 else
381 ps%conf%p = ps_xml%lmax + 1
382 end if
383
384 do ll = 0, ps_xml%lmax
385 ps%conf%l(ll + 1) = ll
386 end do
387
388 ps%kbc = ps_xml%nchannels
389 ps%lmax = ps_xml%lmax
390
391 ps%projectors_per_l = 0
392 do ll = 0, ps_xml%lmax
393 ps%projectors_per_l(ll+1) = pseudo_nprojectors_per_l(ps_xml%pseudo, ll)
394 end do
395
396 if (ps_xml%kleinman_bylander) then
397 ps%llocal = ps_xml%llocal
398 else
399 ! we have several options
400 ps%llocal = 0 ! the default
401 if (ps_xml%llocal >= 0) ps%llocal = ps_xml%llocal ! the one given in the pseudopotential file
402 if (user_llocal /= invalid_l) ps%llocal = user_llocal ! user supplied local component
403 assert(ps%llocal >= 0)
404 assert(ps%llocal <= ps%lmax)
405 end if
406
407 ps%g%nrval = ps_xml%grid_size
408
409 safe_allocate(ps%g%rofi(1:ps%g%nrval))
410 safe_allocate(ps%g%r2ofi(1:ps%g%nrval))
411
412 do ii = 1, ps%g%nrval
413 ps%g%rofi(ii) = ps_xml%grid(ii)
414 ps%g%r2ofi(ii) = ps_xml%grid(ii)**2
415 end do
416
417 end select
418
419 ps%local = (ps%lmax == 0 .and. ps%llocal == 0 ) .or. (ps%lmax == -1 .and. ps%llocal == -1)
420
421 ! We allocate all the stuff
422 safe_allocate(ps%kb (0:ps%lmax, 1:ps%kbc))
423 safe_allocate(ps%dkb (0:ps%lmax, 1:ps%kbc))
424 safe_allocate(ps%ur (1:ps%conf%p, 1:ps%ispin))
425 safe_allocate(ps%bound(1:ps%conf%p, 1:ps%ispin))
426 safe_allocate(ps%h (0:ps%lmax, 1:ps%kbc, 1:ps%kbc))
427 safe_allocate(ps%density(1:ps%ispin))
428 safe_allocate(ps%density_der(1:ps%ispin))
429
430 call spline_init(ps%kb)
431 call spline_init(ps%dkb)
432 call spline_init(ps%vl)
433 call spline_init(ps%core)
434 call spline_init(ps%core_der)
435 call spline_init(ps%density)
436 call spline_init(ps%density_der)
437
438 safe_allocate(eigen(1:ps%conf%p, 1:ps%ispin))
439 eigen = m_zero
440
441 ! Now we load the necessary information.
442 select case (ps%file_format)
443 case (pseudo_format_psf)
444 call ps_psf_get_eigen(ps_psf, eigen)
445 call ps_grid_load(ps, ps_psf%ps_grid)
446 call ps_psf_end(ps_psf)
447 case (pseudo_format_cpi)
448 call ps_grid_load(ps, ps_cpi%ps_grid)
449 call ps_cpi_end(ps_cpi)
450 case (pseudo_format_fhi)
451 call ps_grid_load(ps, ps_fhi%ps_grid)
452 call ps_fhi_end(ps_fhi)
453 case (pseudo_format_hgh)
454 call hgh_get_eigen(ps_hgh, eigen)
455 safe_allocate(ps%k (0:ps%lmax, 1:ps%kbc, 1:ps%kbc))
456 call hgh_load(ps, ps_hgh)
457 call hgh_end(ps_hgh)
458 case (pseudo_format_qso, pseudo_format_upf1, pseudo_format_upf2, pseudo_format_psml, pseudo_format_psp8)
459 call ps_xml_load(ps, ps_xml, namespace)
460 call ps_xml_end(ps_xml)
461 end select
462
463 if (ps_has_density(ps)) then
464 do is = 1, ps%ispin
465 call spline_der(ps%density(is), ps%density_der(is), ps%projectors_sphere_threshold)
466 end do
467 end if
468
469 if (ps_has_nlcc(ps)) then
470 call spline_der(ps%core, ps%core_der, ps%projectors_sphere_threshold)
471 end if
472
473 call ps_check_bound(ps, eigen)
474
475 ps%has_long_range = .true.
476 ps%is_separated = .false.
477
478 call ps_info(ps, filename, namespace)
479
480 safe_deallocate_a(eigen)
481
482 pop_sub(ps_init)
483 end subroutine ps_init
484
485 !------------------------------------------------------------------------
486
487 subroutine ps_info(ps, filename, namespace)
488 type(ps_t), intent(in) :: ps
489 character(len=*), intent(in) :: filename
490 type(namespace_t), intent(in) :: namespace
491
492 call messages_write(" Species '"//trim(ps%label)//"'", new_line = .true.)
493 call messages_write(" type : pseudopotential", new_line = .true.)
494 call messages_write(" file : '"//trim(filename)//"'")
495 call messages_info(namespace=namespace)
496
497 call messages_write(" file format :")
498 select case (ps%file_format)
499 case (pseudo_format_upf1)
500 call messages_write(" UPF1")
501 case (pseudo_format_upf2)
502 call messages_write(" UPF2")
503 case (pseudo_format_qso)
504 call messages_write(" QSO")
505 case (pseudo_format_psml)
506 call messages_write(" PSML")
507 case (pseudo_format_psp8)
508 call messages_write(" PSP8")
509 case (pseudo_format_psf)
510 call messages_write(" PSF")
511 case (pseudo_format_cpi)
512 call messages_write(" CPI")
513 case (pseudo_format_fhi)
514 call messages_write(" FHI")
515 case (pseudo_format_hgh)
516 call messages_write(" HGH")
517 end select
518 call messages_new_line()
519
520 call messages_write(" valence charge :")
521 call messages_write(ps%z_val, align_left = .true., fmt = '(f4.1)')
522 call messages_info(namespace=namespace)
523
524 call messages_write(" atomic number :")
525 call messages_write(nint(ps%z), fmt = '(i4)')
526 call messages_info(namespace=namespace)
527
528 call messages_write(" form on file :")
529 select case (ps%pseudo_type)
531 call messages_write(" ultrasoft")
533 call messages_write(" semilocal")
535 call messages_write(" kleinman-bylander")
536 case (pseudo_type_paw)
537 call messages_write(" paw")
538 end select
539 call messages_info(namespace=namespace)
540
541 if (ps%pseudo_type == pseudo_type_semilocal) then
542 call messages_write(" orbital origin :")
543 select case (ps%file_format)
544 case (pseudo_format_psf)
545 call messages_write(" calculated")
546 case default
547 call messages_write(" from file")
548 end select
549 call messages_info(namespace=namespace)
550 end if
551
552 call messages_write(" lmax :")
553 call messages_write(ps%lmax, fmt = '(i2)')
554 call messages_info(namespace=namespace)
555
556 call messages_write(" llocal :")
557 if (ps%llocal >= 0) then
558 call messages_write(ps%llocal, fmt = '(i2)')
559 else
560 call messages_write(ps%llocal, fmt = '(i3)')
561 end if
562 call messages_info(namespace=namespace)
563
564 call messages_write(" projectors per l :")
565 call messages_write(ps%kbc, fmt = '(i2)')
566 call messages_info(namespace=namespace)
567
568 call messages_write(" total projectors :")
569 if (ps%llocal < 0) then
570 call messages_write(ps%kbc*(ps%lmax + 1), fmt = '(i2)')
571 else
572 call messages_write(ps%kbc*ps%lmax, fmt = '(i2)')
573 end if
574 call messages_info(namespace=namespace)
575
576 if (ps%local) then
577 call messages_write(" application form : local")
578 else
579 call messages_write(" application form : kleinman-bylander")
580 end if
581 call messages_info(namespace=namespace)
583 call messages_write(" orbitals :")
584 call messages_write(ps_niwfs(ps), fmt='(i3)')
585 call messages_info(namespace=namespace)
586 call messages_write(" bound orbitals :")
587 call messages_write(ps_bound_niwfs(ps), fmt='(i3)')
588 call messages_info(namespace=namespace)
589
590 call messages_info(namespace=namespace)
591
592 end subroutine ps_info
593
594
595 ! ---------------------------------------------------------
598 subroutine ps_separate(ps)
599 type(ps_t), intent(inout) :: ps
600
601 real(real64), allocatable :: vsr(:), vlr(:), nlr(:)
602 real(real64) :: r, exp_arg, arg, max_val_vsr
603 integer :: ii
604
605 push_sub(ps_separate)
606
607 assert(ps%g%nrval > 0)
608
609 safe_allocate(vsr(1:ps%g%nrval))
610 safe_allocate(vlr(1:ps%g%nrval))
611 safe_allocate(nlr(1:ps%g%nrval))
612
613 ps%no_vl = .false.
614
615 max_val_vsr = m_zero
616
617 do ii = 1, ps%g%nrval
618 r = ps%g%rofi(ii)
619 arg = r/(ps%sigma_erf*sqrt(m_two))
620 ! Note that the threshold is taken such that the Taylor expansion produces an error smaller than \epsilon
621 ! given a third order Taylor expansion
622 if(arg > 1e-3_real64) then
623 vlr(ii) = -ps%z_val*loct_erf(arg)/r
624 else
625 vlr(ii) = -ps%z_val*m_two/(ps%sigma_erf*sqrt(m_two*m_pi)) * (m_one - arg**2 / m_three)
626 end if
627
628 vsr(ii) = spline_eval(ps%vl, r) - vlr(ii)
629 if(abs(vsr(ii)) < 1.0e-14_real64) vsr(ii) = m_zero
630 max_val_vsr = max(max_val_vsr, abs(vsr(ii)))
631
632 exp_arg = -m_half*r**2/ps%sigma_erf**2
633 if (exp_arg > m_min_exp_arg) then
634 nlr(ii) = -ps%z_val/(ps%sigma_erf*sqrt(m_two*m_pi))**3*exp(exp_arg)
635 else
636 nlr(ii) = m_zero
637 end if
638 end do
639
640 call spline_init(ps%vlr)
641 call spline_fit(ps%g%nrval, ps%g%rofi, vlr, ps%vlr, ps%projectors_sphere_threshold)
642
643 call spline_init(ps%vlr_sq)
644 call spline_fit(ps%g%nrval, ps%g%r2ofi, vlr, ps%vlr_sq, ps%projectors_sphere_threshold)
645
646 call spline_init(ps%nlr)
647 call spline_fit(ps%g%nrval, ps%g%rofi, nlr, ps%nlr, ps%projectors_sphere_threshold)
648
649 !overwrite vl
650 call spline_end(ps%vl)
651 call spline_init(ps%vl)
652 call spline_fit(ps%g%nrval, ps%g%rofi, vsr, ps%vl, ps%projectors_sphere_threshold)
653 if (max_val_vsr < 1.0e-12_real64) ps%no_vl = .true.
654
655 safe_deallocate_a(vsr)
656 safe_deallocate_a(vlr)
657 safe_deallocate_a(nlr)
658
659 ps%is_separated = .true.
660
661 pop_sub(ps_separate)
662 end subroutine ps_separate
663
664
665 ! ---------------------------------------------------------
666 subroutine ps_getradius(ps)
667 type(ps_t), intent(inout) :: ps
668 integer :: l, j
669
670 push_sub(ps_getradius)
671
672 ps%rc_max = m_zero
673
674 do l = 0, ps%lmax
675 if (l == ps%llocal) cycle
676 do j = 1, ps%kbc
677 ps%rc_max = max(ps%rc_max, ps%kb(l, j)%x_threshold)
678 end do
679 end do
680
681 pop_sub(ps_getradius)
682 end subroutine ps_getradius
683
684
685 ! ---------------------------------------------------------
686 subroutine ps_derivatives(ps)
687 type(ps_t), intent(inout) :: ps
688 integer :: l, j
689
690 push_sub(ps_derivatives)
691
692 do l = 0, ps%lmax
693 do j = 1, ps%kbc
694 call spline_der(ps%kb(l, j), ps%dkb(l, j), ps%projectors_sphere_threshold)
695 end do
696 end do
697
698 pop_sub(ps_derivatives)
699 end subroutine ps_derivatives
700
701
702 ! ---------------------------------------------------------
703 subroutine ps_filter(ps, filter, gmax)
704 type(ps_t), intent(inout) :: ps
705 integer, intent(in) :: filter
706 real(real64), intent(in) :: gmax
707
708 integer :: l, k, ispin
709
710 real(real64) :: alpha, beta_fs, rmax, rcut, gamma, beta_rs
711
712 push_sub(ps_filter)
713 call profiling_in("PS_FILTER")
714
715 select case (filter)
716 case (ps_filter_none)
717
718 case (ps_filter_ts)
719 alpha = 1.1_real64
720 gamma = m_two
721
722 if(.not. ps%no_vl) then
723 rmax = ps%vl%x_threshold
724 ! For the special case l=-1, we use l=0 when doing the filtering, but it is not clear this is correct
725 call spline_filter_mask(ps%vl, max(0, ps%llocal), rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
726 end if
727
728 do l = 0, ps%lmax
729 if (l == ps%llocal) cycle
730 do k = 1, ps%kbc
731 call spline_filter_mask(ps%kb(l, k), l, ps%rc_max, gmax, alpha, gamma, ps%projectors_sphere_threshold)
732 end do
733 end do
734
735 if (ps_has_nlcc(ps)) then
736 rmax = ps%core%x_threshold
737 call spline_filter_mask(ps%core, 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
738 end if
739
740 if (ps_has_density(ps)) then
741 do ispin = 1, ps%ispin
742 if (abs(spline_integral(ps%density(ispin))) > 1.0e-12_real64) then
743 rmax = ps%density(ispin)%x_threshold
744 call spline_filter_mask(ps%density(ispin), 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
745 call spline_force_pos(ps%density(ispin), ps%projectors_sphere_threshold)
746 end if
747
748 if (abs(spline_integral(ps%density_der(ispin))) > 1.0e-12_real64) then
749 rmax = ps%density_der(ispin)%x_threshold
750 call spline_filter_mask(ps%density_der(ispin), 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
751 end if
752 end do
753 end if
754
755 case (ps_filter_bsb)
756 alpha = 0.7_real64 ! The original was M_FOUR/7.0_real64
757 beta_fs = 18.0_real64
758 rcut = 2.5_real64
759 beta_rs = 0.4_real64
760
761 ! For the special case l=-1, we use l=0 when doing the filtering, but it is not clear this is correct
762 call spline_filter_bessel(ps%vl, max(0, ps%llocal), gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
763 do l = 0, ps%lmax
764 if (l == ps%llocal) cycle
765 do k = 1, ps%kbc
766 call spline_filter_bessel(ps%kb(l, k), l, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
767 end do
768 end do
769
770 if (ps_has_nlcc(ps)) then
771 call spline_filter_bessel(ps%core, 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
772 end if
773
774 if (ps_has_density(ps)) then
775 do ispin = 1, ps%ispin
776 call spline_filter_bessel(ps%density(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
777 call spline_force_pos(ps%density(ispin), ps%projectors_sphere_threshold)
778 call spline_filter_bessel(ps%density_der(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
779 end do
780 end if
782 end select
783
784 call profiling_out("PS_FILTER")
785 pop_sub(ps_filter)
786 end subroutine ps_filter
787
788 ! ---------------------------------------------------------
789 subroutine ps_check_bound(ps, eigen)
790 type(ps_t), intent(inout) :: ps
791 real(real64), intent(in) :: eigen(:,:)
792
793 integer :: i, is, ir
794 real(real64) :: ur1, ur2
795
796 push_sub(ps_check_bound)
797
798 ! Unbound states have positive eigenvalues
799 where(eigen > m_zero)
800 ps%bound = .false.
801 elsewhere
802 ps%bound = .true.
803 end where
804
805 ! We might not have information about the eigenvalues, so we need to check the wavefunctions
806 do i = 1, ps%conf%p
807 do is = 1, ps%ispin
808 if (.not. ps%bound(i, is)) cycle
809
810 do ir = ps%g%nrval, 3, -1
811 ! First we look for the outmost value that is not zero
812 if (abs(spline_eval(ps%ur(i, is), ps%g%rofi(ir))*ps%g%rofi(ir)) > m_zero) then
813 ! Usually bound states have exponentially decaying wavefunctions,
814 ! while unbound states have exponentially diverging
815 ! wavefunctions. Therefore we check if the wavefunctions
816 ! value is increasing with increasing radius. The fact
817 ! that we do not use the wavefunctions outmost value that
818 ! is not zero is on purpose, as some pseudopotential
819 ! generators do funny things with that point.
820 ur1 = spline_eval(ps%ur(i, is), ps%g%rofi(ir-2))*ps%g%rofi(ir-2)
821 ur2 = spline_eval(ps%ur(i, is), ps%g%rofi(ir-1))*ps%g%rofi(ir-1)
822 if ((ur1*ur2 > m_zero) .and. (abs(ur2) > abs(ur1))) ps%bound(i, is) = .false.
823 exit
824 end if
825 end do
826 end do
827 end do
828
829 pop_sub(ps_check_bound)
830 end subroutine ps_check_bound
831
832
833 ! ---------------------------------------------------------
834 subroutine ps_debug(ps, dir, namespace, gmax)
835 type(ps_t), intent(in) :: ps
836 character(len=*), intent(in) :: dir
837 type(namespace_t), intent(in) :: namespace
838 real(real64), intent(in) :: gmax
839
840 ! We will plot also some Fourier transforms.
841 type(spline_t), allocatable :: fw(:, :)
842
843 integer :: iunit
844 integer :: j, k, l
845
846 push_sub(ps_debug)
847
848 ! A text file with some basic data.
849 iunit = io_open(trim(dir)//'/pseudo-info', namespace, action='write')
850 write(iunit,'(a,/)') ps%label
851 write(iunit,'(a,a,/)') 'Format : ', ps_name(ps%file_format)
852 write(iunit,'(a,f6.3)') 'z : ', ps%z
853 write(iunit,'(a,f6.3,/)') 'zval : ', ps%z_val
854 write(iunit,'(a,i4)') 'lmax : ', ps%lmax
855 write(iunit,'(a,i4)') 'lloc : ', ps%llocal
856 write(iunit,'(a,i4,/)') 'kbc : ', ps%kbc
857 write(iunit,'(a,f9.5,/)') 'rcmax : ', ps%rc_max
858 write(iunit,'(a,/)') 'h matrix:'
859 do l = 0, ps%lmax
860 do k = 1, ps%kbc
861 write(iunit,'(10f9.5)') (ps%h(l, k, j), j = 1, ps%kbc)
862 end do
863 end do
864 if (allocated(ps%k)) then
865 write(iunit,'(/,a,/)') 'k matrix:'
866 do l = 0, ps%lmax
867 do k = 1, ps%kbc
868 write(iunit,'(10f9.5)') (ps%k(l, k, j), j = 1, ps%kbc)
869 end do
870 end do
871 end if
872
873 write(iunit,'(/,a)') 'orbitals:'
874 do j = 1, ps%conf%p
875 if (ps%ispin == 2) then
876 write(iunit,'(1x,a,i2,3x,a,i2,3x,a,f5.1,3x,a,l1,3x,a,f5.1,3x,a,f5.1)') 'n = ', ps%conf%n(j), 'l = ', ps%conf%l(j), &
877 'j = ', ps%conf%j(j), 'bound = ', all(ps%bound(j,:)), &
878 'occ(1) = ', ps%conf%occ(j, 1), 'occ(2) = ', ps%conf%occ(j, 2)
879 else
880 write(iunit,'(1x,a,i2,3x,a,i2,3x,a,f5.1,3x,a,l1,3x,a,f5.1)') 'n = ', ps%conf%n(j), 'l = ', ps%conf%l(j), &
881 'j = ', ps%conf%j(j), 'bound = ', all(ps%bound(j,:)), &
882 'occ = ', ps%conf%occ(j, 1)
883 end if
884 end do
885
886
887 call io_close(iunit)
888
889 ! Local part of the pseudopotential
890 iunit = io_open(trim(dir)//'/local', namespace, action='write')
891 call spline_print(ps%vl, iunit)
892 call io_close(iunit)
893
894 ! Local part of the pseudopotential
895 iunit = io_open(trim(dir)//'/local_long_range', namespace, action='write')
896 call spline_print(ps%vlr, iunit)
897 call io_close(iunit)
898
899 ! Local part of the pseudopotential
900 iunit = io_open(trim(dir)//'/local_long_range_density', namespace, action='write')
901 call spline_print(ps%nlr, iunit)
902 call io_close(iunit)
903
904 ! Fourier transform of the local part
905 iunit = io_open(trim(dir)//'/local_ft', namespace, action='write')
906 safe_allocate(fw(1, 1))
907 call spline_init(fw(1, 1))
908 call spline_3dft(ps%vl, fw(1, 1), ps%projectors_sphere_threshold, gmax = gmax)
909 call spline_print(fw(1, 1), iunit)
910 call spline_end(fw(1, 1))
911 safe_deallocate_a(fw)
912 call io_close(iunit)
913
914 ! Kleinman-Bylander projectors
915 iunit = io_open(trim(dir)//'/nonlocal', namespace, action='write')
916 call spline_print(ps%kb, iunit)
917 call io_close(iunit)
918
919 iunit = io_open(trim(dir)//'/nonlocal_derivative', namespace, action='write')
920 call spline_print(ps%dkb, iunit)
921 call io_close(iunit)
922
923 iunit = io_open(trim(dir)//'/nonlocal_ft', namespace, action='write')
924 safe_allocate(fw(0:ps%lmax, 1:ps%kbc))
925 call spline_init(fw)
926 do k = 0, ps%lmax
927 do j = 1, ps%kbc
928 call spline_besselft(ps%kb(k, j), fw(k, j), k, threshold=ps%projectors_sphere_threshold, &
929 gmax=m_four*gmax)
930 end do
931 end do
932 call spline_print(fw, iunit)
933 call spline_end(fw)
934 safe_deallocate_a(fw)
935 call io_close(iunit)
936
937 ! Pseudo-wavefunctions
938 iunit = io_open(trim(dir)//'/wavefunctions', namespace, action='write')
939 call spline_print(ps%ur, iunit)
940 call io_close(iunit)
941
942 ! Density
943 if (ps%has_density) then
944 iunit = io_open(trim(dir)//'/density', namespace, action='write')
945 call spline_print(ps%density, iunit)
946 call io_close(iunit)
947
948 iunit = io_open(trim(dir)//'/density_derivative', namespace, action='write')
949 call spline_print(ps%density_der, iunit)
950 call io_close(iunit)
951 end if
952
953 ! Non-linear core-corrections
954 if (ps_has_nlcc(ps)) then
955 iunit = io_open(trim(dir)//'/nlcc', namespace, action='write')
956 call spline_print(ps%core, iunit)
957 call io_close(iunit)
958 end if
959
960 pop_sub(ps_debug)
961 end subroutine ps_debug
962
963
964 ! ---------------------------------------------------------
965 subroutine ps_end(ps)
966 type(ps_t), intent(inout) :: ps
967
968 if (.not. allocated(ps%kb)) return
969
970 push_sub(ps_end)
971
972 if (ps%is_separated) then
973 call spline_end(ps%vlr)
974 call spline_end(ps%vlr_sq)
975 call spline_end(ps%nlr)
976 end if
977
978 call spline_end(ps%kb)
979 call spline_end(ps%dkb)
980 call spline_end(ps%ur)
981
982 call spline_end(ps%vl)
983 call spline_end(ps%core)
984 call spline_end(ps%core_der)
985
986 if (allocated(ps%density)) call spline_end(ps%density)
987 if (allocated(ps%density_der)) call spline_end(ps%density_der)
988
989 call logrid_end(ps%g)
990
991 safe_deallocate_a(ps%kb)
992 safe_deallocate_a(ps%dkb)
993 safe_deallocate_a(ps%ur)
994 safe_deallocate_a(ps%bound)
995 safe_deallocate_a(ps%h)
996 safe_deallocate_a(ps%k)
997 safe_deallocate_a(ps%density)
998 safe_deallocate_a(ps%density_der)
999
1000 pop_sub(ps_end)
1001 end subroutine ps_end
1002
1003
1004 ! ---------------------------------------------------------
1005 subroutine hgh_load(ps, ps_hgh)
1006 type(ps_t), intent(inout) :: ps
1007 type(ps_hgh_t), intent(inout) :: ps_hgh
1008
1009 push_sub(hgh_load)
1010
1011 ! Fixes some components of ps
1012 ps%nlcc = .false.
1013 if (ps%lmax >= 0) then
1014 ps%rc_max = 1.1_real64 * maxval(ps_hgh%kbr(0:ps%lmax)) ! Increase a little.
1015 else
1016 ps%rc_max = m_zero
1017 end if
1018 ps%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
1019 ps%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
1020
1021 ! now we fit the splines
1022 call get_splines()
1023
1024 pop_sub(hgh_load)
1025
1026 contains
1027
1028 ! ---------------------------------------------------------
1029 subroutine get_splines()
1030 integer :: l, is, j, ip
1031 real(real64), allocatable :: hato(:), dens(:)
1032
1033 push_sub(hgh_load.get_splines)
1034
1035 safe_allocate(hato(1:ps_hgh%g%nrval))
1036 safe_allocate(dens(1:ps_hgh%g%nrval))
1037
1038 ! Interpolate the KB-projection functions
1039 do l = 0, ps_hgh%l_max
1040 do j = 1, 3
1041 hato(:) = ps_hgh%kb(:, l, j)
1042 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%kb(l, j), ps%projectors_sphere_threshold)
1043 end do
1044 end do
1045
1046 ! Now the part corresponding to the local pseudopotential
1047 ! where the asymptotic part is subtracted
1048 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, ps_hgh%vlocal, ps%vl, ps%projectors_sphere_threshold)
1049
1050 ! Define the table for the pseudo-wavefunction components (using splines)
1051 ! with a correct normalization function
1052 do is = 1, ps%ispin
1053 dens = m_zero
1054 do l = 1, ps%conf%p
1055 hato(2:ps_hgh%g%nrval) = ps_hgh%rphi(2:ps_hgh%g%nrval, l)/ps_hgh%g%rofi(2:ps_hgh%g%nrval)
1056 hato(1) = hato(2)
1057
1058 do ip = 1, ps_hgh%g%nrval
1059 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(m_four*m_pi)
1060 end do
1061
1062 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
1063 end do
1064 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1065 end do
1066
1067 safe_deallocate_a(hato)
1068 safe_deallocate_a(dens)
1069
1070 pop_sub(hgh_load.get_splines)
1071 end subroutine get_splines
1072 end subroutine hgh_load
1073
1074
1075 ! ---------------------------------------------------------
1076 subroutine ps_grid_load(ps, ps_grid)
1077 type(ps_t), intent(inout) :: ps
1078 type(ps_in_grid_t), intent(in) :: ps_grid
1079
1080 push_sub(ps_grid_load)
1081
1082 ! Fixes some components of ps, read in ps_grid
1083 ps%z_val = ps_grid%zval
1084
1085 ps%nlcc = ps_grid%core_corrections
1086
1087 ps%h(0:ps%lmax, 1, 1) = ps_grid%dkbcos(1:ps%lmax+1)
1088
1089 ! Increasing radius a little, just in case.
1090 ! I have hard-coded a larger increase of the cutoff for the filtering.
1091 ps%rc_max = maxval(ps_grid%kb_radius(1:ps%lmax+1)) * 1.5_real64
1092
1093 ! now we fit the splines
1094 call get_splines(ps_grid%g)
1095
1096 ! Passes from Rydbergs to Hartrees.
1097 ps%h(0:ps%lmax,:,:) = ps%h(0:ps%lmax,:,:) / m_two
1098
1099 pop_sub(ps_grid_load)
1101 contains
1102
1103 subroutine get_splines(g)
1104 type(logrid_t), intent(in) :: g
1105
1106 real(real64), allocatable :: hato(:), dens(:)
1107 integer :: is, l, ir, nrc, ip
1108
1109 push_sub(ps_grid_load.get_splines)
1110
1111 safe_allocate(hato(1:g%nrval))
1112 safe_allocate(dens(1:g%nrval))
1113
1114 ! the wavefunctions
1115 do is = 1, ps%ispin
1116
1117 dens = m_zero
1118
1119 do l = 1, ps_grid%no_l_channels
1120 hato(2:) = ps_grid%rphi(2:, l, 1+is)/g%rofi(2:)
1121 hato(1) = first_point_extrapolate(g%rofi, hato)
1122
1123 if(ps%conf%occ(l, is) > m_epsilon) then
1124 do ip = 1, g%nrval
1125 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(m_four*m_pi)
1126 end do
1127 end if
1128
1129 call spline_fit(g%nrval, g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
1130
1131 end do
1132 call spline_fit(g%nrval, g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1133 end do
1134
1135
1136 ! the Kleinman-Bylander projectors
1137 do l = 1, ps%lmax+1
1138 nrc = logrid_index(g, ps_grid%kb_radius(l)) + 1
1139 hato(1:nrc) = ps_grid%KB(1:nrc, l)
1140 hato(nrc+1:g%nrval) = m_zero
1141
1142 call spline_fit(g%nrval, g%rofi, hato, ps%kb(l-1, 1), ps%projectors_sphere_threshold)
1143 end do
1144
1145 ! Now the part corresponding to the local pseudopotential
1146 ! where the asymptotic part is subtracted
1147 hato(:) = ps_grid%vlocal(:)/m_two
1148 call spline_fit(g%nrval, g%rofi, hato, ps%vl, ps%projectors_sphere_threshold)
1149
1150 if (ps_grid%core_corrections) then
1151 ! find cutoff radius
1152 hato(2:) = ps_grid%chcore(2:)/(m_four*m_pi*g%r2ofi(2:))
1153
1154 do ir = g%nrval-1, 2, -1
1155 if (hato(ir) > eps) then
1156 nrc = ir + 1
1157 exit
1158 end if
1159 end do
1160
1161 hato(nrc:g%nrval) = m_zero
1162 hato(1) = first_point_extrapolate(g%rofi, hato)
1163
1164 call spline_fit(g%nrval, g%rofi, hato, ps%core, ps%projectors_sphere_threshold)
1165 end if
1166
1167 safe_deallocate_a(hato)
1168 safe_deallocate_a(dens)
1169
1170 pop_sub(ps_grid_load.get_splines)
1171 end subroutine get_splines
1172 end subroutine ps_grid_load
1173
1174 ! ---------------------------------------------------------
1178 subroutine ps_xml_load(ps, ps_xml, namespace)
1179 type(ps_t), intent(inout) :: ps
1180 type(ps_xml_t), intent(in) :: ps_xml
1181 type(namespace_t), intent(in) :: namespace
1182
1183 integer :: ll, ip, is, ic, jc, ir, nrc, ii
1184 real(real64) :: rr, kbcos, kbnorm, dnrm, avgv, volume_element
1185 real(real64), allocatable :: vlocal(:), kbprojector(:), wavefunction(:), nlcc_density(:), dens(:,:)
1186 integer, allocatable :: cmap(:, :)
1187 real(real64), allocatable :: matrix(:, :), eigenvalues(:)
1188 logical :: density_is_known
1189 integer, allocatable :: order(:)
1190 real(real64), allocatable :: occ_tmp(:,:)
1191 real(real64), parameter :: tol_diagonal=1.0e-10_real64 ! tolerance for taking d_ij as a diagonal matrix
1192
1193 push_sub(ps_xml_load)
1194
1195 ! Vanderbilt-Kleinman-Bylander, i.e. more than one projector per l. Hamann uses 2 for ONCVPSP.
1196 if (pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
1197 if (ps%file_format == pseudo_format_psp8) then
1198 ps%relativistic_treatment = proj_j_average
1199 message(1) = "SOC from PSP8 is not currently supported."
1200 message(2) = "Only scalar relativistic effects will be considered."
1201 call messages_warning(2, namespace=namespace)
1202 else
1203 ps%relativistic_treatment = proj_j_dependent
1204 end if
1205 end if
1206
1207 ps%nlcc = ps_xml%nlcc
1208
1209 ps%z_val = ps_xml%valence_charge
1210
1211 ps%has_density = ps_xml%has_density
1212
1213 ! We start with the local potential
1214
1215 ! Outside of the grid, we use the analytical long-range tail
1216 ! TODO: check usefulness: ps%g%nrval is initialized as ps_xml%grid_size
1217 ! This makes no sense, ps_xml%grid can only go up to ps_xml%grid_size
1218 ! We should just fit ps_xml%potential directly
1219 ! Also, if ip > ps_xml%grid_size, most likely ps_xml%grid is not defined
1220 safe_allocate(vlocal(1:ps%g%nrval))
1221 do ip = 1, ps_xml%grid_size
1222 rr = ps_xml%grid(ip)
1223 if (ip <= ps_xml%grid_size) then
1224 vlocal(ip) = ps_xml%potential(ip, ps%llocal)
1225 else
1226 vlocal(ip) = -ps_xml%valence_charge/rr
1227 end if
1228 end do
1229
1230 ! We then fit this by a spline
1231 call spline_fit(ps%g%nrval, ps%g%rofi, vlocal, ps%vl, ps%projectors_sphere_threshold)
1232
1233 safe_deallocate_a(vlocal)
1234
1235 safe_allocate(kbprojector(1:ps%g%nrval))
1236 safe_allocate(wavefunction(1:ps%g%nrval))
1237
1238 kbprojector = m_zero
1239 wavefunction = m_zero
1240
1241 density_is_known = .false.
1242
1243 ! We then proceed with the nonlocal projectors and the orbitals
1244 if (ps_xml%kleinman_bylander) then
1245
1246 safe_allocate(cmap(0:ps_xml%lmax, 1:ps_xml%nchannels))
1247
1248 ! We order the different projectors and create some mappings
1249 ! the order of the channels is determined by spin orbit and the j value
1250 do ll = 0, ps_xml%lmax
1251 do ic = 1, ps_xml%nchannels
1252 cmap(ll, ic) = ic
1253
1254 if (ll == 0) cycle
1255 if (ll == ps_xml%llocal) cycle
1256 if (ps%relativistic_treatment /= proj_j_dependent) cycle
1257 if (ic > pseudo_nprojectors_per_l(ps_xml%pseudo, ll)) cycle ! This occurs for O and F for pseudodojo FR PBE for instance
1258
1259 ! This is Octopus convention:
1260 ! We treat l+1/2 first and l-1/2 after, so we order the projectors this way
1261
1262 assert(mod(ps_xml%nchannels, 2)==0)
1263 if (pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll - 1) then
1264 cmap(ll, ic) = int((ic-1)/2)*2 + 2
1265 else
1266 assert(pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll + 1)
1267 cmap(ll, ic) = int((ic-1)/2)*2 + 1
1268 end if
1269
1270 end do
1271
1272 ! check that all numbers are present for each l
1273 assert(sum(cmap(ll, 1:ps_xml%nchannels)) == (ps_xml%nchannels + 1)*ps_xml%nchannels/2)
1274 end do
1275
1276 assert(all(cmap >= 0 .and. cmap <= ps_xml%nchannels))
1277
1278
1279 ! Weight of the projectors
1280 ! This is a diagonal matrix after the following treatment
1281 ps%h = m_zero
1282
1283 ! We now take the matrix dij and we diagonalize it
1284 ! The nonlocal projectors are changed accordingly, and then fitted by splines
1285 if (pseudo_nprojectors(ps_xml%pseudo) > 0) then
1286 safe_allocate(matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels))
1287 safe_allocate(eigenvalues(1:ps_xml%nchannels))
1288
1289 do ll = 0, ps_xml%lmax
1290
1291 if (is_diagonal(ps_xml%nchannels, ps_xml%dij(ll, :, :), tol_diagonal) .or. &
1292 pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
1293 matrix = m_zero
1294 do ic = 1, ps_xml%nchannels
1295 eigenvalues(ic) = ps_xml%dij(ll, ic, ic)
1296 matrix(ic, ic) = m_one
1297 end do
1298 else
1299 ! diagonalize the coefficient matrix
1300 matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels) = ps_xml%dij(ll, 1:ps_xml%nchannels, 1:ps_xml%nchannels)
1301 call lalg_eigensolve(ps_xml%nchannels, matrix, eigenvalues)
1302 end if
1303
1304 do ic = 1, ps_xml%nchannels
1305 kbprojector = m_zero
1306 do jc = 1, ps_xml%nchannels
1307 call lalg_axpy(ps_xml%grid_size, matrix(jc, ic), ps_xml%projector(:, ll, jc), kbprojector)
1308 end do
1309
1310 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, cmap(ll, ic)), ps%projectors_sphere_threshold)
1311
1312 ps%h(ll, cmap(ll, ic), cmap(ll, ic)) = eigenvalues(ic)
1313
1314 end do
1315 end do
1316
1317 safe_deallocate_a(matrix)
1318 safe_deallocate_a(eigenvalues)
1319 end if
1320
1321 ps%conf%p = ps_xml%nwavefunctions
1322
1323 ! If we do not have a density but we have wavefunctions, we compute the
1324 ! pseudo-atomic density, from the pseudo-atomic wavefunctions
1325 ! In case we are doing a spin-polarized calculation, it is better to use the
1326 ! wavefunctions with spin-dependent occupations, instead the spin unresolved density
1327 ! given in the pseudopotential
1328 if ((.not. ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0) then
1329 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1330 dens = m_zero
1331 end if
1332
1333 do ii = 1, ps_xml%nwavefunctions
1334
1335 ps%conf%n(ii) = ps_xml%wf_n(ii)
1336 ps%conf%l(ii) = ps_xml%wf_l(ii)
1337
1338 if (ps%ispin == 2) then
1339 ps%conf%occ(ii, 1) = min(ps_xml%wf_occ(ii), m_two*ps_xml%wf_l(ii) + m_one)
1340 ps%conf%occ(ii, 2) = ps_xml%wf_occ(ii) - ps%conf%occ(ii, 1)
1341 else
1342 ps%conf%occ(ii, 1) = ps_xml%wf_occ(ii)
1343 end if
1344
1345 ps%conf%j(ii) = m_zero
1346 if (pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
1347 ps%conf%j(ii) = m_half*pseudo_wavefunction_2j(ps_xml%pseudo, ii)
1348 end if
1349
1350 do ip = 1, ps%g%nrval
1351 if (ip <= ps_xml%grid_size) then
1352 wavefunction(ip) = ps_xml%wavefunction(ip, ii)
1353 else
1354 wavefunction(ip) = m_zero
1355 end if
1356 end do
1357
1358 do is = 1, ps%ispin
1359 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ii, is), ps%projectors_sphere_threshold)
1360 end do
1361
1362 if (.not. ps_has_density(ps) .or. ps%ispin == 2) then
1363 do is = 1, ps%ispin
1364 do ip = 1, ps_xml%grid_size
1365 dens(ip, is) = dens(ip, is) + ps%conf%occ(ii, is)*wavefunction(ip)**2/(m_four*m_pi)
1366 end do
1367 end do
1368 end if
1369
1370 end do
1371
1372 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1373 if ((.not. ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0) then
1374 do is = 1, ps%ispin
1375 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1376 end do
1377 safe_deallocate_a(dens)
1378 density_is_known = .true.
1379 ps%has_density = .true.
1380 end if
1381
1382 safe_deallocate_a(cmap)
1383
1384 else !Not ps_xml%kleinman_bylander
1385
1386 !Get the occupations from the valence charge of the atom
1387 ps%conf%occ = m_zero
1388 ps%conf%symbol = ps%label(1:2)
1389 call ps_guess_atomic_occupations(namespace, ps%z, ps%z_val, ps%ispin, ps%conf)
1390
1391 ! In order to work in the following, we need to sort the occupations by angular momentum
1392 safe_allocate(order(1:ps_xml%lmax+1))
1393 safe_allocate(occ_tmp(1:ps_xml%lmax+1, 1:2))
1394 occ_tmp(:,:) = ps%conf%occ(1:ps_xml%lmax+1,1:2)
1395 call sort(ps%conf%l(1:ps_xml%lmax+1), order)
1396 do ll = 0, ps_xml%lmax
1397 ps%conf%occ(ll+1, 1:2) = occ_tmp(order(ll+1), 1:2)
1398 end do
1399 safe_deallocate_a(order)
1400 safe_deallocate_a(occ_tmp)
1401
1402 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1403 if (abs(sum(ps%conf%occ) - ps%z_val ) < m_epsilon) then
1404 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1405 dens = m_zero
1406 end if
1407
1408 do ll = 0, ps_xml%lmax
1409 ! we need to build the KB projectors
1410 ! the procedure was copied from ps_in_grid.F90 (r12967)
1411 dnrm = m_zero
1412 avgv = m_zero
1413 do ip = 1, ps_xml%grid_size
1414 rr = ps_xml%grid(ip)
1415 volume_element = rr**2*ps_xml%weights(ip)
1416 kbprojector(ip) = (ps_xml%potential(ip, ll) - ps_xml%potential(ip, ps%llocal))*ps_xml%wavefunction(ip, ll)
1417 dnrm = dnrm + kbprojector(ip)**2*volume_element
1418 avgv = avgv + kbprojector(ip)*ps_xml%wavefunction(ip, ll)*volume_element
1419 end do
1420
1421 kbcos = dnrm/(safe_tol(avgv,1.0e-20_real64))
1422 kbnorm = m_one/(safe_tol(sqrt(dnrm),1.0e-20_real64))
1423
1424 if (ll /= ps%llocal) then
1425 ps%h(ll, 1, 1) = kbcos
1426 kbprojector = kbprojector*kbnorm
1427 else
1428 ps%h(ll, 1, 1) = m_zero
1429 end if
1430
1431 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, 1), ps%projectors_sphere_threshold)
1432
1433 ! wavefunctions, for the moment we pad them with zero
1434 call lalg_copy(ps_xml%grid_size, ps_xml%wavefunction(:, ll), wavefunction)
1435 wavefunction(ps_xml%grid_size+1:ps%g%nrval) = m_zero
1436
1437 do is = 1, ps%ispin
1438 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ll + 1, is), ps%projectors_sphere_threshold)
1439 end do
1440
1441 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1442 if (abs(sum(ps%conf%occ) - ps%z_val) < m_epsilon) then
1443 do is = 1, ps%ispin
1444 do ip = 1, ps_xml%grid_size
1445 dens(ip, is) = dens(ip, is) + ps%conf%occ(ll+1, is)*wavefunction(ip)**2/(m_four*m_pi)
1446 end do
1447 end do
1448 end if
1449 end do
1450
1451 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1452 if (abs(sum(ps%conf%occ) - ps%z_val) < m_epsilon) then
1453 do is = 1, ps%ispin
1454 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1455 end do
1456 safe_deallocate_a(dens)
1457 density_is_known = .true.
1458 ps%has_density = .true.
1459 end if
1460 end if ! ps_xml%kleinman_bylander
1461
1462
1463 if (ps_has_density(ps) .and. .not. density_is_known) then
1464
1465 safe_allocate(dens(1:ps%g%nrval, 1))
1466
1467 dens(1:ps_xml%grid_size, 1) = ps_xml%density(1:ps_xml%grid_size)/ps%ispin
1468 dens(ps_xml%grid_size + 1:ps%g%nrval, 1) = m_zero
1469
1470 do is = 1, ps%ispin
1471 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:, 1), ps%density(is), ps%projectors_sphere_threshold)
1472 end do
1473
1474 safe_deallocate_a(dens)
1475 end if
1476
1477 ! Non-linear core-corrections
1478 ! We truncate the NLCC when below eps
1479 if (ps_xml%nlcc) then
1480
1481 safe_allocate(nlcc_density(1:ps%g%nrval))
1482
1483 call lalg_copy(ps_xml%grid_size, ps_xml%nlcc_density, nlcc_density)
1484
1485 ! find cutoff radius
1486 do ir = ps_xml%grid_size - 1, 1, -1
1487 if (nlcc_density(ir) > eps) then
1488 nrc = ir + 1
1489 exit
1490 end if
1491 end do
1492
1493 nlcc_density(nrc:ps%g%nrval) = m_zero
1494
1495 call spline_fit(ps%g%nrval, ps%g%rofi, nlcc_density, ps%core, ps%projectors_sphere_threshold)
1496
1497 safe_deallocate_a(nlcc_density)
1498 end if
1499
1500 call ps_getradius(ps)
1501
1502 !To be consistent with the other pseudopotentials, we are increasing the radius here
1503 ps%rc_max = ps%rc_max*1.05_real64
1504
1505 safe_deallocate_a(kbprojector)
1506 safe_deallocate_a(wavefunction)
1507
1508 pop_sub(ps_xml_load)
1509 end subroutine ps_xml_load
1510
1511
1512 ! ---------------------------------------------------------
1514 pure integer function ps_niwfs(ps)
1515 type(ps_t), intent(in) :: ps
1516
1517 integer :: i, l
1518
1519 ps_niwfs = 0
1520 do i = 1, ps%conf%p
1521 l = ps%conf%l(i)
1522 ps_niwfs = ps_niwfs + (2*l+1)
1523 end do
1524
1525 end function ps_niwfs
1526
1527 ! ---------------------------------------------------------
1529 pure integer function ps_bound_niwfs(ps)
1530 type(ps_t), intent(in) :: ps
1531
1532 integer :: i, l
1533
1534 ps_bound_niwfs = 0
1535 do i = 1, ps%conf%p
1536 l = ps%conf%l(i)
1537 if (any(.not. ps%bound(i,:))) cycle
1538 ps_bound_niwfs = ps_bound_niwfs + (2*l+1)
1539 end do
1540
1541 end function ps_bound_niwfs
1542
1543 !---------------------------------------
1544
1545 pure logical function ps_has_density(ps) result(has_density)
1546 type(ps_t), intent(in) :: ps
1547
1548 has_density = ps%has_density
1549
1550 end function ps_has_density
1551
1552 !---------------------------------------
1553
1554 pure logical function ps_has_nlcc(ps) result(has_nlcc)
1555 type(ps_t), intent(in) :: ps
1556
1557 has_nlcc = ps%nlcc
1558
1559 end function ps_has_nlcc
1560
1561 !---------------------------------------
1562 real(real64) function ps_density_volume(ps, namespace) result(volume)
1563 type(ps_t), intent(in) :: ps
1564 type(namespace_t), intent(in) :: namespace
1565
1566 integer :: ip, ispin
1567 real(real64) :: rr
1568 real(real64), allocatable ::vol(:)
1569 type(spline_t) :: volspl
1570
1571 push_sub(ps_density_volume)
1572
1573 if (.not. ps_has_density(ps)) then
1574 message(1) = "The pseudopotential does not contain an atomic density"
1575 call messages_fatal(1, namespace=namespace)
1576 end if
1577
1578 safe_allocate(vol(1:ps%g%nrval))
1579
1580 do ip = 1, ps%g%nrval
1581 rr = ps%g%rofi(ip)
1582 vol(ip) = m_zero
1583 do ispin = 1, ps%ispin
1584 vol(ip) = vol(ip) + spline_eval(ps%density(ispin), rr)*m_four*m_pi*rr**5
1585 end do
1586 end do
1587
1588 call spline_init(volspl)
1589 call spline_fit(ps%g%nrval, ps%g%rofi, vol, volspl, ps%projectors_sphere_threshold)
1590 volume = spline_integral(volspl)
1591 call spline_end(volspl)
1592
1593 safe_deallocate_a(vol)
1594
1595 pop_sub(ps_density_volume)
1596 end function ps_density_volume
1597
1601 subroutine ps_guess_atomic_occupations(namespace, zz, valcharge, ispin, conf)
1602 type(namespace_t), intent(in) :: namespace
1603 real(real64), intent(in) :: zz
1604 real(real64), intent(in) :: valcharge
1605 integer, intent(in) :: ispin
1606 type(valconf_t), intent(inout) :: conf
1607
1608 real(real64) :: val
1611
1612 val = valcharge
1613 conf%p = 0
1614 conf%l = 0
1615 conf%j = m_zero
1616 conf%n = 0
1617
1618 write(message(1), '(a,a,a)') 'Debug: Guessing the atomic occupations for ', trim(conf%symbol), "."
1619 call messages_info(1, namespace=namespace, debug_only=.true.)
1620
1621 assert(valcharge <= zz)
1622
1623 ! Here we populate the core states
1624 ! 1s state for all atoms after He
1625 if(int(zz) > 2 .and. val > zz - 2) then
1626 call fill_s_orbs(val, 2, 1)
1627 end if
1628 ! 2s state for all atoms after Be
1629 if(int(zz) > 4 .and. val > zz - 4) then
1630 call fill_s_orbs(val, 2, 2)
1631 end if
1632 ! 2p state for all atoms after Ne
1633 ! For pseudopotentials Al-Ar, we fill the 2s but not the 2p
1634 if(int(zz) > 18 .and. val > zz - 10) then
1635 call fill_p_orbs(val, 6, 2)
1636 end if
1637 ! 3s state for all atoms after Mg
1638 if(int(zz) > 12 .and. val > zz - 12) then
1639 call fill_s_orbs(val, 2, 3)
1640 end if
1641 ! 3p state for all atoms after Ar
1642 if(int(zz) > 18 .and. val > zz - 18) then
1643 call fill_p_orbs(val, 6, 3)
1644 end if
1645 ! 3d states for all atoms after Ni
1646 if(int(zz) > 28 .and. val > zz - 28) then
1647 call fill_d_orbs(val, 10, 3)
1648 end if
1649 ! 4s states for all atoms after Zn
1650 if(int(zz) > 30 .and. val > zz - 30) then
1651 call fill_s_orbs(val, 2, 4)
1652 end if
1653 ! 4p states for all atoms after Kr
1654 if(int(zz) > 36 .and. val > zz - 36) then
1655 call fill_p_orbs(val, 6, 4)
1656 end if
1657 ! 4d states for all atoms after Pd
1658 if(int(zz) > 46 .and. val > zz - 46) then
1659 call fill_d_orbs(val, 10, 4)
1660 end if
1661 ! 5s states for all atoms after Cd
1662 ! For Z=71 ot Z=80, the 4f is filled before the 5s/5p
1663 if((int(zz) > 48 .and. val > zz - 48) .or. &
1664 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 62)) then
1665 call fill_s_orbs(val, 2, 5)
1666 end if
1667 ! 5p states for all atoms after Xe
1668 ! For Z=71 ot Z=80, the 4f is filled before the 5s/5p
1669 if((int(zz) > 54 .and. val > zz - 54) .or. &
1670 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 68) ) then
1671 call fill_p_orbs(val, 6, 5)
1672 end if
1673 ! 4f states for all atoms after Yb
1674 ! Only after Z=80 for pseudopotentials
1675 if(int(zz) > 80 .and. val > zz - 68) then
1676 call fill_f_orbs(val, 14, 4)
1677 end if
1678
1679
1680 ! We leave here the valence states
1681 select case (int(zz))
1682 case (1)
1683 call fill_s_orbs(val, 1, 1) ! H 1s^1
1684 case (2)
1685 call fill_s_orbs(val, 2, 1) ! He 1s^2
1686 case (3)
1687 call fill_s_orbs(val, 1, 2) ! Li 1s^2 2s^1
1688 case (4)
1689 call fill_s_orbs(val, 2, 2) ! Be 1s^2 ; 2s^2
1690 case (5)
1691 call fill_p_orbs(val, 1, 2) ! B 1s^2 ; 2s^2 2p^1
1692 case (6)
1693 call fill_p_orbs(val, 2, 2) ! C 1s^2 ; 2s^2 2p^2
1694 case (7)
1695 call fill_p_orbs(val, 3, 2) ! N 1s^2 ; 2s^2 2
1696 case (8)
1697 call fill_p_orbs(val, 4, 2) ! O 1s^2 ; 2s^2 2p^4
1698 case (9)
1699 call fill_p_orbs(val, 5, 2) ! F 1s^2 ; 2s^2 2p^5
1700 case (10)
1701 call fill_p_orbs(val, 6, 2) ! Ne 1s^2 ; 2s^2 2p^6
1702 case (11)
1703 if(val > 6) call fill_p_orbs(val, 6, 2)
1704 call fill_s_orbs(val, 1, 3) ! Na 1s^2 2s^2 2p^6 ; 3s^1
1705 case (12)
1706 if(val > 6) call fill_p_orbs(val, 6, 2)
1707 call fill_s_orbs(val, 2, 3) ! Mg 1s^2 2s^2 2p^6 ; 3s^2
1708 case (13)
1709 if(val > 6) call fill_p_orbs(val, 6, 2)
1710 call fill_p_orbs(val, 1, 3) ! Al 1s^2 2s^2 2p^6 ; 3s^2 3p^1
1711 case (14)
1712 if(val > 6) call fill_p_orbs(val, 6, 2)
1713 call fill_p_orbs(val, 2, 3) ! Si 1s^2 2s^2 2p^6 ; 3s^2 3p^2
1714 case (15)
1715 if(val > 6) call fill_p_orbs(val, 6, 2)
1716 call fill_p_orbs(val, 3, 3) ! P 1s^2 2s^2 2p^6 ; 3s^2 3p^3
1717 case (16)
1718 if(val > 6) call fill_p_orbs(val, 6, 2)
1719 call fill_p_orbs(val, 4, 3) ! S 1s^2 2s^2 2p^6 ; 3s^2 3p^4
1720 case (17)
1721 if(val > 6) call fill_p_orbs(val, 6, 2)
1722 call fill_p_orbs(val, 5, 3) ! Cl 1s^2 2s^2 2p^6 ; 3s^2 3p^5
1723 case (18)
1724 if(val > 6) call fill_p_orbs(val, 6, 2)
1725 call fill_p_orbs(val, 6, 3) ! Ar 1s^2 2s^2 2p^6 ; 3s^2 3p^6
1726 case (19)
1727 call fill_s_orbs(val, 1, 4) ! K 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^1
1728 case (20)
1729 call fill_s_orbs(val, 2, 4) ! Ca 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2
1730 case (21)
1731 if (val > 1) call fill_s_orbs(val, 2, 4)
1732 call fill_d_orbs(val, 1, 3) ! Sc 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^1
1733 case (22)
1734 if (val > 2) call fill_s_orbs(val, 2, 4)
1735 call fill_d_orbs(val, 2, 3) ! Ti 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^2
1736 case (23)
1737 if (val > 3) call fill_s_orbs(val, 2, 4)
1738 call fill_d_orbs(val, 3, 3) ! V 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^3
1739 case (24)
1740 if (val > 4) call fill_s_orbs(val, 2, 4)
1741 call fill_d_orbs(val, 4, 3) ! Cr 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^4
1742 case (25)
1743 if (val > 5) call fill_s_orbs(val, 2, 4)
1744 call fill_d_orbs(val, 5, 3) ! Mn 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^5
1745 case (26)
1746 if (val > 6) call fill_s_orbs(val, 2, 4)
1747 call fill_d_orbs(val, 6, 3) ! Fe 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^6
1748 case (27)
1749 if (val > 7) call fill_s_orbs(val, 2, 4)
1750 call fill_d_orbs(val, 7, 3) ! Co 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^7
1751 case (28)
1752 if (val > 8 ) call fill_s_orbs(val, 2, 4)
1753 call fill_d_orbs(val, 8, 3) ! Ni 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^8
1754 case (29)
1755 call fill_s_orbs(val, 1, 4) ! Cu 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^1
1756 case (30)
1757 call fill_s_orbs(val, 2, 4) ! Zn 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2
1758 case (31)
1759 call fill_p_orbs(val, 1, 4) ! Ga 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^1
1760 case (32)
1761 call fill_p_orbs(val, 2, 4) ! Ge 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^2
1762 case (33)
1763 call fill_p_orbs(val, 3, 4) ! As 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^3
1764 case (34)
1765 call fill_p_orbs(val, 4, 4) ! Se 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^4
1766 case (35)
1767 call fill_p_orbs(val, 5, 4) ! Br 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^5
1768 case (36)
1769 call fill_p_orbs(val, 6, 4) ! Kr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^6
1770 case (37)
1771 call fill_s_orbs(val, 1, 5) ! Rb 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 5s^1
1772 case (38)
1773 call fill_s_orbs(val, 2, 5) ! Sr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 5s^2
1774 case (39)
1775 if (val > 2) call fill_d_orbs(val, 1, 4)
1776 call fill_s_orbs(val, 2, 5) ! Y 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^1 5s^2
1777 case (40)
1778 if (val > 2) call fill_d_orbs(val, 2, 4)
1779 call fill_s_orbs(val, 2, 5) ! Zr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^2 5s^2
1780 case (41)
1781 if (val > 1) call fill_d_orbs(val, 4, 4)
1782 call fill_s_orbs(val, 1, 5) ! Nb 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^4 5s^1
1783 case (42)
1784 if (val > 1) call fill_d_orbs(val, 5, 4)
1785 call fill_s_orbs(val, 1, 5) ! Mo 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^5 5s^1
1786 case (43)
1787 if (val > 2) call fill_d_orbs(val, 5, 4)
1788 call fill_s_orbs(val, 2, 5) ! Tc 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^5 5s^2
1789 case (44)
1790 if (val > 1) call fill_d_orbs(val, 7, 4)
1791 call fill_s_orbs(val, 1, 5) ! Ru 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^7 5s^1
1792 case (45)
1793 if (val > 1) call fill_d_orbs(val, 8, 4)
1794 call fill_s_orbs(val, 1, 5) ! Rh 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^8 5s^1
1795 case (46)
1796 call fill_d_orbs(val, 10, 4) ! Pd 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 ; 4d^10
1797 case (47)
1798 call fill_s_orbs(val, 1, 5) ! Ag 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^1
1799 case (48)
1800 call fill_s_orbs(val, 2, 5) ! Cd 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2
1801 case (49)
1802 call fill_p_orbs(val, 1, 5) ! In 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^1
1803 case (50)
1804 call fill_p_orbs(val, 2, 5) ! Sn 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^2
1805 case (51)
1806 call fill_p_orbs(val, 3, 5) ! Sb 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^3
1807 case (52)
1808 call fill_p_orbs(val, 4, 5) ! Te 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^4
1809 case (53)
1810 call fill_p_orbs(val, 5, 5) ! I 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^5
1811 case (54)
1812 call fill_p_orbs(val, 6, 5) ! Xe 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 ; 5s^2 5p^6
1813 case (55)
1814 call fill_s_orbs(val, 1, 6) ! Cs 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 6s^1
1815 case (56)
1816 call fill_s_orbs(val, 2, 6) ! Ba 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 6s^2
1817 case (57)
1818 if (val > 2) call fill_d_orbs(val, 1, 5)
1819 call fill_s_orbs(val, 2, 6) ! La 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 5d^1 6s^2
1820 case (58)
1821 if (val > 3) call fill_f_orbs(val, 1, 4)
1822 if (val > 2) call fill_d_orbs(val, 1, 5)
1823 call fill_s_orbs(val, 2, 6) ! Ce 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^1 5d^1 6s^2
1824 case (59)
1825 if (val > 2) call fill_f_orbs(val, 3, 4)
1826 call fill_s_orbs(val, 2, 6) ! Pr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^3 6s^2
1827 case (60)
1828 if (val > 2) call fill_f_orbs(val, 4, 4)
1829 call fill_s_orbs(val, 2, 6) ! Nd 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^4 6s^2
1830 case (61)
1831 if (val > 2) call fill_f_orbs(val, 5, 4)
1832 call fill_s_orbs(val, 2, 6) ! Pm 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^5 6s^2
1833 case (62)
1834 if (val > 2) call fill_f_orbs(val, 6, 4)
1835 call fill_s_orbs(val, 2, 6) ! Sm 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^6 6s^2
1836 case (63)
1837 if (val > 2) call fill_f_orbs(val, 7, 4)
1838 call fill_s_orbs(val, 2, 6) ! Eu 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^7 6s^2
1839 case (64)
1840 if (val > 3) call fill_f_orbs(val, 7, 4)
1841 if (val > 2) call fill_d_orbs(val, 1, 5)
1842 call fill_s_orbs(val, 2, 6) ! Gd 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^7 5d^1 6s^2
1843 case (65)
1844 if (val > 2) call fill_f_orbs(val, 9, 4)
1845 call fill_s_orbs(val, 2, 6) ! Tb 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^9 6s^2
1846 case (66)
1847 if (val > 2) call fill_f_orbs(val, 10, 4)
1848 call fill_s_orbs(val, 2, 6) ! Dy 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^10 6s^2
1849 case (67)
1850 if (val > 2) call fill_f_orbs(val, 11, 4)
1851 call fill_s_orbs(val, 2, 6) ! Ho 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^11 6s^2
1852 case (68)
1853 if (val > 2) call fill_f_orbs(val, 12, 4)
1854 call fill_s_orbs(val, 2, 6) ! Er 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^12 6s^2
1855 case (69)
1856 if (val > 2) call fill_f_orbs(val, 13, 4)
1857 call fill_s_orbs(val, 2, 6) ! Tm 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^13 6s^2
1858 case (70)
1859 if (val > 2) call fill_f_orbs(val, 14, 4)
1860 call fill_s_orbs(val, 2, 6) ! Yb 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 ; 4f^14 6s^2
1861 case (71)
1862 if (val > 3) call fill_f_orbs(val, 14, 4)
1863 if (val > 2) call fill_d_orbs(val, 1, 5)
1864 call fill_s_orbs(val, 2, 6) ! Lu 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^1 6s^2
1865 case (72)
1866 if (val > 4) call fill_f_orbs(val, 14, 4)
1867 if (val > 2) call fill_d_orbs(val, 2, 5)
1868 call fill_s_orbs(val, 2, 6) ! Hf 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^2 6s^2
1869 case (73)
1870 if (val > 5) call fill_f_orbs(val, 14, 4)
1871 if (val > 2) call fill_d_orbs(val, 3, 5)
1872 call fill_s_orbs(val, 2, 6) ! Ta 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^3 6s^2
1873 case (74)
1874 if (val > 6) call fill_f_orbs(val, 14, 4)
1875 if (val > 2) call fill_d_orbs(val, 4, 5)
1876 call fill_s_orbs(val, 2, 6) ! W 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^4 6s^2
1877 case (75)
1878 if (val > 7) call fill_f_orbs(val, 14, 4)
1879 if (val > 2) call fill_d_orbs(val, 5, 5)
1880 call fill_s_orbs(val, 2, 6) ! Re 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^5 6s^2
1881 case (76)
1882 if (val > 8) call fill_f_orbs(val, 14, 4)
1883 if (val > 2) call fill_d_orbs(val, 6, 5)
1884 call fill_s_orbs(val, 2, 6) ! Os 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^6 6s^2
1885 case (77)
1886 if (val > 9) call fill_f_orbs(val, 14, 4)
1887 if (val > 2) call fill_d_orbs(val, 7, 5)
1888 call fill_s_orbs(val, 2, 6) ! Ir 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^7 6s^2
1889 case (78)
1890 if (val > 10) call fill_f_orbs(val, 14, 4)
1891 if (val > 1) call fill_d_orbs(val, 9, 5)
1892 call fill_s_orbs(val, 1, 6) ! Pt 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 ; 5d^9 6s^1
1893 case (79)
1894 if (val > 21) call fill_f_orbs(val, 14, 4)
1895 if (val > 1) call fill_d_orbs(val, 10, 5)
1896 call fill_s_orbs(val, 1, 6) ! Au 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^1
1897 case (80)
1898 if (val > 12) call fill_f_orbs(val, 14, 4)
1899 if (val > 2) call fill_d_orbs(val, 10, 5)
1900 call fill_s_orbs(val, 2, 6) ! Hg 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2
1901 case (81)
1902 if (val > 10) call fill_d_orbs(val, 10, 5)
1903 if (val > 1) call fill_s_orbs(val, 2, 6)
1904 call fill_p_orbs(val, 1, 6) ! Tl 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^1
1905 case (82)
1906 if (val > 10) call fill_d_orbs(val, 10, 5)
1907 if (val > 2) call fill_s_orbs(val, 2, 6)
1908 call fill_p_orbs(val, 2, 6) ! Pb 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^2
1909 case (83)
1910 if (val > 10) call fill_d_orbs(val, 10, 5)
1911 if (val > 3) call fill_s_orbs(val, 2, 6)
1912 call fill_p_orbs(val, 3, 6) ! Bi 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^3
1913 case (84)
1914 if (val > 10) call fill_d_orbs(val, 10, 5)
1915 if (val > 4) call fill_s_orbs(val, 2, 6)
1916 call fill_p_orbs(val, 4, 6) ! Po 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^4
1917 case (85)
1918 if (val > 10) call fill_d_orbs(val, 10, 5)
1919 if (val > 5) call fill_s_orbs(val, 2, 6)
1920 call fill_p_orbs(val, 5, 6) ! At 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^5
1921 case (86)
1922 if (val > 10) call fill_d_orbs(val, 10, 5)
1923 if (val > 6) call fill_s_orbs(val, 2, 6)
1924 call fill_p_orbs(val, 6, 6) ! Rn 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 ; 6s^2 6p^6
1925 case (87)
1926 if (val > 1) call fill_p_orbs(val, 6, 6)
1927 call fill_s_orbs(val, 1, 7) ! Fr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 7s1
1928 case (88)
1929 if (val > 2) call fill_p_orbs(val, 6, 6)
1930 call fill_s_orbs(val, 2, 7) ! Ra 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 7s1
1931 case (89)
1932 if (val > 3) call fill_p_orbs(val, 6, 6)
1933 if (val > 2) call fill_d_orbs(val, 1, 6)
1934 call fill_s_orbs(val, 2, 7) ! Ac 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 6d1 7s2
1935 case (90)
1936 if (val > 4) call fill_p_orbs(val, 6, 6)
1937 if (val > 2) call fill_d_orbs(val, 2, 6)
1938 call fill_s_orbs(val, 2, 7) ! Th 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 6d2 7s2
1939 case (91)
1940 if (val > 5) call fill_p_orbs(val, 6, 6)
1941 if (val > 3) call fill_f_orbs(val, 2, 5)
1942 if (val > 2) call fill_d_orbs(val, 1, 6)
1943 call fill_s_orbs(val, 2, 7) ! Pa 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f2 6d1 7s2
1944 case (92)
1945 if (val > 6) call fill_p_orbs(val, 6, 6)
1946 if (val > 3) call fill_f_orbs(val, 3, 5)
1947 if (val > 2) call fill_d_orbs(val, 1, 6)
1948 call fill_s_orbs(val, 2, 7) ! U 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f3 6d1 7s2
1949 case (93)
1950 if (val > 7) call fill_p_orbs(val, 6, 6)
1951 if (val > 3) call fill_f_orbs(val, 4, 5)
1952 if (val > 2) call fill_d_orbs(val, 1, 6)
1953 call fill_s_orbs(val, 2, 7) ! Np 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f4 6d1 7s2
1954 case (94)
1955 if (val > 8) call fill_p_orbs(val, 6, 6)
1956 if (val > 2) call fill_f_orbs(val, 6, 5)
1957 call fill_s_orbs(val, 2, 7) ! Pu 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f6 7s2
1958 case (95)
1959 if (val > 9) call fill_p_orbs(val, 6, 6)
1960 if (val > 2) call fill_f_orbs(val, 7, 5)
1961 call fill_s_orbs(val, 2, 7) ! Am 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f7 7s2
1962 case (96)
1963 if (val > 10) call fill_p_orbs(val, 6, 6)
1964 if (val > 3) call fill_f_orbs(val, 7, 5)
1965 if (val > 2) call fill_d_orbs(val, 1, 6)
1966 call fill_s_orbs(val, 2, 7) ! Cm 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f7 6d1 7s2
1967 case (97)
1968 if (val > 11) call fill_p_orbs(val, 6, 6)
1969 if (val > 2) call fill_f_orbs(val, 9, 5)
1970 call fill_s_orbs(val, 2, 7) ! Bk 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f9 7s2
1971 case (98)
1972 if (val > 12) call fill_p_orbs(val, 6, 6)
1973 if (val > 2) call fill_f_orbs(val, 10, 5)
1974 call fill_s_orbs(val, 2, 7) ! Cf 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f10 7s2
1975 case (99)
1976 if (val > 13) call fill_p_orbs(val, 6, 6)
1977 if (val > 2) call fill_f_orbs(val, 11, 5)
1978 call fill_s_orbs(val, 2, 7) ! Es 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f11 7s2
1979 case (100)
1980 if (val > 14) call fill_p_orbs(val, 6, 6)
1981 if (val > 2) call fill_f_orbs(val, 12, 5)
1982 call fill_s_orbs(val, 2, 7) ! Fm 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f12 7s2
1983 case (101)
1984 if (val > 15) call fill_p_orbs(val, 6, 6)
1985 if (val > 2) call fill_f_orbs(val, 13, 5)
1986 call fill_s_orbs(val, 2, 7) ! Md 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f13 7s2
1987 case (102)
1988 if (val > 16) call fill_p_orbs(val, 6, 6)
1989 if (val > 2) call fill_f_orbs(val, 14, 5)
1990 call fill_s_orbs(val, 2, 7) ! No 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f14 7s2
1991 case (103)
1992 if (val > 3) call fill_f_orbs(val, 14, 5)
1993 if (val > 1) call fill_s_orbs(val, 2, 7)
1994 call fill_p_orbs(val, 1, 7) ! Lr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^14 5d^10 6s^2 6p^6 ; 5f14 7s2 7p1
1995 end select
1996
1997 !If we attributed all the electrons, everything went fine
1998 if (val < m_epsilon) then
1999 !In case of spin-polarized calculations, we properly distribute the electrons
2000 if (ispin == 2) then
2001 call valconf_unpolarized_to_polarized(conf)
2002 end if
2003 else
2004 conf%occ = m_zero
2005 message(1) = "Error in attributing atomic occupations"
2006 call messages_warning(1, namespace=namespace)
2007 end if
2008
2010
2011 contains
2012 subroutine fill_s_orbs(val, max_occ, nn)
2013 real(real64), intent(inout) :: val
2014 integer, intent(in) :: max_occ, nn
2015
2016 conf%p = conf%p + 1
2017 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2018 val = val - conf%occ(conf%p,1)
2019 conf%l(conf%p) = 0
2020 conf%n(conf%p) = nn
2021 end subroutine fill_s_orbs
2022
2023 subroutine fill_p_orbs(val, max_occ, nn)
2024 real(real64), intent(inout) :: val
2025 integer, intent(in) :: max_occ, nn
2026
2027 conf%p = conf%p + 1
2028 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2029 val = val - conf%occ(conf%p,1)
2030 conf%l(conf%p) = 1
2031 conf%n(conf%p) = nn
2032 end subroutine fill_p_orbs
2033
2034 subroutine fill_d_orbs(val, max_occ, nn)
2035 real(real64), intent(inout) :: val
2036 integer, intent(in) :: max_occ, nn
2037
2038 conf%p = conf%p + 1
2039 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2040 val = val - conf%occ(conf%p,1)
2041 conf%l(conf%p) = 2
2042 conf%n(conf%p) = nn
2043 end subroutine fill_d_orbs
2044
2045 subroutine fill_f_orbs(val, max_occ, nn)
2046 real(real64), intent(inout) :: val
2047 integer, intent(in) :: max_occ, nn
2048
2049 conf%p = conf%p + 1
2050 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2051 val = val - conf%occ(conf%p,1)
2052 conf%l(conf%p) = 3
2053 conf%n(conf%p) = nn
2054 end subroutine fill_f_orbs
2055
2056 end subroutine ps_guess_atomic_occupations
2057
2058#include "ps_pspio_inc.F90"
2059
2060end module ps_oct_m
2061
2062!! Local Variables:
2063!! mode: f90
2064!! coding: utf-8
2065!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:151
Some operations may be done for one spline-function, or for an array of them.
Definition: splines.F90:181
The integral may be done with or without integration limits, but we want the interface to be common.
Definition: splines.F90:175
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
subroutine, public valconf_unpolarized_to_polarized(conf)
Definition: atomic.F90:238
subroutine, public valconf_copy(cout, cin)
Definition: atomic.F90:168
real(real64), parameter, public m_two
Definition: global.F90:193
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
real(real64), parameter, public m_min_exp_arg
Definition: global.F90:210
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_three
Definition: global.F90:194
Definition: io.F90:116
subroutine, public io_close(iunit, grp)
Definition: io.F90:466
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:401
subroutine, public logrid_copy(grid_in, grid_out)
Definition: logrid.F90:230
integer function, public logrid_index(grid, rofi)
Definition: logrid.F90:258
subroutine, public logrid_end(grid)
Definition: logrid.F90:215
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
logical pure function, public is_diagonal(dim, matrix, tol)
Returns true is the matrix of size dim x dim is diagonal, given a relative tolerance.
Definition: math.F90:1511
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1023
subroutine, public messages_new_line()
Definition: messages.F90:1112
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
subroutine, public ps_cpi_end(ps_cpi)
Definition: ps_cpi.F90:185
subroutine, public ps_cpi_init(ps_cpi, filename, namespace)
Definition: ps_cpi.F90:150
subroutine, public ps_cpi_process(ps_cpi, lloc, namespace)
Definition: ps_cpi.F90:223
subroutine, public ps_fhi_end(ps_fhi)
Definition: ps_fhi.F90:185
subroutine, public ps_fhi_init(ps_fhi, filename, namespace)
Definition: ps_fhi.F90:151
subroutine, public ps_fhi_process(ps_fhi, lmax, lloc, namespace)
Definition: ps_fhi.F90:204
subroutine, public hgh_end(psp)
Definition: ps_hgh.F90:230
subroutine, public hgh_get_eigen(psp, eigen)
Definition: ps_hgh.F90:291
subroutine, public hgh_process(psp, namespace)
Definition: ps_hgh.F90:249
subroutine, public hgh_init(psp, filename, namespace)
Definition: ps_hgh.F90:185
real(real64) function, public first_point_extrapolate(x, y, high_order)
Definition: ps_in_grid.F90:439
Definition: ps.F90:116
pure logical function, public ps_has_density(ps)
Definition: ps.F90:1641
integer, parameter, public ps_filter_ts
Definition: ps.F90:166
subroutine, public ps_end(ps)
Definition: ps.F90:1061
subroutine hgh_load(ps, ps_hgh)
Definition: ps.F90:1101
subroutine, public ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
Definition: ps.F90:264
subroutine ps_info(ps, filename, namespace)
Definition: ps.F90:583
subroutine, public ps_separate(ps)
Separate the local potential into (soft) long-ranged and (hard) short-ranged parts.
Definition: ps.F90:694
subroutine ps_check_bound(ps, eigen)
Definition: ps.F90:885
subroutine ps_grid_load(ps, ps_grid)
Definition: ps.F90:1172
subroutine, public ps_debug(ps, dir, namespace, gmax)
Definition: ps.F90:930
integer, parameter, public proj_hgh
Definition: ps.F90:171
pure integer function, public ps_bound_niwfs(ps)
Returns the number of bound atomic orbitals taking into account then m quantum number multiplicity.
Definition: ps.F90:1625
subroutine, public ps_getradius(ps)
Definition: ps.F90:762
subroutine, public ps_guess_atomic_occupations(namespace, zz, valcharge, ispin, conf)
This routines provides, given Z and the number of valence electron the occupations of the orbitals....
Definition: ps.F90:1697
pure logical function, public ps_has_nlcc(ps)
Definition: ps.F90:1650
pure integer function, public ps_niwfs(ps)
Returns the number of atomic orbitals taking into account then m quantum number multiplicity.
Definition: ps.F90:1610
subroutine, public ps_pspio_init(ps, namespace, label, z, lmax, ispin, filename)
Definition: ps.F90:2174
integer, parameter, public proj_j_average
Fully-relativistic pseudopotentials with separate j-average and SOC terms.
Definition: ps.F90:177
integer, parameter, public proj_rkb
Definition: ps.F90:171
integer, parameter, public proj_j_dependent
Fully-relativistic j-dependent pseudopotentials.
Definition: ps.F90:177
subroutine, public ps_derivatives(ps)
Definition: ps.F90:782
subroutine, public ps_filter(ps, filter, gmax)
Definition: ps.F90:799
real(real64) function, public ps_density_volume(ps, namespace)
Definition: ps.F90:1658
integer, parameter, public ps_filter_bsb
Definition: ps.F90:166
integer, parameter, public proj_kb
Definition: ps.F90:171
subroutine ps_xml_load(ps, ps_xml, namespace)
Loads XML files for QSO, UPF1, UPF2, PSML, and PSP8 formats.
Definition: ps.F90:1274
subroutine, public ps_psf_init(pstm, ispin, filename, namespace)
Definition: ps_psf.F90:153
subroutine, public ps_psf_process(ps_psf, namespace, lmax, lloc)
Definition: ps_psf.F90:303
subroutine, public ps_psf_get_eigen(ps_psf, eigen)
Definition: ps_psf.F90:337
subroutine, public ps_psf_end(ps_psf)
Definition: ps_psf.F90:204
subroutine, public ps_xml_end(this)
Definition: ps_xml.F90:340
subroutine, public ps_xml_init(this, namespace, filename, fmt, ierr)
Definition: ps_xml.F90:166
integer, parameter, public pseudo_format_file_not_found
Definition: pseudo.F90:175
integer, parameter, public pseudo_type_kleinman_bylander
Definition: pseudo.F90:160
integer, parameter, public pseudo_format_qso
Definition: pseudo.F90:175
integer, parameter, public pseudo_format_hgh
Definition: pseudo.F90:175
integer, parameter, public pseudo_format_upf2
Definition: pseudo.F90:175
integer, parameter, public pseudo_type_paw
Definition: pseudo.F90:160
integer, parameter, public pseudo_format_cpi
Definition: pseudo.F90:175
integer, parameter, public pseudo_exchange_unknown
Definition: pseudo.F90:190
integer, parameter, public pseudo_type_semilocal
Definition: pseudo.F90:160
integer, parameter, public pseudo_type_ultrasoft
Definition: pseudo.F90:160
logical pure function, public pseudo_has_total_angular_momentum(pseudo)
Definition: pseudo.F90:552
integer, parameter, public pseudo_correlation_unknown
Definition: pseudo.F90:194
integer, parameter, public pseudo_format_psf
Definition: pseudo.F90:175
integer, parameter, public pseudo_format_fhi
Definition: pseudo.F90:175
integer, parameter, public pseudo_format_unknown
Definition: pseudo.F90:175
integer, parameter, public pseudo_format_psml
Definition: pseudo.F90:175
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:119
subroutine, public spline_filter_mask(spl, l, rmax, qmax, alpha, gamma, threshold)
subroutine, public spline_filter_bessel(spl, l, qmax, alpha, beta_fs, rcut, beta_rs, threshold)
subroutine, public spline_der(spl, dspl, threshold)
Definition: splines.F90:923
subroutine, public spline_force_pos(spl, threshold)
Definition: splines.F90:843
subroutine, public spline_fit(nrc, rofi, ffit, spl, threshold)
Definition: splines.F90:415
subroutine, public spline_3dft(spl, splw, threshold, gmax)
Definition: splines.F90:601
subroutine, public spline_besselft(spl, splw, l, threshold, gmax)
Definition: splines.F90:679
real(real64) function, public spline_eval(spl, x)
Definition: splines.F90:443
character(kind=c_char, len=1) function, dimension(:), allocatable, public string_f_to_c(f_string)
convert a Fortran string to a C string
Definition: string.F90:273
subroutine fill_s_orbs(val, max_occ, nn)
Definition: ps.F90:2108
subroutine fill_p_orbs(val, max_occ, nn)
Definition: ps.F90:2119
subroutine fill_d_orbs(val, max_occ, nn)
Definition: ps.F90:2130
subroutine get_splines()
Definition: ps.F90:1125
subroutine fill_f_orbs(val, max_occ, nn)
Definition: ps.F90:2141
remember that the FHI format is basically the CPI format with a header
Definition: ps_fhi.F90:139
The following data type contains: (a) the pseudopotential parameters, as read from a *....
Definition: ps_hgh.F90:146
A type storing the information and data about a pseudopotential.
Definition: ps.F90:188
the basic spline datatype
Definition: splines.F90:158
int true(void)