26 use,
intrinsic :: iso_fortran_env
56 integer,
public :: method
57 real(real64),
public :: dsmear
58 real(real64) :: dsmear_cond
59 real(real64),
public :: e_fermi
60 real(real64) :: e_fermi_cond
62 integer,
public :: el_per_state
63 real(real64),
public :: ef_occ
64 real(real64) :: ef_occ_cond
65 logical,
public :: integral_occs
67 integer :: fermi_count
72 logical,
public:: photodop
73 real(real64) :: nephotodop
74 integer :: photodop_bandmin
78 integer,
parameter,
public :: &
79 SMEAR_SEMICONDUCTOR = 1, &
86 real(real64),
parameter :: TOL_SMEAR = 1e-6_real64
91 subroutine smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
92 type(smear_t),
intent(out) :: this
93 type(namespace_t),
intent(in) :: namespace
94 integer,
intent(in) :: ispin
95 logical,
intent(in) :: fixed_occ
96 logical,
intent(in) :: integral_occs
97 type(kpoints_t),
intent(in) :: kpoints
101 this%integral_occs = integral_occs
129 call parse_variable(namespace,
'SmearingFunction', smear_semiconductor, this%method)
143 this%dsmear = 1e-14_real64
144 if (this%method /= smear_semiconductor .and. this%method /=
smear_fixed_occ)
then
157 this%dsmear_cond = this%dsmear
175 this%photodop=.false.
186 call parse_variable(namespace,
'PhotoDopingBand', 1, this%photodop_bandmin)
192 this%el_per_state = 2
194 this%el_per_state = 1
203 if (this%method == smear_semiconductor)
then
206 if (this%nik_factor == 0)
then
207 message(1) =
"k-point weights in KPoints or KPointsReduced blocks must be rational numbers for semiconducting smearing."
230 type(
smear_t),
intent(out) :: to
231 type(
smear_t),
intent(in) :: from
235 to%method = from%method
236 to%dsmear = from%dsmear
237 to%e_fermi = from%e_fermi
238 to%el_per_state = from%el_per_state
239 to%ef_occ = from%ef_occ
241 to%fermi_count = from%fermi_count
242 to%nik_factor = from%nik_factor
250 qtot, nik, nst, kweights)
251 type(
smear_t),
intent(inout) :: this
253 real(real64),
intent(in) :: eigenvalues(:,:), occupations(:,:)
254 real(real64),
intent(in) :: qtot, kweights(:)
255 integer,
intent(in) :: nik, nst
257 real(real64),
parameter :: tol = 1.0e-10_real64
258 integer :: ist, ik, iter, maxq, weight, sumq_int, sum_weight
259 real(real64) :: sumq_frac
261 real(real64),
allocatable :: eigenval_list(:)
262 integer,
allocatable :: k_list(:), reorder(:)
263 integer :: fermi_count_up, fermi_count_down
267 maxq = this%el_per_state * nst * this%nspins
268 if (maxq - qtot <= -tol)
then
269 message(1) =
'Not enough states'
270 write(
message(2),
'(6x,a,f12.6,a,i10)')
'(total charge = ', qtot, &
271 ' max charge = ', maxq
277 ist_cycle:
do ist = nst, 1, -1
278 if (any(occupations(ist, :) >
m_min_occ))
then
279 this%e_fermi = eigenvalues(ist, 1)
280 this%ef_occ = occupations(ist, 1) / this%el_per_state
282 if (eigenvalues(ist, ik) > this%e_fermi .and. occupations(ist, ik) >
m_min_occ)
then
283 this%e_fermi = eigenvalues(ist, ik)
284 this%ef_occ = occupations(ist, ik) / this%el_per_state
291 else if (this%method == smear_semiconductor)
then
293 safe_allocate(eigenval_list(1:nst * nik))
294 safe_allocate( k_list(1:nst * nik))
295 safe_allocate( reorder(1:nst * nik))
300 eigenval_list(iter) = eigenvalues(ist, ik)
307 call sort(eigenval_list, reorder)
309 sumq_int =
floor(qtot) * this%nik_factor
310 sumq_frac = qtot * this%nik_factor - sumq_int
311 if ( sumq_frac + tol > this%nik_factor )
then
312 sumq_int = sumq_int + this%nik_factor
313 sumq_frac = sumq_frac - this%nik_factor
315 if (sumq_frac < tol) sumq_frac =
m_zero
317 do iter = 1, nst * nik
318 weight = int(kweights(k_list(reorder(iter))) * this%nik_factor +
m_half)
319 if (.not. weight > 0) cycle
320 this%e_fermi = eigenval_list(iter)
321 this%ef_occ = (sumq_int + sumq_frac) / (weight * this%el_per_state)
323 if (sumq_int - weight * this%el_per_state <= 0)
then
331 if (iter - fermi_count_down < 1)
exit
332 if (abs(this%e_fermi - eigenval_list(iter - fermi_count_down)) > tol_smear)
exit
333 weight = int(kweights(k_list(reorder(iter - fermi_count_down))) * this%nik_factor +
m_half)
335 sumq_int = sumq_int + weight * this%el_per_state
336 sum_weight = sum_weight + weight
338 fermi_count_down = fermi_count_down + 1
341 if (iter + fermi_count_up > nst*nik)
exit
342 if (abs(this%e_fermi - eigenval_list(iter + fermi_count_up)) > tol_smear)
exit
343 weight = int(kweights(k_list(reorder(iter + fermi_count_up))) * this%nik_factor +
m_half)
345 sum_weight = sum_weight + weight
347 fermi_count_up = fermi_count_up + 1
349 this%fermi_count = fermi_count_up + fermi_count_down - 1
350 this%ef_occ = (sumq_int + sumq_frac) / (sum_weight * this%el_per_state)
354 sumq_int = sumq_int - weight * this%el_per_state
358 safe_deallocate_a(eigenval_list)
359 safe_deallocate_a(k_list)
360 safe_deallocate_a(reorder)
363 if (this%photodop)
then
366 eigenvalues, kweights, nik, qtot-this%nephotodop, 1, this%photodop_bandmin-1, &
367 this%e_fermi, this%ef_occ)
371 eigenvalues, kweights, nik, this%nephotodop, this%photodop_bandmin, nst, &
372 this%e_fermi_cond, this%ef_occ_cond)
375 eigenvalues, kweights, nik, qtot, 1, nst, this%e_fermi, this%ef_occ)
384 eigenvalues, kweights, nik, q_in, start_band, end_band, e_fermi, ef_occ)
385 type(
smear_t),
intent(inout) :: this
387 real(real64),
intent(in) :: dsmear_in
388 real(real64),
intent(in) :: tol
389 integer,
intent(in) :: nik
390 real(real64),
intent(in) :: eigenvalues(:,:)
391 real(real64),
intent(in) :: kweights(:), q_in
392 integer,
intent(in) :: start_band, end_band
393 real(real64),
intent(out) :: e_fermi, ef_occ
395 integer,
parameter :: nitmax = 200
396 integer :: ist, ik, iter
397 real(real64) :: drange, dsmear, emin, emax, xx, sumq
402 dsmear = max(1e-14_real64, dsmear_in)
403 drange = dsmear *
sqrt(-
log(tol * 0.01_real64))
405 emin = minval(eigenvalues) - drange
406 emax = maxval(eigenvalues) + drange
409 e_fermi =
m_half * (emin + emax)
413 do ist = start_band, end_band
414 xx = (e_fermi - eigenvalues(ist, ik))/dsmear
415 sumq = sumq + kweights(ik) * this%el_per_state * &
420 conv = (abs(sumq - q_in) <= tol)
423 if (sumq <= q_in) emin = e_fermi
424 if (sumq >= q_in) emax = e_fermi
430 message(1) =
'Fermi: did not converge.'
439 type(
smear_t),
intent(in) :: this
440 real(real64),
intent(in) :: eigenvalues(:,:)
441 real(real64),
intent(inout) :: occupations(:,:)
442 integer,
intent(in) :: nik, nst
444 integer :: ik, ist, ifermi
445 real(real64) :: dsmear, xx, dsmear_cond
451 else if (this%method == smear_semiconductor)
then
452 assert(this%fermi_count > 0 .and. this%fermi_count <= nik*nst)
457 xx = eigenvalues(ist, ik) - this%e_fermi
458 if (xx < -tol_smear)
then
459 occupations(ist, ik) = this%el_per_state
460 else if (abs(xx) <= tol_smear .and. ifermi < this%fermi_count)
then
461 occupations(ist, ik) = this%ef_occ * this%el_per_state
464 occupations(ist, ik) =
m_zero
471 dsmear = max(1e-14_real64, this%dsmear)
472 if (this%photodop)
then
473 dsmear_cond = max(1e-14_real64, this%dsmear_cond)
476 if (ist < this%photodop_bandmin)
then
478 xx = (this%e_fermi - eigenvalues(ist, ik))/dsmear
481 xx = (this%e_fermi_cond - eigenvalues(ist, ik))/dsmear_cond
489 xx = (this%e_fermi - eigenvalues(ist, ik))/dsmear
501 real(real64) function smear_calc_entropy(this, eigenvalues, &
502 nik, nst, kweights, occ) result(entropy)
503 type(
smear_t),
intent(inout) :: this
504 real(real64),
intent(in) :: eigenvalues(:,:)
505 real(real64),
intent(in) :: kweights(:)
506 integer,
intent(in) :: nik, nst
507 real(real64),
intent(in) :: occ(:, :)
510 real(real64) :: dsmear, xx, term, ff
512 push_sub(smear_calc_entropy)
516 if (this%method ==
smear_fixed_occ .or. this%method == smear_semiconductor)
then
522 ff = occ(ist, ik) / this%el_per_state
524 term = ff *
log(ff) + (1 - ff) *
log(1 - ff)
528 entropy = entropy - kweights(ik) * this%el_per_state * term
532 dsmear = max(1e-14_real64, this%dsmear)
536 if (eigenvalues(ist, ik) < huge(
m_one))
then
537 xx = (this%e_fermi - eigenvalues(ist, ik)) / dsmear
538 entropy = entropy - kweights(ik) * this%el_per_state * &
545 pop_sub(smear_calc_entropy)
551 type(
smear_t),
intent(in) :: this
552 real(real64),
intent(in) :: xx
554 real(real64),
parameter :: maxarg = 200.0_real64
555 real(real64) :: xp, arg, hd, hp, aa
564 select case (this%method)
566 if (abs(xx) <= m_epsilon)
then
571 if (abs(xx) <= 36.0_real64)
then
572 deltaf = m_one / (m_two +
exp(-xx) +
exp(xx))
576 xp = xx - m_one /
sqrt(m_two)
577 arg = min(maxarg, xp**2)
579 deltaf =
exp(-arg) /
sqrt(m_pi) * (m_two -
sqrt(m_two) * xx)
582 arg = min(maxarg, xx**2)
583 deltaf =
exp(-arg) /
sqrt(m_pi)
585 if (this%MP_n > 0)
then
589 aa = m_one /
sqrt(m_pi)
591 hd = m_two * xx * hp - m_two * ni * hd
593 aa = -aa / (m_four * ii)
594 hp = m_two * xx * hd - m_two * ni * hp
596 deltaf = deltaf + aa * hp
601 xp = abs(xx) + m_one /
sqrt(m_two)
602 deltaf =
sqrt(m_e) * xp *
exp(-xp * xp)
611 type(
smear_t),
intent(in) :: this
612 real(real64),
intent(in) :: xx
614 real(real64),
parameter :: maxarg = 200.0_real64
615 real(real64) :: xp, arg, hd, hp, aa
624 select case (this%method)
626 if (xx > m_zero)
then
628 else if (abs(xx) <= m_epsilon)
then
633 if (xx > maxarg)
then
635 else if (xx > -maxarg)
then
636 stepf = m_one / (m_one +
exp(-xx))
640 xp = xx - m_one /
sqrt(m_two)
641 arg = min(maxarg, xp**2)
643 stepf = m_half * loct_erf(xp) + &
644 m_one /
sqrt(m_two * m_pi) *
exp(-arg) + m_half
647 stepf = m_half * loct_erfc(-xx)
649 if (this%MP_n > 0)
then
651 arg = min(maxarg, xx**2)
654 aa = m_one /
sqrt(m_pi)
656 hd = m_two * xx * hp - m_two * ni * hd
658 aa = -aa / (m_four * ii)
659 stepf = stepf - aa * hd
660 hp = m_two * xx * hd - m_two * ni * hp
666 if (xx <= m_zero)
then
667 xp = xx - m_one /
sqrt(m_two)
668 stepf = m_half *
sqrt(m_e) *
exp(-xp * xp)
670 xp = xx + m_one /
sqrt(m_two)
671 stepf = m_one - m_half *
sqrt(m_e) *
exp(-xp * xp)
683 type(
smear_t),
intent(in) :: this
684 real(real64),
intent(in) :: xx
686 real(real64),
parameter :: maxarg = 200.0_real64
687 real(real64) :: xp, arg, hd, hp, hpm1, aa
696 select case (this%method)
700 if (abs(xx) <= 36.0_real64)
then
701 xp = m_one / (m_one +
exp(-xx))
702 entropyf = xp *
log(xp) + (m_one - xp) *
log(m_one - xp)
706 xp = xx - m_one /
sqrt(m_two)
707 arg = min(maxarg, xp**2)
709 entropyf = m_one /
sqrt(m_two * m_pi) * xp *
exp(-arg)
712 arg = min(maxarg, xx**2)
713 entropyf = -m_half *
exp(-arg) /
sqrt(m_pi)
715 if (this%MP_n > 0)
then
719 aa = m_one /
sqrt(m_pi)
721 hd = m_two * xx * hp - m_two * ni * hd
724 hp = m_two * xx * hd - m_two * ni * hp
726 aa = -aa / (m_four * ii)
727 entropyf = entropyf - aa * (m_half * hp + hpm1 * ni)
732 xp = abs(xx) + m_one /
sqrt(m_two)
733 entropyf = -
sqrt(m_e) * (abs(xx) *
exp(-xp * xp) / m_two +
sqrt(m_pi) / m_four * loct_erfc(xp))
743 type(
smear_t),
intent(in) :: this
750 type(
smear_t),
intent(in) :: this
751 type(namespace_t),
intent(in) :: namespace
752 integer,
optional,
intent(in) :: iunit
756 if (this%photodop)
then
757 write(message(1),
'(a,f12.6,1x,a)')
"Fermi energy (valence ) = ", &
758 units_from_atomic(units_out%energy, this%e_fermi), units_abbrev(units_out%energy)
759 write(message(2),
'(a,f12.6,1x,a)')
"Fermi energy (conduction) = ", &
760 units_from_atomic(units_out%energy, this%e_fermi_cond), units_abbrev(units_out%energy)
761 call messages_info(2, iunit=iunit, namespace=namespace)
763 write(message(1),
'(a,f12.6,1x,a)')
"Fermi energy = ", &
764 units_from_atomic(units_out%energy, this%e_fermi), units_abbrev(units_out%energy)
765 call messages_info(1, iunit=iunit, namespace=namespace)
This is the common interface to a sorting routine. It performs the shell algorithm,...
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double floor(double __x) __attribute__((__nothrow__
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public p_ry
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public m_one
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
integer function, public kpoints_kweight_denominator(this)
subroutine, public messages_obsolete_variable(namespace, name, rep)
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)
real(real64) function, public smear_calc_entropy(this, eigenvalues, nik, nst, kweights, occ)
subroutine, public smear_fill_occupations(this, eigenvalues, occupations, nik, nst)
subroutine, public smear_find_fermi_energy(this, namespace, eigenvalues, occupations, qtot, nik, nst, kweights)
real(real64) function, public smear_entropy_function(this, xx)
This function is defined as .
real(real64) function, public smear_delta_function(this, xx)
integer, parameter, public smear_fermi_dirac
integer, parameter, public smear_methfessel_paxton
subroutine, public smear_write_info(this, namespace, iunit)
real(real64) function, public smear_step_function(this, xx)
integer, parameter, public smear_semiconductor
subroutine, public smear_copy(to, from)
integer, parameter, public smear_fixed_occ
integer, parameter, public smear_spline
subroutine bisection_find_fermi_energy(this, namespace, dsmear_in, tol, eigenvalues, kweights, nik, q_in, start_band, end_band, e_fermi, ef_occ)
subroutine, public smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
integer, parameter, public smear_cold
logical pure function, public smear_is_semiconducting(this)
This module is intended to contain "only mathematical" functions and procedures.
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing