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 implicit none
50
51 private
52 public :: &
53 ps_t, &
54 ps_init, &
57 ps_filter, &
60 ps_debug, &
61 ps_niwfs, &
63 ps_end, &
68
69 integer, parameter, public :: &
70 PS_FILTER_NONE = 0, &
71 ps_filter_ts = 2, &
73
74 integer, public, parameter :: &
75 PROJ_NONE = 0, &
76 proj_hgh = 1, &
77 proj_kb = 2, &
78 proj_rkb = 3
79
80 integer, public, parameter :: &
81 PROJ_J_INDEPENDENT = 0, & !< Non-relativisitic or scalar-relativistic pseudopotentials
82 proj_j_dependent = 1, &
84
85 integer, parameter, public :: INVALID_L = 333
86
87 character(len=4), parameter :: ps_name(PSEUDO_FORMAT_UPF1:PSEUDO_FORMAT_PSP8) = &
88 (/"upf1", "upf2", "qso ", "psml", "psf ", "cpi ", "fhi ", "hgh ", "psp8"/)
89
91 type ps_t
92 ! Components are public by default
93 integer :: projector_type
94 integer :: relativistic_treatment
96
97 character(len=10), private :: label
98
99 integer, private :: ispin
100
101 real(real64), private :: z
102 real(real64) :: z_val
103 type(valconf_t) :: conf
104
105 type(logrid_t), private :: g
106
107 type(spline_t), allocatable :: ur(:, :)
108 type(spline_t), allocatable, private :: ur_sq(:, :)
109 logical, allocatable :: bound(:, :)
110
111 ! Kleinman-Bylander projectors stuff
112 integer :: lmax
113 integer :: llocal
115 type(spline_t) :: vl
116 logical :: no_vl = .false.
117
118 real(real64) :: projectors_sphere_threshold
125 real(real64) :: rc_max
126
127 integer :: kbc
128 integer :: projectors_per_l(5)
129 real(real64), allocatable :: h(:,:,:), k(:, :, :)
130 type(spline_t), allocatable :: kb(:, :)
131 type(spline_t), allocatable :: dkb(:, :)
132
133 logical :: nlcc
134 type(spline_t) :: core
135 type(spline_t) :: core_der
136
137
138 !LONG-RANGE PART OF THE LOCAL POTENTIAL
139
140 logical, private :: has_long_range
141 logical, private :: is_separated
142
143 type(spline_t), private :: vlr
144 type(spline_t) :: vlr_sq
146 type(spline_t) :: nlr
147
148 real(real64) :: sigma_erf
149
150 logical, private :: has_density
151 type(spline_t), allocatable :: density(:)
152 type(spline_t), allocatable :: density_der(:)
153
154 logical :: local
155 integer :: file_format
156 integer, private :: pseudo_type
157 integer :: exchange_functional
158 integer :: correlation_functional
159 end type ps_t
160
161 real(real64), parameter :: eps = 1.0e-8_real64
163contains
164
165
166 ! ---------------------------------------------------------
167 subroutine ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
168 type(ps_t), intent(out) :: ps
169 type(namespace_t), intent(in) :: namespace
170 character(len=10), intent(in) :: label
171 integer, intent(in) :: user_lmax
172 integer, intent(in) :: user_llocal
173 integer, intent(in) :: ispin
174 real(real64), intent(in) :: z
175 character(len=*), intent(in) :: filename
176
177 integer :: l, ii, ll, is, ierr
178 type(ps_psf_t) :: ps_psf
179 type(ps_cpi_t) :: ps_cpi
180 type(ps_fhi_t) :: ps_fhi
181 type(ps_hgh_t) :: ps_hgh
182 type(ps_xml_t) :: ps_xml
183 real(real64), allocatable :: eigen(:, :)
185 push_sub(ps_init)
187 ps%exchange_functional = pseudo_exchange_unknown
188 ps%correlation_functional = pseudo_correlation_unknown
189
190 ! Fix the threshold to calculate the radius of the projector-function localization spheres:
191
192 call messages_obsolete_variable(namespace, 'SpecieProjectorSphereThreshold', 'SpeciesProjectorSphereThreshold')
193
194 !%Variable SpeciesProjectorSphereThreshold
195 !%Type float
196 !%Default 0.001
197 !%Section System::Species
198 !%Description
199 !% The pseudopotentials may be composed of a local part, and a linear combination of nonlocal
200 !% operators. These nonlocal projectors have "projector" form, <math> \left| v \right> \left< v \right| </math>
201 !% (or, more generally speaking, <math> \left| u \right> \left< v \right| </math>).
202 !% These projectors are localized in real space -- that is, the function <math>v</math>
203 !% has a finite support around the nucleus. This region where the projectors are localized should
204 !% be small or else the computation time required to operate with them will be very large.
205 !%
206 !% In practice, this localization is fixed by requiring the definition of the projectors to be
207 !% contained in a sphere of a certain radius. This radius is computed by making sure that the
208 !% absolute value of the projector functions, at points outside the localization sphere, is
209 !% below a certain threshold. This threshold is set by <tt>SpeciesProjectorSphereThreshold</tt>.
210 !%End
211 call parse_variable(namespace, 'SpeciesProjectorSphereThreshold', 0.001_real64, ps%projectors_sphere_threshold)
212 if (ps%projectors_sphere_threshold <= m_zero) call messages_input_error(namespace, 'SpeciesProjectorSphereThreshold')
213
214 ps%file_format = pseudo_detect_format(filename)
215
216 if (ps%file_format == pseudo_format_file_not_found) then
217 call messages_write("Cannot open pseudopotential file '"//trim(filename)//"'.")
218 call messages_fatal(namespace=namespace)
219 end if
221 if (ps%file_format == pseudo_format_unknown) then
222 call messages_write("Cannot determine the pseudopotential type for species '"//trim(label)//"' from", &
223 new_line = .true.)
224 call messages_write("file '"//trim(filename)//"'.")
225 call messages_fatal(namespace=namespace)
226 end if
228 ps%label = label
229 ps%ispin = ispin
230 ps%relativistic_treatment = proj_j_independent
231 ps%projector_type = proj_kb
232 ps%sigma_erf = 0.625_real64 ! This is hard-coded to a reasonable value
234 ps%projectors_per_l = 0
235
236 select case (ps%file_format)
238 ps%has_density = .true.
239 case default
240 ps%has_density = .false.
241 end select
242
243 select case (ps%file_format)
245 ps%pseudo_type = pseudo_type_semilocal
246
247 call ps_psf_init(ps_psf, ispin, filename, namespace)
249 call valconf_copy(ps%conf, ps_psf%conf)
250 ps%z = z
251 ps%conf%z = nint(z) ! atomic number
252 ps%kbc = 1 ! only one projector per angular momentum
253
254 ps%lmax = ps_psf%ps_grid%no_l_channels - 1
255
256 if (user_lmax /= invalid_l) then
257 ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components.
258 if (user_lmax /= ps%lmax) then
259 message(1) = "lmax in Species block for " // trim(label) // &
260 " is larger than number available in pseudopotential."
261 call messages_fatal(1, namespace=namespace)
262 end if
263 end if
264
265 ps%conf%p = ps_psf%ps_grid%no_l_channels
266 if (ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0.
267
268 ! the local part of the pseudo
269 if (user_llocal == invalid_l) then
270 ps%llocal = 0
271 else
272 ps%llocal = user_llocal
273 end if
274
275 ps%projectors_per_l(1:ps%lmax+1) = 1
276 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
277
278 call ps_psf_process(ps_psf, namespace, ps%lmax, ps%llocal)
279 call logrid_copy(ps_psf%ps_grid%g, ps%g)
280
282 ps%pseudo_type = pseudo_type_semilocal
283
284 if (ps%file_format == pseudo_format_cpi) then
285 call ps_cpi_init(ps_cpi, trim(filename), namespace)
286 ps%conf%p = ps_cpi%ps_grid%no_l_channels
287 else
288 call ps_fhi_init(ps_fhi, trim(filename), namespace)
289 ps%conf%p = ps_fhi%ps_grid%no_l_channels
290 end if
291
292 ps%conf%z = nint(z)
293 ps%conf%symbol = label(1:2)
294 ps%conf%type = 1
295 do l = 1, ps%conf%p
296 ps%conf%l(l) = l - 1
297 end do
298
299 ps%z = z
300 ps%kbc = 1 ! only one projector per angular momentum
301
302 ps%lmax = ps%conf%p - 1
303
304 if (user_lmax /= invalid_l) then
305 ps%lmax = min(ps%lmax, user_lmax) ! Maybe the file does not have enough components.
306 if (user_lmax /= ps%lmax) then
307 message(1) = "lmax in Species block for " // trim(label) // &
308 " is larger than number available in pseudopotential."
309 call messages_fatal(1, namespace=namespace)
310 end if
311 end if
312
313 if (ps%lmax == 0) ps%llocal = 0 ! Vanderbilt is not acceptable if ps%lmax == 0.
314
315 ! the local part of the pseudo
316 if (user_llocal == invalid_l) then
317 ps%llocal = 0
318 else
319 ps%llocal = user_llocal
320 end if
321
322 ps%projectors_per_l(1:ps%lmax+1) = 1
323 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
324
325 if (ps%file_format == pseudo_format_cpi) then
326 call ps_cpi_process(ps_cpi, ps%llocal, namespace)
327 call logrid_copy(ps_cpi%ps_grid%g, ps%g)
328 else
329 call ps_fhi_process(ps_fhi, ps%lmax, ps%llocal, namespace)
330 call logrid_copy(ps_fhi%ps_grid%g, ps%g)
331 end if
332
333 case (pseudo_format_hgh)
334 ps%pseudo_type = pseudo_type_kleinman_bylander
335 ps%projector_type = proj_hgh
336
337 call hgh_init(ps_hgh, trim(filename), namespace)
338 call valconf_copy(ps%conf, ps_hgh%conf)
339
340 ps%z = z
341 ps%z_val = ps_hgh%z_val
342 ps%kbc = 3
343 ps%llocal = -1
344 ps%lmax = ps_hgh%l_max
345 ps%conf%symbol = label(1:2)
346 ps%sigma_erf = ps_hgh%rlocal ! We use the correct value
347
348 ps%projectors_per_l(1:ps%lmax+1) = 1
349
350 ! Get the occupations from the valence charge of the atom
351 ps%conf%occ = m_zero
352 ! We impose here non-spin-polarized occupations, to preserve the behavior of the code
353 ! We might want to change this to get a better LCAO guess
354 call ps_guess_atomic_occupations(namespace, ps%z, ps%z_val, 1, ps%conf)
355 ! We need the information to solve the Schrodinder equation
356 call valconf_copy(ps_hgh%conf, ps%conf)
357
358 call hgh_process(ps_hgh, namespace)
359 call logrid_copy(ps_hgh%g, ps%g)
360
361 ! In case of spin-polarized calculations, we properly distribute the electrons
362 if (ispin == 2) then
363 call valconf_unpolarized_to_polarized(ps%conf)
364 end if
365
366 case (pseudo_format_qso, pseudo_format_upf1, pseudo_format_upf2, pseudo_format_psml, pseudo_format_psp8)
367
368 call ps_xml_init(ps_xml, namespace, trim(filename), ps%file_format, ierr)
369
370 ps%pseudo_type = pseudo_type(ps_xml%pseudo)
371 ps%exchange_functional = pseudo_exchange(ps_xml%pseudo)
372 ps%correlation_functional = pseudo_correlation(ps_xml%pseudo)
373
374 ps%z = z
375 ps%conf%z = nint(z)
376
377 if (ps_xml%kleinman_bylander) then
378 ps%conf%p = ps_xml%nwavefunctions
379 else
380 ps%conf%p = ps_xml%lmax + 1
381 end if
382
383 do ll = 0, ps_xml%lmax
384 ps%conf%l(ll + 1) = ll
385 end do
386
387 ps%kbc = ps_xml%nchannels
388 ps%lmax = ps_xml%lmax
389
390 ps%projectors_per_l = 0
391 do ll = 0, ps_xml%lmax
392 ps%projectors_per_l(ll+1) = pseudo_nprojectors_per_l(ps_xml%pseudo, ll)
393 end do
394
395 if (ps_xml%kleinman_bylander) then
396 ps%llocal = ps_xml%llocal
397 else
398 ! we have several options
399 ps%llocal = 0 ! the default
400 if (ps_xml%llocal >= 0) ps%llocal = ps_xml%llocal ! the one given in the pseudopotential file
401 if (user_llocal /= invalid_l) ps%llocal = user_llocal ! user supplied local component
402 assert(ps%llocal >= 0)
403 assert(ps%llocal <= ps%lmax)
404 end if
405
406 ps%g%nrval = ps_xml%grid_size
407
408 safe_allocate(ps%g%rofi(1:ps%g%nrval))
409 safe_allocate(ps%g%r2ofi(1:ps%g%nrval))
410
411 do ii = 1, ps%g%nrval
412 ps%g%rofi(ii) = ps_xml%grid(ii)
413 ps%g%r2ofi(ii) = ps_xml%grid(ii)**2
414 end do
415
416 end select
417
418 ps%local = (ps%lmax == 0 .and. ps%llocal == 0 ) .or. (ps%lmax == -1 .and. ps%llocal == -1)
419
420 ! We allocate all the stuff
421 safe_allocate(ps%kb (0:ps%lmax, 1:ps%kbc))
422 safe_allocate(ps%dkb (0:ps%lmax, 1:ps%kbc))
423 safe_allocate(ps%ur (1:ps%conf%p, 1:ps%ispin))
424 safe_allocate(ps%ur_sq(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)
582
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)
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
781
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
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 call spline_end(ps%ur_sq)
982
983 call spline_end(ps%vl)
984 call spline_end(ps%core)
985 call spline_end(ps%core_der)
986
987 if (allocated(ps%density)) call spline_end(ps%density)
988 if (allocated(ps%density_der)) call spline_end(ps%density_der)
989
990 call logrid_end(ps%g)
991
992 safe_deallocate_a(ps%kb)
993 safe_deallocate_a(ps%dkb)
994 safe_deallocate_a(ps%ur)
995 safe_deallocate_a(ps%ur_sq)
996 safe_deallocate_a(ps%bound)
997 safe_deallocate_a(ps%h)
998 safe_deallocate_a(ps%k)
999 safe_deallocate_a(ps%density)
1000 safe_deallocate_a(ps%density_der)
1001
1002 pop_sub(ps_end)
1003 end subroutine ps_end
1004
1005
1006 ! ---------------------------------------------------------
1007 subroutine hgh_load(ps, ps_hgh)
1008 type(ps_t), intent(inout) :: ps
1009 type(ps_hgh_t), intent(inout) :: ps_hgh
1010
1011 push_sub(hgh_load)
1012
1013 ! Fixes some components of ps
1014 ps%nlcc = .false.
1015 if (ps%lmax >= 0) then
1016 ps%rc_max = 1.1_real64 * maxval(ps_hgh%kbr(0:ps%lmax)) ! Increase a little.
1017 else
1018 ps%rc_max = m_zero
1019 end if
1020 ps%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
1021 ps%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
1022
1023 ! now we fit the splines
1024 call get_splines()
1025
1026 pop_sub(hgh_load)
1027
1028 contains
1029
1030 ! ---------------------------------------------------------
1031 subroutine get_splines()
1032 integer :: l, is, j, ip
1033 real(real64), allocatable :: hato(:), dens(:)
1034
1035 push_sub(hgh_load.get_splines)
1036
1037 safe_allocate(hato(1:ps_hgh%g%nrval))
1038 safe_allocate(dens(1:ps_hgh%g%nrval))
1039
1040 ! Interpolate the KB-projection functions
1041 do l = 0, ps_hgh%l_max
1042 do j = 1, 3
1043 hato = ps_hgh%kb(:, l, j)
1044 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%kb(l, j), ps%projectors_sphere_threshold)
1045 end do
1046 end do
1047
1048 ! Now the part corresponding to the local pseudopotential
1049 ! where the asymptotic part is subtracted
1050 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, ps_hgh%vlocal, ps%vl, ps%projectors_sphere_threshold)
1051
1052 ! Define the table for the pseudo-wavefunction components (using splines)
1053 ! with a correct normalization function
1054 do is = 1, ps%ispin
1055 dens = m_zero
1056 do l = 1, ps%conf%p
1057 hato(2:ps_hgh%g%nrval) = ps_hgh%rphi(2:ps_hgh%g%nrval, l)/ps_hgh%g%rofi(2:ps_hgh%g%nrval)
1058 hato(1) = hato(2)
1059
1060 do ip = 1, ps_hgh%g%nrval
1061 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(m_four*m_pi)
1062 end do
1063
1064 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
1065 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%r2ofi, hato, ps%ur_sq(l, is), ps%projectors_sphere_threshold)
1066 end do
1067 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1068 end do
1069
1070 safe_deallocate_a(hato)
1071 safe_deallocate_a(dens)
1072
1073 pop_sub(hgh_load.get_splines)
1074 end subroutine get_splines
1075 end subroutine hgh_load
1076
1077
1078 ! ---------------------------------------------------------
1079 subroutine ps_grid_load(ps, ps_grid)
1080 type(ps_t), intent(inout) :: ps
1081 type(ps_in_grid_t), intent(in) :: ps_grid
1082
1083 push_sub(ps_grid_load)
1084
1085 ! Fixes some components of ps, read in ps_grid
1086 ps%z_val = ps_grid%zval
1087
1088 ps%nlcc = ps_grid%core_corrections
1089
1090 ps%h(0:ps%lmax, 1, 1) = ps_grid%dkbcos(1:ps%lmax+1)
1091
1092 ! Increasing radius a little, just in case.
1093 ! I have hard-coded a larger increase of the cutoff for the filtering.
1094 ps%rc_max = maxval(ps_grid%kb_radius(1:ps%lmax+1)) * 1.5_real64
1095
1096 ! now we fit the splines
1097 call get_splines(ps_grid%g)
1098
1099 ! Passes from Rydbergs to Hartrees.
1100 ps%h(0:ps%lmax,:,:) = ps%h(0:ps%lmax,:,:) / m_two
1101
1102 pop_sub(ps_grid_load)
1103
1104 contains
1105
1106 subroutine get_splines(g)
1107 type(logrid_t), intent(in) :: g
1108
1109 real(real64), allocatable :: hato(:), dens(:)
1110 integer :: is, l, ir, nrc, ip
1111
1112 push_sub(ps_grid_load.get_splines)
1113
1114 safe_allocate(hato(1:g%nrval))
1115 safe_allocate(dens(1:g%nrval))
1116
1117 ! the wavefunctions
1118 do is = 1, ps%ispin
1119
1120 dens = m_zero
1121
1122 do l = 1, ps_grid%no_l_channels
1123 hato(2:) = ps_grid%rphi(2:, l, 1+is)/g%rofi(2:)
1124 hato(1) = first_point_extrapolate(g%rofi, hato)
1125
1126 if(ps%conf%occ(l, is) > m_epsilon) then
1127 do ip = 1, g%nrval
1128 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(m_four*m_pi)
1129 end do
1130 end if
1131
1132 call spline_fit(g%nrval, g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
1133 call spline_fit(g%nrval, g%r2ofi, hato, ps%ur_sq(l, is), ps%projectors_sphere_threshold)
1134
1135 end do
1136 call spline_fit(g%nrval, g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1137 end do
1138
1139
1140 ! the Kleinman-Bylander projectors
1141 do l = 1, ps%lmax+1
1142 nrc = logrid_index(g, ps_grid%kb_radius(l)) + 1
1143 hato(1:nrc) = ps_grid%KB(1:nrc, l)
1144 hato(nrc+1:g%nrval) = m_zero
1145
1146 call spline_fit(g%nrval, g%rofi, hato, ps%kb(l-1, 1), ps%projectors_sphere_threshold)
1147 end do
1148
1149 ! Now the part corresponding to the local pseudopotential
1150 ! where the asymptotic part is subtracted
1151 hato(:) = ps_grid%vlocal(:)/m_two
1152 call spline_fit(g%nrval, g%rofi, hato, ps%vl, ps%projectors_sphere_threshold)
1153
1154 if (ps_grid%core_corrections) then
1155 ! find cutoff radius
1156 hato(2:) = ps_grid%chcore(2:)/(m_four*m_pi*g%r2ofi(2:))
1157
1158 do ir = g%nrval-1, 2, -1
1159 if (hato(ir) > eps) then
1160 nrc = ir + 1
1161 exit
1162 end if
1163 end do
1164
1165 hato(nrc:g%nrval) = m_zero
1166 hato(1) = first_point_extrapolate(g%rofi, hato)
1167
1168 call spline_fit(g%nrval, g%rofi, hato, ps%core, ps%projectors_sphere_threshold)
1169 end if
1170
1171 safe_deallocate_a(hato)
1172 safe_deallocate_a(dens)
1173
1174 pop_sub(ps_grid_load.get_splines)
1175 end subroutine get_splines
1176 end subroutine ps_grid_load
1177
1178 ! ---------------------------------------------------------
1182 subroutine ps_xml_load(ps, ps_xml, namespace)
1183 type(ps_t), intent(inout) :: ps
1184 type(ps_xml_t), intent(in) :: ps_xml
1185 type(namespace_t), intent(in) :: namespace
1186
1187 integer :: ll, ip, is, ic, jc, ir, nrc, ii
1188 real(real64) :: rr, kbcos, kbnorm, dnrm, avgv, volume_element
1189 real(real64), allocatable :: vlocal(:), kbprojector(:), wavefunction(:), nlcc_density(:), dens(:,:)
1190 integer, allocatable :: cmap(:, :)
1191 real(real64), allocatable :: matrix(:, :), eigenvalues(:)
1192 logical :: density_is_known
1193 integer, allocatable :: order(:)
1194 real(real64), allocatable :: occ_tmp(:,:)
1195 real(real64), parameter :: tol_diagonal=1.0e-10_real64 ! tolerance for taking d_ij as a diagonal matrix
1196
1197 push_sub(ps_xml_load)
1198
1199 ! Vanderbilt-Kleinman-Bylander, i.e. more than one projector per l. Hamann uses 2 for ONCVPSP.
1200 if (pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
1201 if (ps%file_format == pseudo_format_psp8) then
1202 ps%relativistic_treatment = proj_j_average
1203 message(1) = "SOC from PSP8 is not currently supported."
1204 message(2) = "Only scalar relativistic effects will be considered."
1205 call messages_warning(2, namespace=namespace)
1206 else
1207 ps%relativistic_treatment = proj_j_dependent
1208 end if
1209 end if
1210
1211 ps%nlcc = ps_xml%nlcc
1212
1213 ps%z_val = ps_xml%valence_charge
1214
1215 ps%has_density = ps_xml%has_density
1216
1217 ! We start with the local potential
1218
1219 ! Outside of the grid, we use the analytical long-range tail
1220 ! TODO: check usefulness: ps%g%nrval is initialized as ps_xml%grid_size
1221 ! This makes no sense, ps_xml%grid can only go up to ps_xml%grid_size
1222 ! We should just fit ps_xml%potential directly
1223 ! Also, if ip > ps_xml%grid_size, most likely ps_xml%grid is not defined
1224 safe_allocate(vlocal(1:ps%g%nrval))
1225 do ip = 1, ps_xml%grid_size
1226 rr = ps_xml%grid(ip)
1227 if (ip <= ps_xml%grid_size) then
1228 vlocal(ip) = ps_xml%potential(ip, ps%llocal)
1229 else
1230 vlocal(ip) = -ps_xml%valence_charge/rr
1231 end if
1232 end do
1233
1234 ! We then fit this by a spline
1235 call spline_fit(ps%g%nrval, ps%g%rofi, vlocal, ps%vl, ps%projectors_sphere_threshold)
1236
1237 safe_deallocate_a(vlocal)
1238
1239 safe_allocate(kbprojector(1:ps%g%nrval))
1240 safe_allocate(wavefunction(1:ps%g%nrval))
1241
1242 kbprojector = m_zero
1243 wavefunction = m_zero
1244
1245 density_is_known = .false.
1246
1247 ! We then proceed with the nonlocal projectors and the orbitals
1248 if (ps_xml%kleinman_bylander) then
1249
1250 safe_allocate(cmap(0:ps_xml%lmax, 1:ps_xml%nchannels))
1251
1252 ! We order the different projectors and create some mappings
1253 ! the order of the channels is determined by spin orbit and the j value
1254 do ll = 0, ps_xml%lmax
1255 do ic = 1, ps_xml%nchannels
1256 cmap(ll, ic) = ic
1257
1258 if (ll == 0) cycle
1259 if (ll == ps_xml%llocal) cycle
1260 if (ps%relativistic_treatment /= proj_j_dependent) cycle
1261 if (ic > pseudo_nprojectors_per_l(ps_xml%pseudo, ll)) cycle ! This occurs for O and F for pseudodojo FR PBE for instance
1262
1263 ! This is Octopus convention:
1264 ! We treat l+1/2 first and l-1/2 after, so we order the projectors this way
1265
1266 assert(mod(ps_xml%nchannels, 2)==0)
1267 if (pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll - 1) then
1268 cmap(ll, ic) = int((ic-1)/2)*2 + 2
1269 else
1270 assert(pseudo_projector_2j(ps_xml%pseudo, ll, ic) == 2*ll + 1)
1271 cmap(ll, ic) = int((ic-1)/2)*2 + 1
1272 end if
1273
1274 end do
1276 ! check that all numbers are present for each l
1277 assert(sum(cmap(ll, 1:ps_xml%nchannels)) == (ps_xml%nchannels + 1)*ps_xml%nchannels/2)
1278 end do
1279
1280 assert(all(cmap >= 0 .and. cmap <= ps_xml%nchannels))
1281
1282
1283 ! Weight of the projectors
1284 ! This is a diagonal matrix after the following treatment
1285 ps%h = m_zero
1286
1287 ! We now take the matrix dij and we diagonalize it
1288 ! The nonlocal projectors are changed accordingly, and then fitted by splines
1289 if (pseudo_nprojectors(ps_xml%pseudo) > 0) then
1290 safe_allocate(matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels))
1291 safe_allocate(eigenvalues(1:ps_xml%nchannels))
1292
1293 do ll = 0, ps_xml%lmax
1294
1295 if (is_diagonal(ps_xml%nchannels, ps_xml%dij(ll, :, :), tol_diagonal) .or. &
1296 pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
1297 matrix = m_zero
1298 do ic = 1, ps_xml%nchannels
1299 eigenvalues(ic) = ps_xml%dij(ll, ic, ic)
1300 matrix(ic, ic) = m_one
1301 end do
1302 else
1303 ! diagonalize the coefficient matrix
1304 matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels) = ps_xml%dij(ll, 1:ps_xml%nchannels, 1:ps_xml%nchannels)
1305 call lalg_eigensolve(ps_xml%nchannels, matrix, eigenvalues)
1306 end if
1307
1308 do ic = 1, ps_xml%nchannels
1309 kbprojector = m_zero
1310 do jc = 1, ps_xml%nchannels
1311 call lalg_axpy(ps_xml%grid_size, matrix(jc, ic), ps_xml%projector(:, ll, jc), kbprojector)
1312 end do
1313
1314 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, cmap(ll, ic)), ps%projectors_sphere_threshold)
1315
1316 ps%h(ll, cmap(ll, ic), cmap(ll, ic)) = eigenvalues(ic)
1317
1318 end do
1319 end do
1320
1321 safe_deallocate_a(matrix)
1322 safe_deallocate_a(eigenvalues)
1323 end if
1324
1325 ps%conf%p = ps_xml%nwavefunctions
1326
1327 ! If we do not have a density but we have wavefunctions, we compute the
1328 ! pseudo-atomic density, from the pseudo-atomic wavefunctions
1329 ! In case we are doing a spin-polarized calculation, it is better to use the
1330 ! wavefunctions with spin-dependent occupations, instead the spin unresolved density
1331 ! given in the pseudopotential
1332 if ((.not. ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0) then
1333 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1334 dens = m_zero
1335 end if
1336
1337 do ii = 1, ps_xml%nwavefunctions
1338
1339 ps%conf%n(ii) = ps_xml%wf_n(ii)
1340 ps%conf%l(ii) = ps_xml%wf_l(ii)
1341
1342 if (ps%ispin == 2) then
1343 ps%conf%occ(ii, 1) = min(ps_xml%wf_occ(ii), m_two*ps_xml%wf_l(ii) + m_one)
1344 ps%conf%occ(ii, 2) = ps_xml%wf_occ(ii) - ps%conf%occ(ii, 1)
1345 else
1346 ps%conf%occ(ii, 1) = ps_xml%wf_occ(ii)
1347 end if
1348
1349 ps%conf%j(ii) = m_zero
1350 if (pseudo_has_total_angular_momentum(ps_xml%pseudo)) then
1351 ps%conf%j(ii) = m_half*pseudo_wavefunction_2j(ps_xml%pseudo, ii)
1352 end if
1353
1354 do ip = 1, ps%g%nrval
1355 if (ip <= ps_xml%grid_size) then
1356 wavefunction(ip) = ps_xml%wavefunction(ip, ii)
1357 else
1358 wavefunction(ip) = m_zero
1359 end if
1360 end do
1361
1362 do is = 1, ps%ispin
1363 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ii, is), ps%projectors_sphere_threshold)
1364 call spline_fit(ps%g%nrval, ps%g%r2ofi, wavefunction, ps%ur_sq(ii, is), ps%projectors_sphere_threshold)
1365 end do
1366
1367 if (.not. ps_has_density(ps) .or. ps%ispin == 2) then
1368 do is = 1, ps%ispin
1369 do ip = 1, ps_xml%grid_size
1370 dens(ip, is) = dens(ip, is) + ps%conf%occ(ii, is)*wavefunction(ip)**2/(m_four*m_pi)
1371 end do
1372 end do
1373 end if
1374
1375 end do
1376
1377 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1378 if ((.not. ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0) then
1379 do is = 1, ps%ispin
1380 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1381 end do
1382 safe_deallocate_a(dens)
1383 density_is_known = .true.
1384 ps%has_density = .true.
1385 end if
1386
1387 safe_deallocate_a(cmap)
1388
1389 else !Not ps_xml%kleinman_bylander
1390
1391 !Get the occupations from the valence charge of the atom
1392 ps%conf%occ = m_zero
1393 ps%conf%symbol = ps%label(1:2)
1394 call ps_guess_atomic_occupations(namespace, ps%z, ps%z_val, ps%ispin, ps%conf)
1395
1396 ! In order to work in the following, we need to sort the occupations by angular momentum
1397 safe_allocate(order(1:ps_xml%lmax+1))
1398 safe_allocate(occ_tmp(1:ps_xml%lmax+1, 1:2))
1399 occ_tmp = ps%conf%occ
1400 call sort(ps%conf%l(1:ps_xml%lmax+1), order)
1401 do ll = 0, ps_xml%lmax
1402 ps%conf%occ(ll+1, 1:2) = occ_tmp(order(ll+1), 1:2)
1403 end do
1404 safe_deallocate_a(order)
1405 safe_deallocate_a(occ_tmp)
1406
1407 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1408 if (abs(sum(ps%conf%occ) - ps%z_val ) < m_epsilon) then
1409 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1410 dens = m_zero
1411 end if
1412
1413 do ll = 0, ps_xml%lmax
1414 ! we need to build the KB projectors
1415 ! the procedure was copied from ps_in_grid.F90 (r12967)
1416 dnrm = m_zero
1417 avgv = m_zero
1418 do ip = 1, ps_xml%grid_size
1419 rr = ps_xml%grid(ip)
1420 volume_element = rr**2*ps_xml%weights(ip)
1421 kbprojector(ip) = (ps_xml%potential(ip, ll) - ps_xml%potential(ip, ps%llocal))*ps_xml%wavefunction(ip, ll)
1422 dnrm = dnrm + kbprojector(ip)**2*volume_element
1423 avgv = avgv + kbprojector(ip)*ps_xml%wavefunction(ip, ll)*volume_element
1424 end do
1425
1426 kbcos = dnrm/(safe_tol(avgv,1.0e-20_real64))
1427 kbnorm = m_one/(safe_tol(sqrt(dnrm),1.0e-20_real64))
1428
1429 if (ll /= ps%llocal) then
1430 ps%h(ll, 1, 1) = kbcos
1431 kbprojector = kbprojector*kbnorm
1432 else
1433 ps%h(ll, 1, 1) = m_zero
1434 end if
1435
1436 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, 1), ps%projectors_sphere_threshold)
1437
1438 ! wavefunctions, for the moment we pad them with zero
1439 call lalg_copy(ps_xml%grid_size, ps_xml%wavefunction(:, ll), wavefunction)
1440 wavefunction(ps_xml%grid_size+1:ps%g%nrval) = m_zero
1441
1442 do is = 1, ps%ispin
1443 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ll + 1, is), ps%projectors_sphere_threshold)
1444 call spline_fit(ps%g%nrval, ps%g%r2ofi, wavefunction, ps%ur_sq(ll + 1, is), ps%projectors_sphere_threshold)
1445 end do
1446
1447 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1448 if (abs(sum(ps%conf%occ) - ps%z_val) < m_epsilon) then
1449 do is = 1, ps%ispin
1450 do ip = 1, ps_xml%grid_size
1451 dens(ip, is) = dens(ip, is) + ps%conf%occ(ll+1, is)*wavefunction(ip)**2/(m_four*m_pi)
1452 end do
1453 end do
1454 end if
1455 end do
1456
1457 !If we assigned all the valence electrons, we can compute the (spin-resolved) atomic density
1458 if (abs(sum(ps%conf%occ) - ps%z_val) < m_epsilon) then
1459 do is = 1, ps%ispin
1460 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1461 end do
1462 safe_deallocate_a(dens)
1463 density_is_known = .true.
1464 ps%has_density = .true.
1465 end if
1466 end if ! ps_xml%kleinman_bylander
1467
1468
1469 if (ps_has_density(ps) .and. .not. density_is_known) then
1470
1471 safe_allocate(dens(1:ps%g%nrval, 1))
1472
1473 dens(1:ps_xml%grid_size, 1) = ps_xml%density(1:ps_xml%grid_size)/ps%ispin
1474 dens(ps_xml%grid_size + 1:ps%g%nrval, 1) = m_zero
1475
1476 do is = 1, ps%ispin
1477 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:, 1), ps%density(is), ps%projectors_sphere_threshold)
1478 end do
1479
1480 safe_deallocate_a(dens)
1481 end if
1482
1483 ! Non-linear core-corrections
1484 ! We truncate the NLCC when below eps
1485 if (ps_xml%nlcc) then
1486
1487 safe_allocate(nlcc_density(1:ps%g%nrval))
1488
1489 call lalg_copy(ps_xml%grid_size, ps_xml%nlcc_density, nlcc_density)
1490
1491 ! find cutoff radius
1492 do ir = ps_xml%grid_size - 1, 1, -1
1493 if (nlcc_density(ir) > eps) then
1494 nrc = ir + 1
1495 exit
1496 end if
1497 end do
1498
1499 nlcc_density(nrc:ps%g%nrval) = m_zero
1500
1501 call spline_fit(ps%g%nrval, ps%g%rofi, nlcc_density, ps%core, ps%projectors_sphere_threshold)
1502
1503 safe_deallocate_a(nlcc_density)
1504 end if
1505
1506 call ps_getradius(ps)
1507
1508 !To be consistent with the other pseudopotentials, we are increasing the radius here
1509 ps%rc_max = ps%rc_max*1.05_real64
1510
1511 safe_deallocate_a(kbprojector)
1512 safe_deallocate_a(wavefunction)
1513
1514 pop_sub(ps_xml_load)
1515 end subroutine ps_xml_load
1516
1517
1518 ! ---------------------------------------------------------
1520 pure integer function ps_niwfs(ps)
1521 type(ps_t), intent(in) :: ps
1522
1523 integer :: i, l
1524
1525 ps_niwfs = 0
1526 do i = 1, ps%conf%p
1527 l = ps%conf%l(i)
1528 ps_niwfs = ps_niwfs + (2*l+1)
1529 end do
1530
1531 end function ps_niwfs
1532
1533 ! ---------------------------------------------------------
1535 pure integer function ps_bound_niwfs(ps)
1536 type(ps_t), intent(in) :: ps
1537
1538 integer :: i, l
1539
1540 ps_bound_niwfs = 0
1541 do i = 1, ps%conf%p
1542 l = ps%conf%l(i)
1543 if (any(.not. ps%bound(i,:))) cycle
1544 ps_bound_niwfs = ps_bound_niwfs + (2*l+1)
1545 end do
1546
1547 end function ps_bound_niwfs
1548
1549 !---------------------------------------
1550
1551 pure logical function ps_has_density(ps) result(has_density)
1552 type(ps_t), intent(in) :: ps
1553
1554 has_density = ps%has_density
1555
1556 end function ps_has_density
1557
1558 !---------------------------------------
1559
1560 pure logical function ps_has_nlcc(ps) result(has_nlcc)
1561 type(ps_t), intent(in) :: ps
1562
1563 has_nlcc = ps%nlcc
1564
1565 end function ps_has_nlcc
1566
1567 !---------------------------------------
1568 real(real64) function ps_density_volume(ps, namespace) result(volume)
1569 type(ps_t), intent(in) :: ps
1570 type(namespace_t), intent(in) :: namespace
1571
1572 integer :: ip, ispin
1573 real(real64) :: rr
1574 real(real64), allocatable ::vol(:)
1575 type(spline_t) :: volspl
1576
1577 push_sub(ps_density_volume)
1578
1579 if (.not. ps_has_density(ps)) then
1580 message(1) = "The pseudopotential does not contain an atomic density"
1581 call messages_fatal(1, namespace=namespace)
1582 end if
1583
1584 safe_allocate(vol(1:ps%g%nrval))
1585
1586 do ip = 1, ps%g%nrval
1587 rr = ps%g%rofi(ip)
1588 vol(ip) = m_zero
1589 do ispin = 1, ps%ispin
1590 vol(ip) = vol(ip) + spline_eval(ps%density(ispin), rr)*m_four*m_pi*rr**5
1591 end do
1592 end do
1593
1594 call spline_init(volspl)
1595 call spline_fit(ps%g%nrval, ps%g%rofi, vol, volspl, ps%projectors_sphere_threshold)
1596 volume = spline_integral(volspl)
1597 call spline_end(volspl)
1598
1599 safe_deallocate_a(vol)
1600
1601 pop_sub(ps_density_volume)
1602 end function ps_density_volume
1603
1607 subroutine ps_guess_atomic_occupations(namespace, zz, valcharge, ispin, conf)
1608 type(namespace_t), intent(in) :: namespace
1609 real(real64), intent(in) :: zz
1610 real(real64), intent(in) :: valcharge
1611 integer, intent(in) :: ispin
1612 type(valconf_t), intent(inout) :: conf
1614 real(real64) :: val
1615
1617
1618 val = valcharge
1619 conf%p = 0
1620 conf%l = 0
1621 conf%j = m_zero
1622 conf%n = 0
1623
1624 write(message(1), '(a,a,a)') 'Debug: Guessing the atomic occupations for ', trim(conf%symbol), "."
1625 call messages_info(1, namespace=namespace, debug_only=.true.)
1626
1627 assert(valcharge <= zz)
1629 ! Here we populate the core states
1630 ! 1s state for all atoms after He
1631 if(int(zz) > 2 .and. val > zz - 2) then
1632 call fill_s_orbs(val, 2, 1)
1633 end if
1634 ! 2s state for all atoms after Be
1635 if(int(zz) > 4 .and. val > zz - 4) then
1636 call fill_s_orbs(val, 2, 2)
1637 end if
1638 ! 2p state for all atoms after Ne
1639 ! For pseudopotentials Al-Ar, we fill the 2s but not the 2p
1640 if(int(zz) > 18 .and. val > zz - 10) then
1641 call fill_p_orbs(val, 6, 2)
1642 end if
1643 ! 3s state for all atoms after Mg
1644 if(int(zz) > 12 .and. val > zz - 12) then
1645 call fill_s_orbs(val, 2, 3)
1646 end if
1647 ! 3p state for all atoms after Ar
1648 if(int(zz) > 18 .and. val > zz - 18) then
1649 call fill_p_orbs(val, 6, 3)
1650 end if
1651 ! 3d states for all atoms after Ni
1652 if(int(zz) > 28 .and. val > zz - 28) then
1653 call fill_d_orbs(val, 10, 3)
1654 end if
1655 ! 4s states for all atoms after Zn
1656 if(int(zz) > 30 .and. val > zz - 30) then
1657 call fill_s_orbs(val, 2, 4)
1658 end if
1659 ! 4p states for all atoms after Kr
1660 if(int(zz) > 36 .and. val > zz - 36) then
1661 call fill_p_orbs(val, 6, 4)
1662 end if
1663 ! 4d states for all atoms after Pd
1664 if(int(zz) > 46 .and. val > zz - 46) then
1665 call fill_d_orbs(val, 10, 4)
1666 end if
1667 ! 5s states for all atoms after Cd
1668 ! For Z=71 ot Z=80, the 4f is filled before the 5s/5p
1669 if((int(zz) > 48 .and. val > zz - 48) .or. &
1670 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 62)) then
1671 call fill_s_orbs(val, 2, 5)
1672 end if
1673 ! 5p states for all atoms after Xe
1674 ! For Z=71 ot Z=80, the 4f is filled before the 5s/5p
1675 if((int(zz) > 54 .and. val > zz - 54) .or. &
1676 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 68) ) then
1677 call fill_p_orbs(val, 6, 5)
1678 end if
1679 ! 4f states for all atoms after Yb
1680 ! Only after Z=80 for pseudopotentials
1681 if(int(zz) > 80 .and. val > zz - 68) then
1682 call fill_f_orbs(val, 14, 4)
1683 end if
1684
1685
1686 ! We leave here the valence states
1687 select case (int(zz))
1688 case (1)
1689 call fill_s_orbs(val, 1, 1) ! H 1s^1
1690 case (2)
1691 call fill_s_orbs(val, 2, 1) ! He 1s^2
1692 case (3)
1693 call fill_s_orbs(val, 1, 2) ! Li 1s^2 2s^1
1694 case (4)
1695 call fill_s_orbs(val, 2, 2) ! Be 1s^2 ; 2s^2
1696 case (5)
1697 call fill_p_orbs(val, 1, 2) ! B 1s^2 ; 2s^2 2p^1
1698 case (6)
1699 call fill_p_orbs(val, 2, 2) ! C 1s^2 ; 2s^2 2p^2
1700 case (7)
1701 call fill_p_orbs(val, 3, 2) ! N 1s^2 ; 2s^2 2
1702 case (8)
1703 call fill_p_orbs(val, 4, 2) ! O 1s^2 ; 2s^2 2p^4
1704 case (9)
1705 call fill_p_orbs(val, 5, 2) ! F 1s^2 ; 2s^2 2p^5
1706 case (10)
1707 call fill_p_orbs(val, 6, 2) ! Ne 1s^2 ; 2s^2 2p^6
1708 case (11)
1709 if(val > 6) call fill_p_orbs(val, 6, 2)
1710 call fill_s_orbs(val, 1, 3) ! Na 1s^2 2s^2 2p^6 ; 3s^1
1711 case (12)
1712 if(val > 6) call fill_p_orbs(val, 6, 2)
1713 call fill_s_orbs(val, 2, 3) ! Mg 1s^2 2s^2 2p^6 ; 3s^2
1714 case (13)
1715 if(val > 6) call fill_p_orbs(val, 6, 2)
1716 call fill_p_orbs(val, 1, 3) ! Al 1s^2 2s^2 2p^6 ; 3s^2 3p^1
1717 case (14)
1718 if(val > 6) call fill_p_orbs(val, 6, 2)
1719 call fill_p_orbs(val, 2, 3) ! Si 1s^2 2s^2 2p^6 ; 3s^2 3p^2
1720 case (15)
1721 if(val > 6) call fill_p_orbs(val, 6, 2)
1722 call fill_p_orbs(val, 3, 3) ! P 1s^2 2s^2 2p^6 ; 3s^2 3p^3
1723 case (16)
1724 if(val > 6) call fill_p_orbs(val, 6, 2)
1725 call fill_p_orbs(val, 4, 3) ! S 1s^2 2s^2 2p^6 ; 3s^2 3p^4
1726 case (17)
1727 if(val > 6) call fill_p_orbs(val, 6, 2)
1728 call fill_p_orbs(val, 5, 3) ! Cl 1s^2 2s^2 2p^6 ; 3s^2 3p^5
1729 case (18)
1730 if(val > 6) call fill_p_orbs(val, 6, 2)
1731 call fill_p_orbs(val, 6, 3) ! Ar 1s^2 2s^2 2p^6 ; 3s^2 3p^6
1732 case (19)
1733 call fill_s_orbs(val, 1, 4) ! K 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^1
1734 case (20)
1735 call fill_s_orbs(val, 2, 4) ! Ca 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2
1736 case (21)
1737 if (val > 1) call fill_s_orbs(val, 2, 4)
1738 call fill_d_orbs(val, 1, 3) ! Sc 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^1
1739 case (22)
1740 if (val > 2) call fill_s_orbs(val, 2, 4)
1741 call fill_d_orbs(val, 2, 3) ! Ti 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^2
1742 case (23)
1743 if (val > 3) call fill_s_orbs(val, 2, 4)
1744 call fill_d_orbs(val, 3, 3) ! V 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^3
1745 case (24)
1746 if (val > 4) call fill_s_orbs(val, 2, 4)
1747 call fill_d_orbs(val, 4, 3) ! Cr 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^4
1748 case (25)
1749 if (val > 5) call fill_s_orbs(val, 2, 4)
1750 call fill_d_orbs(val, 5, 3) ! Mn 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^5
1751 case (26)
1752 if (val > 6) call fill_s_orbs(val, 2, 4)
1753 call fill_d_orbs(val, 6, 3) ! Fe 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^6
1754 case (27)
1755 if (val > 7) call fill_s_orbs(val, 2, 4)
1756 call fill_d_orbs(val, 7, 3) ! Co 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^7
1757 case (28)
1758 if (val > 8 ) call fill_s_orbs(val, 2, 4)
1759 call fill_d_orbs(val, 8, 3) ! Ni 1s^2 2s^2 2p^6 3s^2 3p^6 ; 4s^2 3d^8
1760 case (29)
1761 call fill_s_orbs(val, 1, 4) ! Cu 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^1
1762 case (30)
1763 call fill_s_orbs(val, 2, 4) ! Zn 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2
1764 case (31)
1765 call fill_p_orbs(val, 1, 4) ! Ga 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^1
1766 case (32)
1767 call fill_p_orbs(val, 2, 4) ! Ge 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^2
1768 case (33)
1769 call fill_p_orbs(val, 3, 4) ! As 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^3
1770 case (34)
1771 call fill_p_orbs(val, 4, 4) ! Se 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^4
1772 case (35)
1773 call fill_p_orbs(val, 5, 4) ! Br 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^5
1774 case (36)
1775 call fill_p_orbs(val, 6, 4) ! Kr 1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 ; 4s^2 4p^6
1776 case (37)
1777 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
1778 case (38)
1779 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
1780 case (39)
1781 if (val > 2) call fill_d_orbs(val, 1, 4)
1782 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
1783 case (40)
1784 if (val > 2) call fill_d_orbs(val, 2, 4)
1785 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
1786 case (41)
1787 if (val > 1) call fill_d_orbs(val, 4, 4)
1788 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
1789 case (42)
1790 if (val > 1) call fill_d_orbs(val, 5, 4)
1791 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
1792 case (43)
1793 if (val > 2) call fill_d_orbs(val, 5, 4)
1794 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
1795 case (44)
1796 if (val > 1) call fill_d_orbs(val, 7, 4)
1797 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
1798 case (45)
1799 if (val > 1) call fill_d_orbs(val, 8, 4)
1800 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
1801 case (46)
1802 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
1803 case (47)
1804 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
1805 case (48)
1806 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
1807 case (49)
1808 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
1809 case (50)
1810 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
1811 case (51)
1812 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
1813 case (52)
1814 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
1815 case (53)
1816 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
1817 case (54)
1818 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
1819 case (55)
1820 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
1821 case (56)
1822 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
1823 case (57)
1824 if (val > 2) call fill_d_orbs(val, 1, 5)
1825 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
1826 case (58)
1827 if (val > 3) call fill_f_orbs(val, 1, 4)
1828 if (val > 2) call fill_d_orbs(val, 1, 5)
1829 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
1830 case (59)
1831 if (val > 2) call fill_f_orbs(val, 3, 4)
1832 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
1833 case (60)
1834 if (val > 2) call fill_f_orbs(val, 4, 4)
1835 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
1836 case (61)
1837 if (val > 2) call fill_f_orbs(val, 5, 4)
1838 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
1839 case (62)
1840 if (val > 2) call fill_f_orbs(val, 6, 4)
1841 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
1842 case (63)
1843 if (val > 2) call fill_f_orbs(val, 7, 4)
1844 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
1845 case (64)
1846 if (val > 3) call fill_f_orbs(val, 7, 4)
1847 if (val > 2) call fill_d_orbs(val, 1, 5)
1848 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
1849 case (65)
1850 if (val > 2) call fill_f_orbs(val, 9, 4)
1851 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
1852 case (66)
1853 if (val > 2) call fill_f_orbs(val, 10, 4)
1854 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
1855 case (67)
1856 if (val > 2) call fill_f_orbs(val, 11, 4)
1857 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
1858 case (68)
1859 if (val > 2) call fill_f_orbs(val, 12, 4)
1860 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
1861 case (69)
1862 if (val > 2) call fill_f_orbs(val, 13, 4)
1863 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
1864 case (70)
1865 if (val > 2) call fill_f_orbs(val, 14, 4)
1866 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
1867 case (71)
1868 if (val > 3) call fill_f_orbs(val, 14, 4)
1869 if (val > 2) call fill_d_orbs(val, 1, 5)
1870 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
1871 case (72)
1872 if (val > 4) call fill_f_orbs(val, 14, 4)
1873 if (val > 2) call fill_d_orbs(val, 2, 5)
1874 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
1875 case (73)
1876 if (val > 5) call fill_f_orbs(val, 14, 4)
1877 if (val > 2) call fill_d_orbs(val, 3, 5)
1878 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
1879 case (74)
1880 if (val > 6) call fill_f_orbs(val, 14, 4)
1881 if (val > 2) call fill_d_orbs(val, 4, 5)
1882 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
1883 case (75)
1884 if (val > 7) call fill_f_orbs(val, 14, 4)
1885 if (val > 2) call fill_d_orbs(val, 5, 5)
1886 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
1887 case (76)
1888 if (val > 8) call fill_f_orbs(val, 14, 4)
1889 if (val > 2) call fill_d_orbs(val, 6, 5)
1890 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
1891 case (77)
1892 if (val > 9) call fill_f_orbs(val, 14, 4)
1893 if (val > 2) call fill_d_orbs(val, 7, 5)
1894 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
1895 case (78)
1896 if (val > 10) call fill_f_orbs(val, 14, 4)
1897 if (val > 1) call fill_d_orbs(val, 9, 5)
1898 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
1899 case (79)
1900 if (val > 21) call fill_f_orbs(val, 14, 4)
1901 if (val > 1) call fill_d_orbs(val, 10, 5)
1902 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
1903 case (80)
1904 if (val > 12) call fill_f_orbs(val, 14, 4)
1905 if (val > 2) call fill_d_orbs(val, 10, 5)
1906 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
1907 case (81)
1908 if (val > 10) call fill_d_orbs(val, 10, 5)
1909 if (val > 1) call fill_s_orbs(val, 2, 6)
1910 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
1911 case (82)
1912 if (val > 10) call fill_d_orbs(val, 10, 5)
1913 if (val > 2) call fill_s_orbs(val, 2, 6)
1914 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
1915 case (83)
1916 if (val > 10) call fill_d_orbs(val, 10, 5)
1917 if (val > 3) call fill_s_orbs(val, 2, 6)
1918 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
1919 case (84)
1920 if (val > 10) call fill_d_orbs(val, 10, 5)
1921 if (val > 4) call fill_s_orbs(val, 2, 6)
1922 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
1923 case (85)
1924 if (val > 10) call fill_d_orbs(val, 10, 5)
1925 if (val > 5) call fill_s_orbs(val, 2, 6)
1926 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
1927 case (86)
1928 if (val > 10) call fill_d_orbs(val, 10, 5)
1929 if (val > 6) call fill_s_orbs(val, 2, 6)
1930 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
1931 end select
1932
1933 !If we attributed all the electrons, everything went fine
1934 if (val < m_epsilon) then
1935 !In case of spin-polarized calculations, we properly distribute the electrons
1936 if (ispin == 2) then
1937 call valconf_unpolarized_to_polarized(conf)
1938 end if
1939 else
1940 conf%occ = m_zero
1941 message(1) = "Error in attributing atomic occupations"
1942 call messages_warning(1, namespace=namespace)
1943 end if
1944
1946
1947 contains
1948 subroutine fill_s_orbs(val, max_occ, nn)
1949 real(real64), intent(inout) :: val
1950 integer, intent(in) :: max_occ, nn
1951
1952 conf%p = conf%p + 1
1953 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1954 val = val - conf%occ(conf%p,1)
1955 conf%l(conf%p) = 0
1956 conf%n(conf%p) = nn
1957 end subroutine fill_s_orbs
1958
1959 subroutine fill_p_orbs(val, max_occ, nn)
1960 real(real64), intent(inout) :: val
1961 integer, intent(in) :: max_occ, nn
1962
1963 conf%p = conf%p + 1
1964 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1965 val = val - conf%occ(conf%p,1)
1966 conf%l(conf%p) = 1
1967 conf%n(conf%p) = nn
1968 end subroutine fill_p_orbs
1969
1970 subroutine fill_d_orbs(val, max_occ, nn)
1971 real(real64), intent(inout) :: val
1972 integer, intent(in) :: max_occ, nn
1973
1974 conf%p = conf%p + 1
1975 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1976 val = val - conf%occ(conf%p,1)
1977 conf%l(conf%p) = 2
1978 conf%n(conf%p) = nn
1979 end subroutine fill_d_orbs
1980
1981 subroutine fill_f_orbs(val, max_occ, nn)
1982 real(real64), intent(inout) :: val
1983 integer, intent(in) :: max_occ, nn
1984
1985 conf%p = conf%p + 1
1986 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1987 val = val - conf%occ(conf%p,1)
1988 conf%l(conf%p) = 3
1989 conf%n(conf%p) = nn
1990 end subroutine fill_f_orbs
1991
1992 end subroutine ps_guess_atomic_occupations
1993
1994#include "ps_pspio_inc.F90"
1995
1996end module ps_oct_m
1997
1998!! Local Variables:
1999!! mode: f90
2000!! coding: utf-8
2001!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:171
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:186
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:149
Some operations may be done for one spline-function, or for an array of them.
Definition: splines.F90:179
The integral may be done with or without integration limits, but we want the interface to be common.
Definition: splines.F90:173
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
real(real64), parameter, public m_min_exp_arg
Definition: global.F90:207
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
subroutine, public logrid_copy(grid_in, grid_out)
Definition: logrid.F90:228
integer function, public logrid_index(grid, rofi)
Definition: logrid.F90:256
subroutine, public logrid_end(grid)
Definition: logrid.F90:213
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
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:1472
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
subroutine, public messages_new_line()
Definition: messages.F90:1134
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
subroutine, public ps_cpi_end(ps_cpi)
Definition: ps_cpi.F90:183
subroutine, public ps_cpi_init(ps_cpi, filename, namespace)
Definition: ps_cpi.F90:148
subroutine, public ps_cpi_process(ps_cpi, lloc, namespace)
Definition: ps_cpi.F90:221
subroutine, public ps_fhi_end(ps_fhi)
Definition: ps_fhi.F90:183
subroutine, public ps_fhi_init(ps_fhi, filename, namespace)
Definition: ps_fhi.F90:149
subroutine, public ps_fhi_process(ps_fhi, lmax, lloc, namespace)
Definition: ps_fhi.F90:202
subroutine, public hgh_end(psp)
Definition: ps_hgh.F90:228
subroutine, public hgh_get_eigen(psp, eigen)
Definition: ps_hgh.F90:289
subroutine, public hgh_process(psp, namespace)
Definition: ps_hgh.F90:247
subroutine, public hgh_init(psp, filename, namespace)
Definition: ps_hgh.F90:183
real(real64) function, public first_point_extrapolate(x, y, high_order)
Definition: ps_in_grid.F90:437
Definition: ps.F90:114
pure logical function, public ps_has_density(ps)
Definition: ps.F90:1645
integer, parameter, public ps_filter_ts
Definition: ps.F90:162
subroutine, public ps_end(ps)
Definition: ps.F90:1059
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:261
subroutine ps_info(ps, filename, namespace)
Definition: ps.F90:581
subroutine, public ps_separate(ps)
Separate the local potential into (soft) long-ranged and (hard) short-ranged parts.
Definition: ps.F90:692
subroutine ps_check_bound(ps, eigen)
Definition: ps.F90:883
subroutine ps_grid_load(ps, ps_grid)
Definition: ps.F90:1173
subroutine, public ps_debug(ps, dir, namespace, gmax)
Definition: ps.F90:928
integer, parameter, public proj_hgh
Definition: ps.F90:167
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:1629
subroutine, public ps_getradius(ps)
Definition: ps.F90:760
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:1701
pure logical function, public ps_has_nlcc(ps)
Definition: ps.F90:1654
pure integer function, public ps_niwfs(ps)
Returns the number of atomic orbitals taking into account then m quantum number multiplicity.
Definition: ps.F90:1614
subroutine, public ps_pspio_init(ps, namespace, label, z, lmax, ispin, filename)
Definition: ps.F90:2108
integer, parameter, public proj_j_average
Fully-relativistic pseudopotentials with separate j-average and SOC terms.
Definition: ps.F90:173
integer, parameter, public proj_rkb
Definition: ps.F90:167
integer, parameter, public proj_j_dependent
Fully-relativistic j-dependent pseudopotentials.
Definition: ps.F90:173
subroutine, public ps_derivatives(ps)
Definition: ps.F90:780
subroutine, public ps_filter(ps, filter, gmax)
Definition: ps.F90:797
real(real64) function, public ps_density_volume(ps, namespace)
Definition: ps.F90:1662
integer, parameter, public ps_filter_bsb
Definition: ps.F90:162
integer, parameter, public proj_kb
Definition: ps.F90:167
subroutine ps_xml_load(ps, ps_xml, namespace)
Loads XML files for QSO, UPF1, UPF2, PSML, and PSP8 formats.
Definition: ps.F90:1276
subroutine, public ps_psf_init(pstm, ispin, filename, namespace)
Definition: ps_psf.F90:151
subroutine, public ps_psf_process(ps_psf, namespace, lmax, lloc)
Definition: ps_psf.F90:301
subroutine, public ps_psf_get_eigen(ps_psf, eigen)
Definition: ps_psf.F90:335
subroutine, public ps_psf_end(ps_psf)
Definition: ps_psf.F90:202
subroutine, public ps_xml_end(this)
Definition: ps_xml.F90:337
subroutine, public ps_xml_init(this, namespace, filename, fmt, ierr)
Definition: ps_xml.F90:163
integer, parameter, public pseudo_format_file_not_found
Definition: pseudo.F90:173
integer, parameter, public pseudo_type_kleinman_bylander
Definition: pseudo.F90:158
integer, parameter, public pseudo_format_qso
Definition: pseudo.F90:173
integer, parameter, public pseudo_format_hgh
Definition: pseudo.F90:173
integer, parameter, public pseudo_format_upf2
Definition: pseudo.F90:173
integer, parameter, public pseudo_type_paw
Definition: pseudo.F90:158
integer, parameter, public pseudo_format_cpi
Definition: pseudo.F90:173
integer, parameter, public pseudo_exchange_unknown
Definition: pseudo.F90:188
integer, parameter, public pseudo_type_semilocal
Definition: pseudo.F90:158
integer, parameter, public pseudo_type_ultrasoft
Definition: pseudo.F90:158
logical pure function, public pseudo_has_total_angular_momentum(pseudo)
Definition: pseudo.F90:548
integer, parameter, public pseudo_correlation_unknown
Definition: pseudo.F90:192
integer, parameter, public pseudo_format_psf
Definition: pseudo.F90:173
integer, parameter, public pseudo_format_fhi
Definition: pseudo.F90:173
integer, parameter, public pseudo_format_unknown
Definition: pseudo.F90:173
integer, parameter, public pseudo_format_psml
Definition: pseudo.F90:173
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
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:921
subroutine, public spline_force_pos(spl, threshold)
Definition: splines.F90:841
subroutine, public spline_fit(nrc, rofi, ffit, spl, threshold)
Definition: splines.F90:413
subroutine, public spline_3dft(spl, splw, threshold, gmax)
Definition: splines.F90:599
subroutine, public spline_besselft(spl, splw, l, threshold, gmax)
Definition: splines.F90:677
real(real64) function, public spline_eval(spl, x)
Definition: splines.F90:441
subroutine fill_s_orbs(val, max_occ, nn)
Definition: ps.F90:2042
subroutine fill_p_orbs(val, max_occ, nn)
Definition: ps.F90:2053
subroutine fill_d_orbs(val, max_occ, nn)
Definition: ps.F90:2064
subroutine get_splines()
Definition: ps.F90:1125
subroutine fill_f_orbs(val, max_occ, nn)
Definition: ps.F90:2075
remember that the FHI format is basically the CPI format with a header
Definition: ps_fhi.F90:137
The following data type contains: (a) the pseudopotential parameters, as read from a *....
Definition: ps_hgh.F90:144
A type storing the information and data about a pseudopotential.
Definition: ps.F90:184
the basic spline datatype
Definition: splines.F90:156
int true(void)