26 use,
intrinsic :: iso_fortran_env
69 integer,
parameter,
public :: &
74 integer,
public,
parameter :: &
80 integer,
public,
parameter :: &
81 PROJ_J_INDEPENDENT = 0, & !< Non-relativisitic or scalar-relativistic pseudopotentials
85 integer,
parameter,
public :: INVALID_L = 333
87 character(len=4),
parameter :: ps_name(PSEUDO_FORMAT_UPF1:PSEUDO_FORMAT_PSP8) = &
88 (/
"upf1",
"upf2",
"qso ",
"psml",
"psf ",
"cpi ",
"fhi ",
"hgh ",
"psp8"/)
93 integer :: projector_type
94 integer :: relativistic_treatment
97 character(len=10),
private :: label
99 integer,
private :: ispin
101 real(real64),
private :: z
102 real(real64) :: z_val
103 type(valconf_t) :: conf
105 type(logrid_t),
private :: g
107 type(spline_t),
allocatable :: ur(:, :)
108 type(spline_t),
allocatable,
private :: ur_sq(:, :)
109 logical,
allocatable :: bound(:, :)
116 logical :: no_vl = .false.
118 real(real64) :: projectors_sphere_threshold
125 real(real64) :: rc_max
128 integer :: projectors_per_l(5)
129 real(real64),
allocatable :: h(:,:,:), k(:, :, :)
130 type(spline_t),
allocatable :: kb(:, :)
131 type(spline_t),
allocatable :: dkb(:, :)
134 type(spline_t) :: core
135 type(spline_t) :: core_der
140 logical,
private :: has_long_range
141 logical,
private :: is_separated
143 type(spline_t),
private :: vlr
144 type(spline_t) :: vlr_sq
146 type(spline_t) :: nlr
148 real(real64) :: sigma_erf
150 logical,
private :: has_density
151 type(spline_t),
allocatable :: density(:)
152 type(spline_t),
allocatable :: density_der(:)
155 integer :: file_format
156 integer,
private :: pseudo_type
157 integer :: exchange_functional
158 integer :: correlation_functional
161 real(real64),
parameter :: eps = 1.0e-8_real64
167 subroutine ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
168 type(
ps_t),
intent(out) :: ps
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
177 integer :: l, ii, ll, is, ierr
183 real(real64),
allocatable :: eigen(:, :)
211 call parse_variable(namespace,
'SpeciesProjectorSphereThreshold', 0.001_real64, ps%projectors_sphere_threshold)
217 call messages_write(
"Cannot open pseudopotential file '"//trim(filename)//
"'.")
222 call messages_write(
"Cannot determine the pseudopotential type for species '"//trim(label)//
"' from", &
230 ps%relativistic_treatment = proj_j_independent
232 ps%sigma_erf = 0.625_real64
234 ps%projectors_per_l = 0
236 select case (ps%file_format)
238 ps%has_density = .
true.
240 ps%has_density = .false.
243 select case (ps%file_format)
249 call valconf_copy(ps%conf, ps_psf%conf)
254 ps%lmax = ps_psf%ps_grid%no_l_channels - 1
256 if (user_lmax /= invalid_l)
then
257 ps%lmax = min(ps%lmax, user_lmax)
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."
265 ps%conf%p = ps_psf%ps_grid%no_l_channels
266 if (ps%lmax == 0) ps%llocal = 0
269 if (user_llocal == invalid_l)
then
272 ps%llocal = user_llocal
275 ps%projectors_per_l(1:ps%lmax+1) = 1
276 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
285 call ps_cpi_init(ps_cpi, trim(filename), namespace)
286 ps%conf%p = ps_cpi%ps_grid%no_l_channels
288 call ps_fhi_init(ps_fhi, trim(filename), namespace)
289 ps%conf%p = ps_fhi%ps_grid%no_l_channels
293 ps%conf%symbol = label(1:2)
302 ps%lmax = ps%conf%p - 1
304 if (user_lmax /= invalid_l)
then
305 ps%lmax = min(ps%lmax, user_lmax)
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."
313 if (ps%lmax == 0) ps%llocal = 0
316 if (user_llocal == invalid_l)
then
319 ps%llocal = user_llocal
322 ps%projectors_per_l(1:ps%lmax+1) = 1
323 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
337 call hgh_init(ps_hgh, trim(filename), namespace)
338 call valconf_copy(ps%conf, ps_hgh%conf)
341 ps%z_val = ps_hgh%z_val
344 ps%lmax = ps_hgh%l_max
345 ps%conf%symbol = label(1:2)
346 ps%sigma_erf = ps_hgh%rlocal
348 ps%projectors_per_l(1:ps%lmax+1) = 1
356 call valconf_copy(ps_hgh%conf, ps%conf)
363 call valconf_unpolarized_to_polarized(ps%conf)
368 call ps_xml_init(ps_xml, namespace, trim(filename), ps%file_format, ierr)
377 if (ps_xml%kleinman_bylander)
then
378 ps%conf%p = ps_xml%nwavefunctions
380 ps%conf%p = ps_xml%lmax + 1
383 do ll = 0, ps_xml%lmax
384 ps%conf%l(ll + 1) = ll
387 ps%kbc = ps_xml%nchannels
388 ps%lmax = ps_xml%lmax
390 ps%projectors_per_l = 0
391 do ll = 0, ps_xml%lmax
395 if (ps_xml%kleinman_bylander)
then
396 ps%llocal = ps_xml%llocal
400 if (ps_xml%llocal >= 0) ps%llocal = ps_xml%llocal
401 if (user_llocal /= invalid_l) ps%llocal = user_llocal
402 assert(ps%llocal >= 0)
403 assert(ps%llocal <= ps%lmax)
406 ps%g%nrval = ps_xml%grid_size
408 safe_allocate(ps%g%rofi(1:ps%g%nrval))
409 safe_allocate(ps%g%r2ofi(1:ps%g%nrval))
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
418 ps%local = (ps%lmax == 0 .and. ps%llocal == 0 ) .or. (ps%lmax == -1 .and. ps%llocal == -1)
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))
438 safe_allocate(eigen(1:ps%conf%p, 1:ps%ispin))
442 select case (ps%file_format)
455 safe_allocate(ps%k (0:ps%lmax, 1:ps%kbc, 1:ps%kbc))
465 call spline_der(ps%density(is), ps%density_der(is), ps%projectors_sphere_threshold)
470 call spline_der(ps%core, ps%core_der, ps%projectors_sphere_threshold)
475 ps%has_long_range = .
true.
476 ps%is_separated = .false.
478 call ps_info(ps, filename, namespace)
480 safe_deallocate_a(eigen)
487 subroutine ps_info(ps, filename, namespace)
488 type(
ps_t),
intent(in) :: ps
489 character(len=*),
intent(in) :: filename
498 select case (ps%file_format)
499 case (pseudo_format_upf1)
507 case (pseudo_format_psp8)
529 select case (ps%pseudo_type)
543 select case (ps%file_format)
557 if (ps%llocal >= 0)
then
569 if (ps%llocal < 0)
then
599 type(
ps_t),
intent(inout) :: ps
601 real(real64),
allocatable :: vsr(:), vlr(:), nlr(:)
602 real(real64) :: r, exp_arg, arg, max_val_vsr
607 assert(ps%g%nrval > 0)
609 safe_allocate(vsr(1:ps%g%nrval))
610 safe_allocate(vlr(1:ps%g%nrval))
611 safe_allocate(nlr(1:ps%g%nrval))
617 do ii = 1, ps%g%nrval
622 if(arg > 1e-3_real64)
then
629 if(abs(vsr(ii)) < 1.0e-14_real64) vsr(ii) =
m_zero
630 max_val_vsr = max(max_val_vsr, abs(vsr(ii)))
632 exp_arg = -
m_half*r**2/ps%sigma_erf**2
641 call spline_fit(ps%g%nrval, ps%g%rofi, vlr, ps%vlr, ps%projectors_sphere_threshold)
644 call spline_fit(ps%g%nrval, ps%g%r2ofi, vlr, ps%vlr_sq, ps%projectors_sphere_threshold)
647 call spline_fit(ps%g%nrval, ps%g%rofi, nlr, ps%nlr, ps%projectors_sphere_threshold)
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.
655 safe_deallocate_a(vsr)
656 safe_deallocate_a(vlr)
657 safe_deallocate_a(nlr)
659 ps%is_separated = .
true.
667 type(
ps_t),
intent(inout) :: ps
675 if (l == ps%llocal) cycle
677 ps%rc_max = max(ps%rc_max, ps%kb(l, j)%x_threshold)
687 type(
ps_t),
intent(inout) :: ps
694 call spline_der(ps%kb(l, j), ps%dkb(l, j), ps%projectors_sphere_threshold)
704 type(
ps_t),
intent(inout) :: ps
705 integer,
intent(in) :: filter
706 real(real64),
intent(in) :: gmax
708 integer :: l, k, ispin
710 real(real64) :: alpha, beta_fs, rmax, rcut, gamma, beta_rs
716 case (ps_filter_none)
722 if(.not. ps%no_vl)
then
723 rmax = ps%vl%x_threshold
725 call spline_filter_mask(ps%vl, max(0, ps%llocal), rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
729 if (l == ps%llocal) cycle
731 call spline_filter_mask(ps%kb(l, k), l, ps%rc_max, gmax, alpha, gamma, ps%projectors_sphere_threshold)
736 rmax = ps%core%x_threshold
737 call spline_filter_mask(ps%core, 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
741 do ispin = 1, ps%ispin
743 rmax = ps%density(ispin)%x_threshold
744 call spline_filter_mask(ps%density(ispin), 0, rmax, gmax, alpha, gamma, ps%projectors_sphere_threshold)
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)
757 beta_fs = 18.0_real64
762 call spline_filter_bessel(ps%vl, max(0, ps%llocal), gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
764 if (l == ps%llocal) cycle
766 call spline_filter_bessel(ps%kb(l, k), l, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
771 call spline_filter_bessel(ps%core, 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
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)
778 call spline_filter_bessel(ps%density_der(ispin), 0, gmax, alpha, beta_fs, rcut, beta_rs, ps%projectors_sphere_threshold)
790 type(
ps_t),
intent(inout) :: ps
791 real(real64),
intent(in) :: eigen(:,:)
794 real(real64) :: ur1, ur2
808 if (.not. ps%bound(i, is)) cycle
810 do ir = ps%g%nrval, 3, -1
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.
834 subroutine ps_debug(ps, dir, namespace, gmax)
835 type(
ps_t),
intent(in) :: ps
836 character(len=*),
intent(in) :: dir
838 real(real64),
intent(in) :: gmax
841 type(
spline_t),
allocatable :: fw(:, :)
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:'
861 write(iunit,
'(10f9.5)') (ps%h(l, k, j), j = 1, ps%kbc)
864 if (
allocated(ps%k))
then
865 write(iunit,
'(/,a,/)')
'k matrix:'
868 write(iunit,
'(10f9.5)') (ps%k(l, k, j), j = 1, ps%kbc)
873 write(iunit,
'(/,a)')
'orbitals:'
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)
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)
890 iunit =
io_open(trim(dir)//
'/local', namespace, action=
'write')
895 iunit =
io_open(trim(dir)//
'/local_long_range', namespace, action=
'write')
900 iunit =
io_open(trim(dir)//
'/local_long_range_density', namespace, action=
'write')
905 iunit =
io_open(trim(dir)//
'/local_ft', namespace, action=
'write')
906 safe_allocate(fw(1, 1))
908 call spline_3dft(ps%vl, fw(1, 1), ps%projectors_sphere_threshold, gmax = gmax)
911 safe_deallocate_a(fw)
915 iunit =
io_open(trim(dir)//
'/nonlocal', namespace, action=
'write')
919 iunit =
io_open(trim(dir)//
'/nonlocal_derivative', namespace, action=
'write')
923 iunit =
io_open(trim(dir)//
'/nonlocal_ft', namespace, action=
'write')
924 safe_allocate(fw(0:ps%lmax, 1:ps%kbc))
928 call spline_besselft(ps%kb(k, j), fw(k, j), k, threshold=ps%projectors_sphere_threshold, &
934 safe_deallocate_a(fw)
938 iunit =
io_open(trim(dir)//
'/wavefunctions', namespace, action=
'write')
943 if (ps%has_density)
then
944 iunit =
io_open(trim(dir)//
'/density', namespace, action=
'write')
948 iunit =
io_open(trim(dir)//
'/density_derivative', namespace, action=
'write')
955 iunit =
io_open(trim(dir)//
'/nlcc', namespace, action=
'write')
966 type(
ps_t),
intent(inout) :: ps
968 if (.not.
allocated(ps%kb))
return
972 if (ps%is_separated)
then
987 if (
allocated(ps%density))
call spline_end(ps%density)
988 if (
allocated(ps%density_der))
call spline_end(ps%density_der)
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)
1008 type(
ps_t),
intent(inout) :: ps
1009 type(
ps_hgh_t),
intent(inout) :: ps_hgh
1015 if (ps%lmax >= 0)
then
1016 ps%rc_max = 1.1_real64 * maxval(ps_hgh%kbr(0:ps%lmax))
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)
1032 integer :: l, is, j, ip
1033 real(real64),
allocatable :: hato(:), dens(:)
1037 safe_allocate(hato(1:ps_hgh%g%nrval))
1038 safe_allocate(dens(1:ps_hgh%g%nrval))
1041 do l = 0, ps_hgh%l_max
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)
1050 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, ps_hgh%vlocal, ps%vl, ps%projectors_sphere_threshold)
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)
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)
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)
1067 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1070 safe_deallocate_a(hato)
1071 safe_deallocate_a(dens)
1080 type(
ps_t),
intent(inout) :: ps
1086 ps%z_val = ps_grid%zval
1088 ps%nlcc = ps_grid%core_corrections
1090 ps%h(0:ps%lmax, 1, 1) = ps_grid%dkbcos(1:ps%lmax+1)
1094 ps%rc_max = maxval(ps_grid%kb_radius(1:ps%lmax+1)) * 1.5_real64
1100 ps%h(0:ps%lmax,:,:) = ps%h(0:ps%lmax,:,:) /
m_two
1109 real(real64),
allocatable :: hato(:), dens(:)
1110 integer :: is, l, ir, nrc, ip
1114 safe_allocate(hato(1:g%nrval))
1115 safe_allocate(dens(1:g%nrval))
1122 do l = 1, ps_grid%no_l_channels
1123 hato(2:) = ps_grid%rphi(2:, l, 1+is)/g%rofi(2:)
1128 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(
m_four*
m_pi)
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)
1136 call spline_fit(g%nrval, g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1143 hato(1:nrc) = ps_grid%KB(1:nrc, l)
1144 hato(nrc+1:g%nrval) =
m_zero
1146 call spline_fit(g%nrval, g%rofi, hato, ps%kb(l-1, 1), ps%projectors_sphere_threshold)
1151 hato(:) = ps_grid%vlocal(:)/
m_two
1152 call spline_fit(g%nrval, g%rofi, hato, ps%vl, ps%projectors_sphere_threshold)
1154 if (ps_grid%core_corrections)
then
1156 hato(2:) = ps_grid%chcore(2:)/(
m_four*
m_pi*g%r2ofi(2:))
1158 do ir = g%nrval-1, 2, -1
1159 if (hato(ir) > eps)
then
1165 hato(nrc:g%nrval) =
m_zero
1168 call spline_fit(g%nrval, g%rofi, hato, ps%core, ps%projectors_sphere_threshold)
1171 safe_deallocate_a(hato)
1172 safe_deallocate_a(dens)
1183 type(
ps_t),
intent(inout) :: ps
1184 type(
ps_xml_t),
intent(in) :: ps_xml
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
1201 if (ps%file_format == pseudo_format_psp8)
then
1203 message(1) =
"SOC from PSP8 is not currently supported."
1204 message(2) =
"Only scalar relativistic effects will be considered."
1211 ps%nlcc = ps_xml%nlcc
1213 ps%z_val = ps_xml%valence_charge
1215 ps%has_density = ps_xml%has_density
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)
1230 vlocal(ip) = -ps_xml%valence_charge/rr
1235 call spline_fit(ps%g%nrval, ps%g%rofi, vlocal, ps%vl, ps%projectors_sphere_threshold)
1237 safe_deallocate_a(vlocal)
1239 safe_allocate(kbprojector(1:ps%g%nrval))
1240 safe_allocate(wavefunction(1:ps%g%nrval))
1245 density_is_known = .false.
1248 if (ps_xml%kleinman_bylander)
then
1250 safe_allocate(cmap(0:ps_xml%lmax, 1:ps_xml%nchannels))
1254 do ll = 0, ps_xml%lmax
1255 do ic = 1, ps_xml%nchannels
1259 if (ll == ps_xml%llocal) cycle
1266 assert(mod(ps_xml%nchannels, 2)==0)
1268 cmap(ll, ic) = int((ic-1)/2)*2 + 2
1271 cmap(ll, ic) = int((ic-1)/2)*2 + 1
1277 assert(sum(cmap(ll, 1:ps_xml%nchannels)) == (ps_xml%nchannels + 1)*ps_xml%nchannels/2)
1280 assert(all(cmap >= 0 .and. cmap <= ps_xml%nchannels))
1290 safe_allocate(matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels))
1291 safe_allocate(eigenvalues(1:ps_xml%nchannels))
1293 do ll = 0, ps_xml%lmax
1295 if (
is_diagonal(ps_xml%nchannels, ps_xml%dij(ll, :, :), tol_diagonal) .or. &
1298 do ic = 1, ps_xml%nchannels
1299 eigenvalues(ic) = ps_xml%dij(ll, ic, ic)
1300 matrix(ic, ic) =
m_one
1304 matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels) = ps_xml%dij(ll, 1:ps_xml%nchannels, 1:ps_xml%nchannels)
1308 do ic = 1, ps_xml%nchannels
1310 do jc = 1, ps_xml%nchannels
1311 call lalg_axpy(ps_xml%grid_size, matrix(jc, ic), ps_xml%projector(:, ll, jc), kbprojector)
1314 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, cmap(ll, ic)), ps%projectors_sphere_threshold)
1316 ps%h(ll, cmap(ll, ic), cmap(ll, ic)) = eigenvalues(ic)
1321 safe_deallocate_a(matrix)
1322 safe_deallocate_a(eigenvalues)
1325 ps%conf%p = ps_xml%nwavefunctions
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))
1337 do ii = 1, ps_xml%nwavefunctions
1339 ps%conf%n(ii) = ps_xml%wf_n(ii)
1340 ps%conf%l(ii) = ps_xml%wf_l(ii)
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)
1346 ps%conf%occ(ii, 1) = ps_xml%wf_occ(ii)
1354 do ip = 1, ps%g%nrval
1355 if (ip <= ps_xml%grid_size)
then
1356 wavefunction(ip) = ps_xml%wavefunction(ip, ii)
1358 wavefunction(ip) =
m_zero
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)
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)
1378 if ((.not.
ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0)
then
1380 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1382 safe_deallocate_a(dens)
1383 density_is_known = .
true.
1384 ps%has_density = .
true.
1387 safe_deallocate_a(cmap)
1393 ps%conf%symbol = ps%label(1:2)
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)
1404 safe_deallocate_a(order)
1405 safe_deallocate_a(occ_tmp)
1408 if (abs(sum(ps%conf%occ) - ps%z_val ) <
m_epsilon)
then
1409 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1413 do ll = 0, ps_xml%lmax
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
1426 kbcos = dnrm/(safe_tol(avgv,1.0e-20_real64))
1427 kbnorm =
m_one/(safe_tol(
sqrt(dnrm),1.0e-20_real64))
1429 if (ll /= ps%llocal)
then
1430 ps%h(ll, 1, 1) = kbcos
1431 kbprojector = kbprojector*kbnorm
1436 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, 1), ps%projectors_sphere_threshold)
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
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)
1448 if (abs(sum(ps%conf%occ) - ps%z_val) <
m_epsilon)
then
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)
1458 if (abs(sum(ps%conf%occ) - ps%z_val) <
m_epsilon)
then
1460 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1462 safe_deallocate_a(dens)
1463 density_is_known = .
true.
1464 ps%has_density = .
true.
1471 safe_allocate(dens(1:ps%g%nrval, 1))
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
1477 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:, 1), ps%density(is), ps%projectors_sphere_threshold)
1480 safe_deallocate_a(dens)
1485 if (ps_xml%nlcc)
then
1487 safe_allocate(nlcc_density(1:ps%g%nrval))
1489 call lalg_copy(ps_xml%grid_size, ps_xml%nlcc_density, nlcc_density)
1492 do ir = ps_xml%grid_size - 1, 1, -1
1493 if (nlcc_density(ir) > eps)
then
1499 nlcc_density(nrc:ps%g%nrval) =
m_zero
1501 call spline_fit(ps%g%nrval, ps%g%rofi, nlcc_density, ps%core, ps%projectors_sphere_threshold)
1503 safe_deallocate_a(nlcc_density)
1509 ps%rc_max = ps%rc_max*1.05_real64
1511 safe_deallocate_a(kbprojector)
1512 safe_deallocate_a(wavefunction)
1521 type(
ps_t),
intent(in) :: ps
1536 type(
ps_t),
intent(in) :: ps
1543 if (any(.not. ps%bound(i,:))) cycle
1552 type(
ps_t),
intent(in) :: ps
1554 has_density = ps%has_density
1560 pure logical function ps_has_nlcc(ps)
result(has_nlcc)
1561 type(
ps_t),
intent(in) :: ps
1568 real(real64) function ps_density_volume(ps, namespace) result(volume)
1569 type(
ps_t),
intent(in) :: ps
1572 integer :: ip, ispin
1574 real(real64),
allocatable ::vol(:)
1577 push_sub(ps_density_volume)
1580 message(1) =
"The pseudopotential does not contain an atomic density"
1584 safe_allocate(vol(1:ps%g%nrval))
1586 do ip = 1, ps%g%nrval
1589 do ispin = 1, ps%ispin
1595 call spline_fit(ps%g%nrval, ps%g%rofi, vol, volspl, ps%projectors_sphere_threshold)
1599 safe_deallocate_a(vol)
1601 pop_sub(ps_density_volume)
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
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.)
1627 assert(valcharge <= zz)
1631 if(int(zz) > 2 .and. val > zz - 2)
then
1635 if(int(zz) > 4 .and. val > zz - 4)
then
1640 if(int(zz) > 18 .and. val > zz - 10)
then
1644 if(int(zz) > 12 .and. val > zz - 12)
then
1648 if(int(zz) > 18 .and. val > zz - 18)
then
1652 if(int(zz) > 28 .and. val > zz - 28)
then
1656 if(int(zz) > 30 .and. val > zz - 30)
then
1660 if(int(zz) > 36 .and. val > zz - 36)
then
1664 if(int(zz) > 46 .and. val > zz - 46)
then
1669 if((int(zz) > 48 .and. val > zz - 48) .or. &
1670 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 62))
then
1675 if((int(zz) > 54 .and. val > zz - 54) .or. &
1676 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 68) )
then
1681 if(int(zz) > 80 .and. val > zz - 68)
then
1687 select case (int(zz))
1934 if (val < m_epsilon)
then
1936 if (ispin == 2)
then
1937 call valconf_unpolarized_to_polarized(conf)
1941 message(1) =
"Error in attributing atomic occupations"
1942 call messages_warning(1, namespace=namespace)
1949 real(real64),
intent(inout) :: val
1950 integer,
intent(in) :: max_occ, nn
1953 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1954 val = val - conf%occ(conf%p,1)
1960 real(real64),
intent(inout) :: val
1961 integer,
intent(in) :: max_occ, nn
1964 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1965 val = val - conf%occ(conf%p,1)
1971 real(real64),
intent(inout) :: val
1972 integer,
intent(in) :: max_occ, nn
1975 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1976 val = val - conf%occ(conf%p,1)
1982 real(real64),
intent(inout) :: val
1983 integer,
intent(in) :: max_occ, nn
1986 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
1987 val = val - conf%occ(conf%p,1)
1994#include "ps_pspio_inc.F90"
constant times a vector plus a vector
Copies a vector x, to a vector y.
This is the common interface to a sorting routine. It performs the shell algorithm,...
Some operations may be done for one spline-function, or for an array of them.
The integral may be done with or without integration limits, but we want the interface to be common.
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_four
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public m_min_exp_arg
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_three
subroutine, public io_close(iunit, grp)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
subroutine, public logrid_copy(grid_in, grid_out)
integer function, public logrid_index(grid, rofi)
subroutine, public logrid_end(grid)
This module is intended to contain "only mathematical" functions and procedures.
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.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
subroutine, public messages_new_line()
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
subroutine, public ps_cpi_end(ps_cpi)
subroutine, public ps_cpi_init(ps_cpi, filename, namespace)
subroutine, public ps_cpi_process(ps_cpi, lloc, namespace)
subroutine, public ps_fhi_end(ps_fhi)
subroutine, public ps_fhi_init(ps_fhi, filename, namespace)
subroutine, public ps_fhi_process(ps_fhi, lmax, lloc, namespace)
subroutine, public hgh_end(psp)
subroutine, public hgh_get_eigen(psp, eigen)
subroutine, public hgh_process(psp, namespace)
subroutine, public hgh_init(psp, filename, namespace)
real(real64) function, public first_point_extrapolate(x, y, high_order)
pure logical function, public ps_has_density(ps)
integer, parameter, public ps_filter_ts
subroutine, public ps_end(ps)
subroutine hgh_load(ps, ps_hgh)
subroutine, public ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
subroutine ps_info(ps, filename, namespace)
subroutine, public ps_separate(ps)
Separate the local potential into (soft) long-ranged and (hard) short-ranged parts.
subroutine ps_check_bound(ps, eigen)
subroutine ps_grid_load(ps, ps_grid)
subroutine, public ps_debug(ps, dir, namespace, gmax)
integer, parameter, public proj_hgh
pure integer function, public ps_bound_niwfs(ps)
Returns the number of bound atomic orbitals taking into account then m quantum number multiplicity.
subroutine, public ps_getradius(ps)
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....
pure logical function, public ps_has_nlcc(ps)
pure integer function, public ps_niwfs(ps)
Returns the number of atomic orbitals taking into account then m quantum number multiplicity.
subroutine, public ps_pspio_init(ps, namespace, label, z, lmax, ispin, filename)
integer, parameter, public proj_j_average
Fully-relativistic pseudopotentials with separate j-average and SOC terms.
integer, parameter, public proj_rkb
integer, parameter, public proj_j_dependent
Fully-relativistic j-dependent pseudopotentials.
subroutine, public ps_derivatives(ps)
subroutine, public ps_filter(ps, filter, gmax)
real(real64) function, public ps_density_volume(ps, namespace)
integer, parameter, public ps_filter_bsb
integer, parameter, public proj_kb
subroutine ps_xml_load(ps, ps_xml, namespace)
Loads XML files for QSO, UPF1, UPF2, PSML, and PSP8 formats.
subroutine, public ps_psf_init(pstm, ispin, filename, namespace)
subroutine, public ps_psf_process(ps_psf, namespace, lmax, lloc)
subroutine, public ps_psf_get_eigen(ps_psf, eigen)
subroutine, public ps_psf_end(ps_psf)
subroutine, public ps_xml_end(this)
subroutine, public ps_xml_init(this, namespace, filename, fmt, ierr)
integer, parameter, public pseudo_format_file_not_found
integer, parameter, public pseudo_type_kleinman_bylander
integer, parameter, public pseudo_format_qso
integer, parameter, public pseudo_format_hgh
integer, parameter, public pseudo_format_upf2
integer, parameter, public pseudo_type_paw
integer, parameter, public pseudo_format_cpi
integer, parameter, public pseudo_exchange_unknown
integer, parameter, public pseudo_type_semilocal
integer, parameter, public pseudo_type_ultrasoft
logical pure function, public pseudo_has_total_angular_momentum(pseudo)
integer, parameter, public pseudo_correlation_unknown
integer, parameter, public pseudo_format_psf
integer, parameter, public pseudo_format_fhi
integer, parameter, public pseudo_format_unknown
integer, parameter, public pseudo_format_psml
This module is intended to contain "only mathematical" functions and procedures.
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)
subroutine, public spline_force_pos(spl, threshold)
subroutine, public spline_fit(nrc, rofi, ffit, spl, threshold)
subroutine, public spline_3dft(spl, splw, threshold, gmax)
subroutine, public spline_besselft(spl, splw, l, threshold, gmax)
real(real64) function, public spline_eval(spl, x)
subroutine fill_s_orbs(val, max_occ, nn)
subroutine fill_p_orbs(val, max_occ, nn)
subroutine fill_d_orbs(val, max_occ, nn)
subroutine fill_f_orbs(val, max_occ, nn)
remember that the FHI format is basically the CPI format with a header
The following data type contains: (a) the pseudopotential parameters, as read from a *....
A type storing the information and data about a pseudopotential.
the basic spline datatype