26 use,
intrinsic :: iso_fortran_env
71 integer,
parameter,
public :: &
76 integer,
public,
parameter :: &
82 integer,
public,
parameter :: &
83 PROJ_J_INDEPENDENT = 0, & !< Non-relativisitic or scalar-relativistic pseudopotentials
87 integer,
parameter,
public :: INVALID_L = 333
89 character(len=4),
parameter :: ps_name(PSEUDO_FORMAT_UPF1:PSEUDO_FORMAT_PSP8) = &
90 (/
"upf1",
"upf2",
"qso ",
"psml",
"psf ",
"cpi ",
"fhi ",
"hgh ",
"psp8"/)
95 integer :: projector_type
96 integer :: relativistic_treatment
99 character(len=10),
private :: label
101 integer,
private :: ispin
103 real(real64),
private :: z
104 real(real64) :: z_val
105 type(valconf_t) :: conf
107 type(logrid_t),
private :: g
109 type(spline_t),
allocatable :: ur(:, :)
110 logical,
allocatable :: bound(:, :)
117 logical :: no_vl = .false.
119 real(real64) :: projectors_sphere_threshold
126 real(real64) :: rc_max
129 integer :: projectors_per_l(5)
130 real(real64),
allocatable :: h(:,:,:), k(:, :, :)
131 type(spline_t),
allocatable :: kb(:, :)
132 type(spline_t),
allocatable :: dkb(:, :)
135 type(spline_t) :: core
136 type(spline_t) :: core_der
141 logical,
private :: has_long_range
142 logical,
private :: is_separated
144 type(spline_t),
private :: vlr
145 type(spline_t) :: vlr_sq
147 type(spline_t) :: nlr
149 real(real64) :: sigma_erf
151 logical,
private :: has_density
152 type(spline_t),
allocatable :: density(:)
153 type(spline_t),
allocatable :: density_der(:)
156 integer :: file_format
157 integer,
private :: pseudo_type
158 integer :: exchange_functional
159 integer :: correlation_functional
162 real(real64),
parameter :: eps = 1.0e-8_real64
168 subroutine ps_init(ps, namespace, label, z, user_lmax, user_llocal, ispin, filename)
169 type(
ps_t),
intent(out) :: ps
171 character(len=10),
intent(in) :: label
172 integer,
intent(in) :: user_lmax
173 integer,
intent(in) :: user_llocal
174 integer,
intent(in) :: ispin
175 real(real64),
intent(in) :: z
176 character(len=*),
intent(in) :: filename
178 integer :: l, ii, ll, is, ierr
184 real(real64),
allocatable :: eigen(:, :)
212 call parse_variable(namespace,
'SpeciesProjectorSphereThreshold', 0.001_real64, ps%projectors_sphere_threshold)
218 call messages_write(
"Cannot open pseudopotential file '"//trim(filename)//
"'.")
223 call messages_write(
"Cannot determine the pseudopotential type for species '"//trim(label)//
"' from", &
231 ps%relativistic_treatment = proj_j_independent
233 ps%sigma_erf = 0.625_real64
235 ps%projectors_per_l = 0
237 select case (ps%file_format)
241 ps%has_density = .false.
244 select case (ps%file_format)
255 ps%lmax = ps_psf%ps_grid%no_l_channels - 1
257 if (user_lmax /= invalid_l)
then
258 ps%lmax = min(ps%lmax, user_lmax)
259 if (user_lmax /= ps%lmax)
then
260 message(1) =
"lmax in Species block for " // trim(label) // &
261 " is larger than number available in pseudopotential."
266 ps%conf%p = ps_psf%ps_grid%no_l_channels
267 if (ps%lmax == 0) ps%llocal = 0
270 if (user_llocal == invalid_l)
then
273 ps%llocal = user_llocal
276 ps%projectors_per_l(1:ps%lmax+1) = 1
277 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
286 call ps_cpi_init(ps_cpi, trim(filename), namespace)
287 ps%conf%p = ps_cpi%ps_grid%no_l_channels
289 call ps_fhi_init(ps_fhi, trim(filename), namespace)
290 ps%conf%p = ps_fhi%ps_grid%no_l_channels
294 ps%conf%symbol = label(1:2)
303 ps%lmax = ps%conf%p - 1
305 if (user_lmax /= invalid_l)
then
306 ps%lmax = min(ps%lmax, user_lmax)
307 if (user_lmax /= ps%lmax)
then
308 message(1) =
"lmax in Species block for " // trim(label) // &
309 " is larger than number available in pseudopotential."
314 if (ps%lmax == 0) ps%llocal = 0
317 if (user_llocal == invalid_l)
then
320 ps%llocal = user_llocal
323 ps%projectors_per_l(1:ps%lmax+1) = 1
324 if (ps%llocal > -1) ps%projectors_per_l(ps%llocal+1) = 0
338 call hgh_init(ps_hgh, trim(filename), namespace)
342 ps%z_val = ps_hgh%z_val
345 ps%lmax = ps_hgh%l_max
346 ps%conf%symbol = label(1:2)
347 ps%sigma_erf = ps_hgh%rlocal
349 ps%projectors_per_l(1:ps%lmax+1) = 1
369 call ps_xml_init(ps_xml, namespace, trim(filename), ps%file_format, ierr)
378 if (ps_xml%kleinman_bylander)
then
379 ps%conf%p = ps_xml%nwavefunctions
381 ps%conf%p = ps_xml%lmax + 1
384 do ll = 0, ps_xml%lmax
385 ps%conf%l(ll + 1) = ll
388 ps%kbc = ps_xml%nchannels
389 ps%lmax = ps_xml%lmax
391 ps%projectors_per_l = 0
392 do ll = 0, ps_xml%lmax
396 if (ps_xml%kleinman_bylander)
then
397 ps%llocal = ps_xml%llocal
401 if (ps_xml%llocal >= 0) ps%llocal = ps_xml%llocal
402 if (user_llocal /= invalid_l) ps%llocal = user_llocal
403 assert(ps%llocal >= 0)
404 assert(ps%llocal <= ps%lmax)
407 ps%g%nrval = ps_xml%grid_size
409 safe_allocate(ps%g%rofi(1:ps%g%nrval))
410 safe_allocate(ps%g%r2ofi(1:ps%g%nrval))
412 do ii = 1, ps%g%nrval
413 ps%g%rofi(ii) = ps_xml%grid(ii)
414 ps%g%r2ofi(ii) = ps_xml%grid(ii)**2
419 ps%local = (ps%lmax == 0 .and. ps%llocal == 0 ) .or. (ps%lmax == -1 .and. ps%llocal == -1)
422 safe_allocate(ps%kb (0:ps%lmax, 1:ps%kbc))
423 safe_allocate(ps%dkb (0:ps%lmax, 1:ps%kbc))
424 safe_allocate(ps%ur (1:ps%conf%p, 1:ps%ispin))
425 safe_allocate(ps%bound(1:ps%conf%p, 1:ps%ispin))
426 safe_allocate(ps%h (0:ps%lmax, 1:ps%kbc, 1:ps%kbc))
427 safe_allocate(ps%density(1:ps%ispin))
428 safe_allocate(ps%density_der(1:ps%ispin))
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
986 if (
allocated(ps%density))
call spline_end(ps%density)
987 if (
allocated(ps%density_der))
call spline_end(ps%density_der)
991 safe_deallocate_a(ps%kb)
992 safe_deallocate_a(ps%dkb)
993 safe_deallocate_a(ps%ur)
994 safe_deallocate_a(ps%bound)
995 safe_deallocate_a(ps%h)
996 safe_deallocate_a(ps%k)
997 safe_deallocate_a(ps%density)
998 safe_deallocate_a(ps%density_der)
1006 type(
ps_t),
intent(inout) :: ps
1007 type(
ps_hgh_t),
intent(inout) :: ps_hgh
1013 if (ps%lmax >= 0)
then
1014 ps%rc_max = 1.1_real64 * maxval(ps_hgh%kbr(0:ps%lmax))
1018 ps%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%h(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
1019 ps%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc) = ps_hgh%k(0:ps%lmax, 1:ps%kbc, 1:ps%kbc)
1030 integer :: l, is, j, ip
1031 real(real64),
allocatable :: hato(:), dens(:)
1035 safe_allocate(hato(1:ps_hgh%g%nrval))
1036 safe_allocate(dens(1:ps_hgh%g%nrval))
1039 do l = 0, ps_hgh%l_max
1041 hato(:) = ps_hgh%kb(:, l, j)
1042 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%kb(l, j), ps%projectors_sphere_threshold)
1048 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, ps_hgh%vlocal, ps%vl, ps%projectors_sphere_threshold)
1055 hato(2:ps_hgh%g%nrval) = ps_hgh%rphi(2:ps_hgh%g%nrval, l)/ps_hgh%g%rofi(2:ps_hgh%g%nrval)
1058 do ip = 1, ps_hgh%g%nrval
1059 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(
m_four*
m_pi)
1062 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
1064 call spline_fit(ps_hgh%g%nrval, ps_hgh%g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1067 safe_deallocate_a(hato)
1068 safe_deallocate_a(dens)
1077 type(
ps_t),
intent(inout) :: ps
1083 ps%z_val = ps_grid%zval
1085 ps%nlcc = ps_grid%core_corrections
1087 ps%h(0:ps%lmax, 1, 1) = ps_grid%dkbcos(1:ps%lmax+1)
1091 ps%rc_max = maxval(ps_grid%kb_radius(1:ps%lmax+1)) * 1.5_real64
1097 ps%h(0:ps%lmax,:,:) = ps%h(0:ps%lmax,:,:) /
m_two
1106 real(real64),
allocatable :: hato(:), dens(:)
1107 integer :: is, l, ir, nrc, ip
1111 safe_allocate(hato(1:g%nrval))
1112 safe_allocate(dens(1:g%nrval))
1119 do l = 1, ps_grid%no_l_channels
1120 hato(2:) = ps_grid%rphi(2:, l, 1+is)/g%rofi(2:)
1125 dens(ip) = dens(ip) + ps%conf%occ(l, is)*hato(ip)**2/(
m_four*
m_pi)
1129 call spline_fit(g%nrval, g%rofi, hato, ps%ur(l, is), ps%projectors_sphere_threshold)
1132 call spline_fit(g%nrval, g%rofi, dens, ps%density(is), ps%projectors_sphere_threshold)
1139 hato(1:nrc) = ps_grid%KB(1:nrc, l)
1140 hato(nrc+1:g%nrval) =
m_zero
1142 call spline_fit(g%nrval, g%rofi, hato, ps%kb(l-1, 1), ps%projectors_sphere_threshold)
1147 hato(:) = ps_grid%vlocal(:)/
m_two
1148 call spline_fit(g%nrval, g%rofi, hato, ps%vl, ps%projectors_sphere_threshold)
1150 if (ps_grid%core_corrections)
then
1152 hato(2:) = ps_grid%chcore(2:)/(
m_four*
m_pi*g%r2ofi(2:))
1154 do ir = g%nrval-1, 2, -1
1155 if (hato(ir) > eps)
then
1161 hato(nrc:g%nrval) =
m_zero
1164 call spline_fit(g%nrval, g%rofi, hato, ps%core, ps%projectors_sphere_threshold)
1167 safe_deallocate_a(hato)
1168 safe_deallocate_a(dens)
1179 type(
ps_t),
intent(inout) :: ps
1180 type(
ps_xml_t),
intent(in) :: ps_xml
1183 integer :: ll, ip, is, ic, jc, ir, nrc, ii
1184 real(real64) :: rr, kbcos, kbnorm, dnrm, avgv, volume_element
1185 real(real64),
allocatable :: vlocal(:), kbprojector(:), wavefunction(:), nlcc_density(:), dens(:,:)
1186 integer,
allocatable :: cmap(:, :)
1187 real(real64),
allocatable :: matrix(:, :), eigenvalues(:)
1188 logical :: density_is_known
1189 integer,
allocatable :: order(:)
1190 real(real64),
allocatable :: occ_tmp(:,:)
1191 real(real64),
parameter :: tol_diagonal=1.0e-10_real64
1197 if (ps%file_format == pseudo_format_psp8)
then
1199 message(1) =
"SOC from PSP8 is not currently supported."
1200 message(2) =
"Only scalar relativistic effects will be considered."
1207 ps%nlcc = ps_xml%nlcc
1209 ps%z_val = ps_xml%valence_charge
1211 ps%has_density = ps_xml%has_density
1220 safe_allocate(vlocal(1:ps%g%nrval))
1221 do ip = 1, ps_xml%grid_size
1222 rr = ps_xml%grid(ip)
1223 if (ip <= ps_xml%grid_size)
then
1224 vlocal(ip) = ps_xml%potential(ip, ps%llocal)
1226 vlocal(ip) = -ps_xml%valence_charge/rr
1231 call spline_fit(ps%g%nrval, ps%g%rofi, vlocal, ps%vl, ps%projectors_sphere_threshold)
1233 safe_deallocate_a(vlocal)
1235 safe_allocate(kbprojector(1:ps%g%nrval))
1236 safe_allocate(wavefunction(1:ps%g%nrval))
1241 density_is_known = .false.
1244 if (ps_xml%kleinman_bylander)
then
1246 safe_allocate(cmap(0:ps_xml%lmax, 1:ps_xml%nchannels))
1250 do ll = 0, ps_xml%lmax
1251 do ic = 1, ps_xml%nchannels
1255 if (ll == ps_xml%llocal) cycle
1262 assert(mod(ps_xml%nchannels, 2)==0)
1264 cmap(ll, ic) = int((ic-1)/2)*2 + 2
1267 cmap(ll, ic) = int((ic-1)/2)*2 + 1
1273 assert(sum(cmap(ll, 1:ps_xml%nchannels)) == (ps_xml%nchannels + 1)*ps_xml%nchannels/2)
1276 assert(all(cmap >= 0 .and. cmap <= ps_xml%nchannels))
1286 safe_allocate(matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels))
1287 safe_allocate(eigenvalues(1:ps_xml%nchannels))
1289 do ll = 0, ps_xml%lmax
1291 if (
is_diagonal(ps_xml%nchannels, ps_xml%dij(ll, :, :), tol_diagonal) .or. &
1294 do ic = 1, ps_xml%nchannels
1295 eigenvalues(ic) = ps_xml%dij(ll, ic, ic)
1296 matrix(ic, ic) =
m_one
1300 matrix(1:ps_xml%nchannels, 1:ps_xml%nchannels) = ps_xml%dij(ll, 1:ps_xml%nchannels, 1:ps_xml%nchannels)
1304 do ic = 1, ps_xml%nchannels
1306 do jc = 1, ps_xml%nchannels
1307 call lalg_axpy(ps_xml%grid_size, matrix(jc, ic), ps_xml%projector(:, ll, jc), kbprojector)
1310 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, cmap(ll, ic)), ps%projectors_sphere_threshold)
1312 ps%h(ll, cmap(ll, ic), cmap(ll, ic)) = eigenvalues(ic)
1317 safe_deallocate_a(matrix)
1318 safe_deallocate_a(eigenvalues)
1321 ps%conf%p = ps_xml%nwavefunctions
1328 if ((.not.
ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0)
then
1329 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1333 do ii = 1, ps_xml%nwavefunctions
1335 ps%conf%n(ii) = ps_xml%wf_n(ii)
1336 ps%conf%l(ii) = ps_xml%wf_l(ii)
1338 if (ps%ispin == 2)
then
1339 ps%conf%occ(ii, 1) = min(ps_xml%wf_occ(ii),
m_two*ps_xml%wf_l(ii) +
m_one)
1340 ps%conf%occ(ii, 2) = ps_xml%wf_occ(ii) - ps%conf%occ(ii, 1)
1342 ps%conf%occ(ii, 1) = ps_xml%wf_occ(ii)
1350 do ip = 1, ps%g%nrval
1351 if (ip <= ps_xml%grid_size)
then
1352 wavefunction(ip) = ps_xml%wavefunction(ip, ii)
1354 wavefunction(ip) =
m_zero
1359 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ii, is), ps%projectors_sphere_threshold)
1364 do ip = 1, ps_xml%grid_size
1365 dens(ip, is) = dens(ip, is) + ps%conf%occ(ii, is)*wavefunction(ip)**2/(
m_four*
m_pi)
1373 if ((.not.
ps_has_density(ps) .or. ps%ispin == 2) .and. ps_xml%nwavefunctions > 0)
then
1375 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1377 safe_deallocate_a(dens)
1378 density_is_known = .
true.
1379 ps%has_density = .
true.
1382 safe_deallocate_a(cmap)
1388 ps%conf%symbol = ps%label(1:2)
1392 safe_allocate(order(1:ps_xml%lmax+1))
1393 safe_allocate(occ_tmp(1:ps_xml%lmax+1, 1:2))
1394 occ_tmp(:,:) = ps%conf%occ(1:ps_xml%lmax+1,1:2)
1395 call sort(ps%conf%l(1:ps_xml%lmax+1), order)
1396 do ll = 0, ps_xml%lmax
1397 ps%conf%occ(ll+1, 1:2) = occ_tmp(order(ll+1), 1:2)
1399 safe_deallocate_a(order)
1400 safe_deallocate_a(occ_tmp)
1403 if (abs(sum(ps%conf%occ) - ps%z_val ) <
m_epsilon)
then
1404 safe_allocate(dens(1:ps%g%nrval, 1:ps%ispin))
1408 do ll = 0, ps_xml%lmax
1413 do ip = 1, ps_xml%grid_size
1414 rr = ps_xml%grid(ip)
1415 volume_element = rr**2*ps_xml%weights(ip)
1416 kbprojector(ip) = (ps_xml%potential(ip, ll) - ps_xml%potential(ip, ps%llocal))*ps_xml%wavefunction(ip, ll)
1417 dnrm = dnrm + kbprojector(ip)**2*volume_element
1418 avgv = avgv + kbprojector(ip)*ps_xml%wavefunction(ip, ll)*volume_element
1421 kbcos = dnrm/(safe_tol(avgv,1.0e-20_real64))
1422 kbnorm =
m_one/(safe_tol(
sqrt(dnrm),1.0e-20_real64))
1424 if (ll /= ps%llocal)
then
1425 ps%h(ll, 1, 1) = kbcos
1426 kbprojector = kbprojector*kbnorm
1431 call spline_fit(ps%g%nrval, ps%g%rofi, kbprojector, ps%kb(ll, 1), ps%projectors_sphere_threshold)
1434 call lalg_copy(ps_xml%grid_size, ps_xml%wavefunction(:, ll), wavefunction)
1435 wavefunction(ps_xml%grid_size+1:ps%g%nrval) =
m_zero
1438 call spline_fit(ps%g%nrval, ps%g%rofi, wavefunction, ps%ur(ll + 1, is), ps%projectors_sphere_threshold)
1442 if (abs(sum(ps%conf%occ) - ps%z_val) <
m_epsilon)
then
1444 do ip = 1, ps_xml%grid_size
1445 dens(ip, is) = dens(ip, is) + ps%conf%occ(ll+1, is)*wavefunction(ip)**2/(
m_four*
m_pi)
1452 if (abs(sum(ps%conf%occ) - ps%z_val) <
m_epsilon)
then
1454 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:,is), ps%density(is), ps%projectors_sphere_threshold)
1456 safe_deallocate_a(dens)
1457 density_is_known = .
true.
1458 ps%has_density = .
true.
1465 safe_allocate(dens(1:ps%g%nrval, 1))
1467 dens(1:ps_xml%grid_size, 1) = ps_xml%density(1:ps_xml%grid_size)/ps%ispin
1468 dens(ps_xml%grid_size + 1:ps%g%nrval, 1) =
m_zero
1471 call spline_fit(ps%g%nrval, ps%g%rofi, dens(:, 1), ps%density(is), ps%projectors_sphere_threshold)
1474 safe_deallocate_a(dens)
1479 if (ps_xml%nlcc)
then
1481 safe_allocate(nlcc_density(1:ps%g%nrval))
1483 call lalg_copy(ps_xml%grid_size, ps_xml%nlcc_density, nlcc_density)
1486 do ir = ps_xml%grid_size - 1, 1, -1
1487 if (nlcc_density(ir) > eps)
then
1493 nlcc_density(nrc:ps%g%nrval) =
m_zero
1495 call spline_fit(ps%g%nrval, ps%g%rofi, nlcc_density, ps%core, ps%projectors_sphere_threshold)
1497 safe_deallocate_a(nlcc_density)
1503 ps%rc_max = ps%rc_max*1.05_real64
1505 safe_deallocate_a(kbprojector)
1506 safe_deallocate_a(wavefunction)
1515 type(
ps_t),
intent(in) :: ps
1530 type(
ps_t),
intent(in) :: ps
1537 if (any(.not. ps%bound(i,:))) cycle
1546 type(
ps_t),
intent(in) :: ps
1548 has_density = ps%has_density
1554 pure logical function ps_has_nlcc(ps)
result(has_nlcc)
1555 type(
ps_t),
intent(in) :: ps
1562 real(real64) function ps_density_volume(ps, namespace) result(volume)
1563 type(
ps_t),
intent(in) :: ps
1566 integer :: ip, ispin
1568 real(real64),
allocatable ::vol(:)
1571 push_sub(ps_density_volume)
1574 message(1) =
"The pseudopotential does not contain an atomic density"
1578 safe_allocate(vol(1:ps%g%nrval))
1580 do ip = 1, ps%g%nrval
1583 do ispin = 1, ps%ispin
1589 call spline_fit(ps%g%nrval, ps%g%rofi, vol, volspl, ps%projectors_sphere_threshold)
1593 safe_deallocate_a(vol)
1595 pop_sub(ps_density_volume)
1602 type(namespace_t),
intent(in) :: namespace
1603 real(real64),
intent(in) :: zz
1604 real(real64),
intent(in) :: valcharge
1605 integer,
intent(in) :: ispin
1606 type(valconf_t),
intent(inout) :: conf
1618 write(message(1),
'(a,a,a)')
'Debug: Guessing the atomic occupations for ', trim(conf%symbol),
"."
1619 call messages_info(1, namespace=namespace, debug_only=.
true.)
1621 assert(valcharge <= zz)
1625 if(int(zz) > 2 .and. val > zz - 2)
then
1629 if(int(zz) > 4 .and. val > zz - 4)
then
1634 if(int(zz) > 18 .and. val > zz - 10)
then
1638 if(int(zz) > 12 .and. val > zz - 12)
then
1642 if(int(zz) > 18 .and. val > zz - 18)
then
1646 if(int(zz) > 28 .and. val > zz - 28)
then
1650 if(int(zz) > 30 .and. val > zz - 30)
then
1654 if(int(zz) > 36 .and. val > zz - 36)
then
1658 if(int(zz) > 46 .and. val > zz - 46)
then
1663 if((int(zz) > 48 .and. val > zz - 48) .or. &
1664 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 62))
then
1669 if((int(zz) > 54 .and. val > zz - 54) .or. &
1670 (int(zz) > 70 .and. int(zz) <= 81 .and. val > zz - 68) )
then
1675 if(int(zz) > 80 .and. val > zz - 68)
then
1681 select case (int(zz))
1998 if (val < m_epsilon)
then
2000 if (ispin == 2)
then
2001 call valconf_unpolarized_to_polarized(conf)
2005 message(1) =
"Error in attributing atomic occupations"
2006 call messages_warning(1, namespace=namespace)
2013 real(real64),
intent(inout) :: val
2014 integer,
intent(in) :: max_occ, nn
2017 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2018 val = val - conf%occ(conf%p,1)
2024 real(real64),
intent(inout) :: val
2025 integer,
intent(in) :: max_occ, nn
2028 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2029 val = val - conf%occ(conf%p,1)
2035 real(real64),
intent(inout) :: val
2036 integer,
intent(in) :: max_occ, nn
2039 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2040 val = val - conf%occ(conf%p,1)
2046 real(real64),
intent(inout) :: val
2047 integer,
intent(in) :: max_occ, nn
2050 conf%occ(conf%p,1) = min(val, real(max_occ, real64) )
2051 val = val - conf%occ(conf%p,1)
2058#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__
subroutine, public valconf_unpolarized_to_polarized(conf)
subroutine, public valconf_copy(cout, cin)
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)
character(kind=c_char, len=1) function, dimension(:), allocatable, public string_f_to_c(f_string)
convert a Fortran string to a C string
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