Octopus
smear.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
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! Some pieces glued from PW/w(0,1)gauss.f90 from PWSCF
20
21#include "global.h"
22
23module smear_oct_m
24 use debug_oct_m
25 use global_oct_m
26 use, intrinsic :: iso_fortran_env
31 use parser_oct_m
33 use sort_oct_m
34 use unit_oct_m
37
38 implicit none
39
40 private
41 public :: &
42 smear_t, &
43 smear_init, &
44 smear_copy, &
53
54 type smear_t
55 private
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
61
62 integer, public :: el_per_state
63 real(real64), public :: ef_occ
64 real(real64) :: ef_occ_cond
65 logical, public :: integral_occs
66 integer :: MP_n
67 integer :: fermi_count
68 integer :: nik_factor
69 integer :: nspins
70
71 ! PhotonDoping
72 logical, public:: photodop
73 real(real64) :: nephotodop
74 integer :: photodop_bandmin
75
76 end type smear_t
77
78 integer, parameter, public :: &
79 SMEAR_SEMICONDUCTOR = 1, &
81 smear_cold = 3, &
83 smear_spline = 5, &
85
86 real(real64), parameter :: TOL_SMEAR = 1e-6_real64
87
88contains
89
90 !--------------------------------------------------
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
98
99 push_sub(smear_init)
100
101 this%integral_occs = integral_occs
102
103 !%Variable SmearingFunction
104 !%Type integer
105 !%Default semiconducting
106 !%Section States
107 !%Description
108 !% This is the function used to smear the electronic occupations.
109 !% It is ignored if the <tt>Occupations</tt> block is set.
110 !%Option semiconducting 1
111 !% Semiconducting occupations, <i>i.e.</i> the lowest lying states are occupied
112 !% until no more electrons are left.
113 !%Option fermi_dirac 2
114 !% Simple Fermi-Dirac distribution. In this case, <tt>Smearing</tt> has
115 !% the meaning of an electronic temperature. DN Mermin, <i>Phys. Rev.</i> <b>137</b>, A1441 (1965).
116 !%Option cold_smearing 3
117 !% N Marzari, D Vanderbilt, A De Vita, and MC Payne, <i>Phys. Rev. Lett.</i> <b>82</b>, 3296 (1999).
118 !%Option methfessel_paxton 4
119 !% M Methfessel and AT Paxton, <i>Phys. Rev. B</i> <b>40</b>, 3616 (1989).
120 !% In this case, the variable <tt>SmearingMPOrder</tt> sets the order of the smearing.
121 !% Occupations may be negative.
122 !%Option spline_smearing 5
123 !% Nearly identical to Gaussian smearing.
124 !% JM Holender, MJ Gillan, MC Payne, and AD Simpson, <i>Phys. Rev. B</i> <b>52</b>, 967 (1995).
125 !%End
126 if (fixed_occ) then
127 this%method = smear_fixed_occ
128 else
129 call parse_variable(namespace, 'SmearingFunction', smear_semiconductor, this%method)
130 if (.not. varinfo_valid_option('SmearingFunction', this%method)) call messages_input_error(namespace, 'SmearingFunction')
131 call messages_print_var_option('SmearingFunction', this%method, namespace=namespace)
132 end if
133
134 !%Variable Smearing
135 !%Type float
136 !%Default 0.1 eV
137 !%Section States
138 !%Description
139 !% If <tt>Occupations</tt> is not set, <tt>Smearing</tt> is the
140 !% smearing width used in the <tt>SmearingFunction</tt> to distribute the electrons
141 !% among the existing states.
142 !%End
143 this%dsmear = 1e-14_real64
144 if (this%method /= smear_semiconductor .and. this%method /= smear_fixed_occ) then
145 call parse_variable(namespace, 'Smearing', 0.1_real64 / (m_two * p_ry), this%dsmear, units_inp%energy)
146 end if
148 !%Variable PhotoDopingSmearing
149 !%Type float
150 !%Default 0.1 eV
151 !%Section States
152 !%Description
153 !% If <tt>Occupations</tt> is not set, <tt>PhotoDopingSmearing</tt> is the
154 !% smearing width used in the <tt>SmearingFunction</tt> to distribute the electrons in conduction bands
155 !% among the existing states.
156 !%End
157 this%dsmear_cond = this%dsmear
158 if (this%method /= smear_semiconductor .and. this%method /= smear_fixed_occ) then
159 call parse_variable(namespace, 'PhotoDopingSmearing', 0.1_real64 / (m_two * p_ry), this%dsmear_cond, units_inp%energy)
160 end if
162 !%Variable PhotoDopingNumElectrons
163 !%Type float
164 !%Default 0
165 !%Section States
166 !%Description
167 !% This variable allows to set the number of valence electrons taken out to put
168 !% into the conduction bands. The method follows Marini and Calandra, PRB 104, 144103 (2021).
169 !%End
170 call parse_variable(namespace, 'PhotoDopingNumElectrons', m_zero, this%nephotodop)
172 if (this%nephotodop > m_min_occ) then
173 this%photodop=.true.
174 else
175 this%photodop=.false.
176 endif
177
178
179 !%Variable PhotoDopingBand
180 !%Type integer
181 !%Default 1
182 !%Section States
183 !%Description
184 !% This variable specifies where the minium (starting) conduction band index for the photodoping electrons.
185 !%End
186 call parse_variable(namespace, 'PhotoDopingBand', 1, this%photodop_bandmin)
187
188
189 call messages_obsolete_variable(namespace, 'ElectronicTemperature', 'Smearing')
190
191 if (ispin == 1) then ! unpolarized
192 this%el_per_state = 2
193 else
194 this%el_per_state = 1
195 end if
196
197 if (ispin == 2) then
198 this%nspins = 2
199 else
200 this%nspins = 1
201 end if
202
203 if (this%method == smear_semiconductor) then
204 this%nik_factor = kpoints_kweight_denominator(kpoints)
205
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."
208 call messages_fatal(1, namespace=namespace)
209 end if
210 end if
211
212 this%MP_n = 0
213 if (this%method == smear_methfessel_paxton) then
214 !%Variable SmearingMPOrder
215 !%Type integer
216 !%Default 1
217 !%Section States
218 !%Description
219 !% Sets the order of the Methfessel-Paxton smearing function.
220 !%End
221 call parse_variable(namespace, 'SmearingMPOrder', 1, this%MP_n)
222 end if
223
224 pop_sub(smear_init)
225 end subroutine smear_init
226
227
228 !--------------------------------------------------
229 subroutine smear_copy(to, from)
230 type(smear_t), intent(out) :: to
231 type(smear_t), intent(in) :: from
232
233 push_sub(smear_copy)
234
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
240 to%MP_n = from%MP_n
241 to%fermi_count = from%fermi_count
242 to%nik_factor = from%nik_factor
243
244 pop_sub(smear_copy)
245 end subroutine smear_copy
246
247
248 !--------------------------------------------------
249 subroutine smear_find_fermi_energy(this, namespace, eigenvalues, occupations, &
250 qtot, nik, nst, kweights)
251 type(smear_t), intent(inout) :: this
252 type(namespace_t), intent(in) :: namespace
253 real(real64), intent(in) :: eigenvalues(:,:), occupations(:,:)
254 real(real64), intent(in) :: qtot, kweights(:)
255 integer, intent(in) :: nik, nst
256
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
260 logical :: conv
261 real(real64), allocatable :: eigenval_list(:)
262 integer, allocatable :: k_list(:), reorder(:)
263 integer :: fermi_count_up, fermi_count_down
264
266
267 maxq = this%el_per_state * nst * this%nspins
268 if (maxq - qtot <= -tol) then ! not enough states
269 message(1) = 'Not enough states'
270 write(message(2),'(6x,a,f12.6,a,i10)')'(total charge = ', qtot, &
271 ' max charge = ', maxq
272 call messages_fatal(2, namespace=namespace)
273 end if
274
275 conv = .true.
276 if (this%method == smear_fixed_occ) then ! Fermi energy: last of occupied states
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
281 do ik = 2, nik
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
285 end if
286 end do
287 exit ist_cycle
288 end if
289 end do ist_cycle
290
291 else if (this%method == smear_semiconductor) then
292 ! first we sort the eigenvalues
293 safe_allocate(eigenval_list(1:nst * nik))
294 safe_allocate( k_list(1:nst * nik))
295 safe_allocate( reorder(1:nst * nik))
296
297 iter = 1
298 do ist = 1, nst
299 do ik = 1, nik
300 eigenval_list(iter) = eigenvalues(ist, ik)
301 k_list(iter) = ik
302 reorder(iter) = iter
303 iter = iter + 1
304 end do
305 end do
306
307 call sort(eigenval_list, reorder)
308
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
314 end if
315 if (sumq_frac < tol) sumq_frac = m_zero
316
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
324 ! count how many occupied states are at the fermi level,
325 ! this is required later to fill the states
326 this%fermi_count = 1
327 fermi_count_up = 1
328 fermi_count_down = 1
329 sum_weight = weight
330 do
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)
334 if (weight > m_epsilon) then
335 sumq_int = sumq_int + weight * this%el_per_state
336 sum_weight = sum_weight + weight
337 end if
338 fermi_count_down = fermi_count_down + 1
339 end do
340 do
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)
344 if (weight > m_epsilon) then
345 sum_weight = sum_weight + weight
346 end if
347 fermi_count_up = fermi_count_up + 1
348 end do
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)
351 exit
352 end if
353
354 sumq_int = sumq_int - weight * this%el_per_state
355 end do
356 assert(this%ef_occ < m_one + 10.0_real64*m_epsilon)
357
358 safe_deallocate_a(eigenval_list)
359 safe_deallocate_a(k_list)
360 safe_deallocate_a(reorder)
361
362 else ! bisection
363 if (this%photodop) then
364 ! find fermi energy for valence electrons
365 call bisection_find_fermi_energy(this, namespace, this%dsmear, tol, &
366 eigenvalues, kweights, nik, qtot-this%nephotodop, 1, this%photodop_bandmin-1, &
367 this%e_fermi, this%ef_occ)
368
369 ! find fermi energy for conduction electrons
370 call bisection_find_fermi_energy(this, namespace, this%dsmear_cond, tol, &
371 eigenvalues, kweights, nik, this%nephotodop, this%photodop_bandmin, nst, &
372 this%e_fermi_cond, this%ef_occ_cond)
373 else
374 call bisection_find_fermi_energy(this, namespace, this%dsmear, tol, &
375 eigenvalues, kweights, nik, qtot, 1, nst, this%e_fermi, this%ef_occ)
376 end if
377 end if
378
380 end subroutine smear_find_fermi_energy
381
382 ! ---------------------------------------------------------
383 subroutine bisection_find_fermi_energy(this, namespace, dsmear_in, tol, &
384 eigenvalues, kweights, nik, q_in, start_band, end_band, e_fermi, ef_occ)
385 type(smear_t), intent(inout) :: this
386 type(namespace_t), intent(in) :: namespace
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
394
395 integer, parameter :: nitmax = 200
396 integer :: ist, ik, iter
397 real(real64) :: drange, dsmear, emin, emax, xx, sumq
398 logical :: conv
399
401
402 dsmear = max(1e-14_real64, dsmear_in)
403 drange = dsmear * sqrt(-log(tol * 0.01_real64))
404
405 emin = minval(eigenvalues) - drange
406 emax = maxval(eigenvalues) + drange
407
408 do iter = 1, nitmax
409 e_fermi = m_half * (emin + emax)
410 sumq = m_zero
411
412 do ik = 1, nik
413 do ist = start_band, end_band
414 xx = (e_fermi - eigenvalues(ist, ik))/dsmear
415 sumq = sumq + kweights(ik) * this%el_per_state * &
416 smear_step_function(this, xx)
417 end do
418 end do
419
420 conv = (abs(sumq - q_in) <= tol)
421 if (conv) exit
422
423 if (sumq <= q_in) emin = e_fermi
424 if (sumq >= q_in) emax = e_fermi
425
426 ef_occ = smear_step_function(this, m_zero)
427 end do
428
429 if (.not. conv) then
430 message(1) = 'Fermi: did not converge.'
431 call messages_fatal(1, namespace=namespace)
432 end if
433
435 end subroutine bisection_find_fermi_energy
436
437 ! ---------------------------------------------------------
438 subroutine smear_fill_occupations(this, eigenvalues, occupations, nik, nst)
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
443
444 integer :: ik, ist, ifermi
445 real(real64) :: dsmear, xx, dsmear_cond
446
447 push_sub(smear_fill_occupations)
448
449 if (this%method == smear_fixed_occ) then
450 ! do nothing
451 else if (this%method == smear_semiconductor) then
452 assert(this%fermi_count > 0 .and. this%fermi_count <= nik*nst)
453
454 ifermi = 0
455 do ik = 1, nik
456 do ist = 1, 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
462 ifermi = ifermi + 1
463 else
464 occupations(ist, ik) = m_zero
465 end if
466
467 end do
468 end do
469
470 else
471 dsmear = max(1e-14_real64, this%dsmear)
472 if (this%photodop) then
473 dsmear_cond = max(1e-14_real64, this%dsmear_cond)
474 do ik = 1, nik
475 do ist = 1, nst
476 if (ist < this%photodop_bandmin) then
477 ! valence electrons
478 xx = (this%e_fermi - eigenvalues(ist, ik))/dsmear
479 else
480 ! conduction electrons
481 xx = (this%e_fermi_cond - eigenvalues(ist, ik))/dsmear_cond
482 end if
483 occupations(ist, ik) = smear_step_function(this, xx) * this%el_per_state
484 end do
485 end do
486 else
487 do ik = 1, nik
488 do ist = 1, nst
489 xx = (this%e_fermi - eigenvalues(ist, ik))/dsmear
490 occupations(ist, ik) = smear_step_function(this, xx) * this%el_per_state
491 end do
492 end do
493 end if
494 end if
495
497 end subroutine smear_fill_occupations
498
499
500 !--------------------------------------------------
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(:, :)
508
509 integer :: ist, ik
510 real(real64) :: dsmear, xx, term, ff
511
512 push_sub(smear_calc_entropy)
513
514 entropy = m_zero
515
516 if (this%method == smear_fixed_occ .or. this%method == smear_semiconductor) then
517 ! Fermi-Dirac entropy, not quite the same as will be obtained with true smearing
518 ! RM Wentzcovitch, JL Martins, and PB Allen, Phys. Rev. B 45, 11372 (1992) eqn (5)
519 ! also N Marzari PhD thesis p 117, http://quasiamore.mit.edu/phd/Marzari_PhD.pdf
520 do ik = 1, nik
521 do ist = 1, nst
522 ff = occ(ist, ik) / this%el_per_state
523 if (ff > m_zero .and. ff < m_one) then
524 term = ff * log(ff) + (1 - ff) * log(1 - ff)
525 else ! we have semiconducting smearing, or perverse occupations as in Methfessel-Paxton
526 term = m_zero
527 end if
528 entropy = entropy - kweights(ik) * this%el_per_state * term
529 end do
530 end do
531 else
532 dsmear = max(1e-14_real64, this%dsmear)
533
534 do ik = 1, nik
535 do ist = 1, nst
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 * &
539 smear_entropy_function(this, xx)
540 end if
541 end do
542 end do
543 end if
544
545 pop_sub(smear_calc_entropy)
546 end function smear_calc_entropy
547
548
549 ! ---------------------------------------------------------
550 real(real64) function smear_delta_function(this, xx) result(deltaf)
551 type(smear_t), intent(in) :: this
552 real(real64), intent(in) :: xx
553
554 real(real64), parameter :: maxarg = 200.0_real64
555 real(real64) :: xp, arg, hd, hp, aa
556 integer :: ii, ni
557
558 ! no PUSH_SUB, called too often
559
560 ! smear_delta_function is not defined for SMEAR_FIXED_OCC
561 assert(this%method /= smear_fixed_occ)
562
563 deltaf = m_zero
564 select case (this%method)
566 if (abs(xx) <= m_epsilon) then
567 deltaf = this%ef_occ
568 end if
569
570 case (smear_fermi_dirac)
571 if (abs(xx) <= 36.0_real64) then
572 deltaf = m_one / (m_two + exp(-xx) + exp(xx))
573 end if
574
575 case (smear_cold)
576 xp = xx - m_one / sqrt(m_two)
577 arg = min(maxarg, xp**2)
578
579 deltaf = exp(-arg) / sqrt(m_pi) * (m_two - sqrt(m_two) * xx)
580
582 arg = min(maxarg, xx**2)
583 deltaf = exp(-arg) / sqrt(m_pi)
584
585 if (this%MP_n > 0) then ! recursion
586 hd = m_zero
587 hp = exp(-arg)
588 ni = 0
589 aa = m_one / sqrt(m_pi)
590 do ii = 1, this%MP_n
591 hd = m_two * xx * hp - m_two * ni * hd
592 ni = ni + 1
593 aa = -aa / (m_four * ii)
594 hp = m_two * xx * hd - m_two * ni * hp
595 ni = ni + 1
596 deltaf = deltaf + aa * hp
597 end do
598 end if
599
600 case (smear_spline)
601 xp = abs(xx) + m_one / sqrt(m_two)
602 deltaf = sqrt(m_e) * xp * exp(-xp * xp)
603
604 end select
605
606 end function smear_delta_function
607
608
609 ! ---------------------------------------------------------
610 real(real64) function smear_step_function(this, xx) result(stepf)
611 type(smear_t), intent(in) :: this
612 real(real64), intent(in) :: xx
613
614 real(real64), parameter :: maxarg = 200.0_real64
615 real(real64) :: xp, arg, hd, hp, aa
616 integer :: ii, ni
617
618 push_sub(smear_step_function)
619
620 ! smear_step_function is not defined for SMEAR_FIXED_OCC
621 assert(this%method /= smear_fixed_occ)
622
623 stepf = m_zero
624 select case (this%method)
626 if (xx > m_zero) then
627 stepf = m_one
628 else if (abs(xx) <= m_epsilon) then
629 stepf = this%ef_occ
630 end if
631
632 case (smear_fermi_dirac)
633 if (xx > maxarg) then
634 stepf = m_one
635 else if (xx > -maxarg) then
636 stepf = m_one / (m_one + exp(-xx))
637 end if
638
639 case (smear_cold)
640 xp = xx - m_one / sqrt(m_two)
641 arg = min(maxarg, xp**2)
642
643 stepf = m_half * loct_erf(xp) + &
644 m_one / sqrt(m_two * m_pi) * exp(-arg) + m_half
645
647 stepf = m_half * loct_erfc(-xx)
648
649 if (this%MP_n > 0) then ! recursion
650 hd = m_zero
651 arg = min(maxarg, xx**2)
652 hp = exp(-arg)
653 ni = 0
654 aa = m_one / sqrt(m_pi)
655 do ii = 1, this%MP_n
656 hd = m_two * xx * hp - m_two * ni * hd
657 ni = ni + 1
658 aa = -aa / (m_four * ii)
659 stepf = stepf - aa * hd
660 hp = m_two * xx * hd - m_two * ni * hp
661 ni = ni + 1
662 end do
663 end if
664
665 case (smear_spline)
666 if (xx <= m_zero) then
667 xp = xx - m_one / sqrt(m_two)
668 stepf = m_half * sqrt(m_e) * exp(-xp * xp)
669 else
670 xp = xx + m_one / sqrt(m_two)
671 stepf = m_one - m_half * sqrt(m_e) * exp(-xp * xp)
672 end if
673
674 end select
675
676 pop_sub(smear_step_function)
677 end function smear_step_function
678
679
680 ! ---------------------------------------------------------
682 real(real64) function smear_entropy_function(this, xx) result(entropyf)
683 type(smear_t), intent(in) :: this
684 real(real64), intent(in) :: xx
685
686 real(real64), parameter :: maxarg = 200.0_real64
687 real(real64) :: xp, arg, hd, hp, hpm1, aa
688 integer :: ii, ni
689
690 push_sub(smear_entropy_function)
691
692 ! smear_entropy_function is not defined for SMEAR_FIXED_OCC
693 assert(this%method /= smear_fixed_occ)
694
695 entropyf = m_zero
696 select case (this%method)
698
699 case (smear_fermi_dirac)
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)
703 end if
704
705 case (smear_cold)
706 xp = xx - m_one / sqrt(m_two)
707 arg = min(maxarg, xp**2)
708
709 entropyf = m_one / sqrt(m_two * m_pi) * xp * exp(-arg)
710
712 arg = min(maxarg, xx**2)
713 entropyf = -m_half * exp(-arg) / sqrt(m_pi)
714
715 if (this%MP_n > 0) then ! recursion
716 hd = m_zero
717 hp = exp(-arg)
718 ni = 0
719 aa = m_one / sqrt(m_pi)
720 do ii = 1, this%MP_n
721 hd = m_two * xx * hp - m_two * ni * hd
722 ni = ni + 1
723 hpm1 = hp
724 hp = m_two * xx * hd - m_two * ni * hp
725 ni = ni + 1
726 aa = -aa / (m_four * ii)
727 entropyf = entropyf - aa * (m_half * hp + hpm1 * ni)
728 end do
729 end if
730
731 case (smear_spline)
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))
734
735 end select
736
738 end function smear_entropy_function
739
740
741 ! ---------------------------------------------------------
742 logical pure function smear_is_semiconducting(this) result(answer)
743 type(smear_t), intent(in) :: this
744
745 answer = this%method == smear_semiconductor
746
747 end function smear_is_semiconducting
748
749 subroutine smear_write_info(this, namespace, iunit)
750 type(smear_t), intent(in) :: this
751 type(namespace_t), intent(in) :: namespace
752 integer, optional, intent(in) :: iunit
753
754
755 if (this%method /= smear_semiconductor .and. this%method /= smear_fixed_occ) then
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)
762 else
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)
766 end if
767 end if
768
769 end subroutine smear_write_info
770
771
772end module smear_oct_m
773
774!! Local Variables:
775!! mode: f90
776!! coding: utf-8
777!! End:
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:149
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
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public p_ry
Definition: global.F90:221
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_min_occ
Minimal occupation that is considered to be non-zero.
Definition: global.F90:211
integer function, public kpoints_kweight_denominator(this)
Definition: kpoints.F90:1609
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
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
real(real64) function, public smear_calc_entropy(this, eigenvalues, nik, nst, kweights, occ)
Definition: smear.F90:596
subroutine, public smear_fill_occupations(this, eigenvalues, occupations, nik, nst)
Definition: smear.F90:532
subroutine, public smear_find_fermi_energy(this, namespace, eigenvalues, occupations, qtot, nik, nst, kweights)
Definition: smear.F90:344
real(real64) function, public smear_entropy_function(this, xx)
This function is defined as .
Definition: smear.F90:776
real(real64) function, public smear_delta_function(this, xx)
Definition: smear.F90:644
integer, parameter, public smear_fermi_dirac
Definition: smear.F90:171
integer, parameter, public smear_methfessel_paxton
Definition: smear.F90:171
subroutine, public smear_write_info(this, namespace, iunit)
Definition: smear.F90:843
real(real64) function, public smear_step_function(this, xx)
Definition: smear.F90:704
integer, parameter, public smear_semiconductor
Definition: smear.F90:171
subroutine, public smear_copy(to, from)
Definition: smear.F90:323
integer, parameter, public smear_fixed_occ
Definition: smear.F90:171
integer, parameter, public smear_spline
Definition: smear.F90:171
subroutine bisection_find_fermi_energy(this, namespace, dsmear_in, tol, eigenvalues, kweights, nik, q_in, start_band, end_band, e_fermi, ef_occ)
Definition: smear.F90:478
subroutine, public smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
Definition: smear.F90:185
integer, parameter, public smear_cold
Definition: smear.F90:171
logical pure function, public smear_is_semiconducting(this)
Definition: smear.F90:836
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
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
int true(void)