(git:1155b05)
Loading...
Searching...
No Matches
smearing_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Unified smearing module supporting four methods:
10!> smear_fermi_dirac — Fermi-Dirac distribution
11!> smear_gaussian — Gaussian broadening
12!> smear_mp — Methfessel-Paxton first order
13!> smear_mv — Marzari-Vanderbilt (cold smearing)
14!>
15!> All methods share the bisection framework, LFOMO/HOMO logic, and the
16!> analytical rank-1 Jacobian. Only the per-state math (f, kTS, g_i)
17!> differs, selected via a method integer from input_constants.
18!>
19!> \par History
20!> 09.2008: Created (fermi_utils.F)
21!> 02.2026: Extended to more smearing method and renamed
22!> \author Joost VandeVondele
23! **************************************************************************************************
24MODULE smearing_utils
25
26 USE input_constants, ONLY: smear_fermi_dirac,&
27 smear_gaussian,&
28 smear_mp,&
29 smear_mv
30 USE kahan_sum, ONLY: accurate_sum
31 USE kinds, ONLY: dp
32 USE mathconstants, ONLY: rootpi,&
33 sqrt2,&
34 sqrthalf
35#include "base/base_uses.f90"
36
37 IMPLICIT NONE
38
39 PRIVATE
40
41 ! Unified interface (method as parameter)
42 PUBLIC :: smearocc, smearfixed, smearfixedderiv, smearfixedderivmv
43 PUBLIC :: smearkp, smearkp2
44
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'smearing_utils'
46 INTEGER, PARAMETER, PRIVATE :: BISECT_MAX_ITER = 400
47
48CONTAINS
49
50! **************************************************************************************************
51!> \brief Returns occupations and smearing correction for a given set of
52!> energies and chemical potential, using one of four smearing methods.
53!>
54!> Fermi-Dirac: f_i = occ / [1 + exp((e_i - mu)/sigma)]
55!> Gaussian: f_i = (occ/2) * erfc[(e_i - mu)/sigma]
56!> MP-1: f_i = (occ/2) * erfc(x) - occ*x/(2*sqrt(pi)) * exp(-x^2)
57!> MV: f_i = (occ/2) * erfc(u) + occ/(sqrt(2*pi)) * exp(-u^2), u = x + 1/sqrt(2)
58!>
59!> kTS is the smearing correction to the free energy (physically -TS
60!> for Fermi-Dirac; a variational correction term for the other methods).
61!> It enters the total energy and the Gillan extrapolation E(0) = E - kTS/2.
62!>
63!> \param f occupations (output)
64!> \param N total number of electrons (output)
65!> \param kTS smearing correction to the free energy (output)
66!> \param e eigenvalues (input)
67!> \param mu chemical potential (input)
68!> \param sigma smearing width: kT for Fermi-Dirac, sigma for others (input)
69!> \param maxocc maximum occupation of an orbital (input)
70!> \param method smearing method selector from input_constants (input)
71!> \param estate excited state index for core-level spectroscopy (optional)
72!> \param festate occupation of the excited state (optional)
73! **************************************************************************************************
74 SUBROUTINE smearocc(f, N, kTS, e, mu, sigma, maxocc, method, estate, festate)
75
76 REAL(kind=dp), INTENT(OUT) :: f(:), n, kts
77 REAL(kind=dp), INTENT(IN) :: e(:), mu, sigma, maxocc
78 INTEGER, INTENT(IN) :: method
79 INTEGER, INTENT(IN), OPTIONAL :: estate
80 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
81
82 INTEGER :: i, nstate
83 REAL(kind=dp) :: arg, expu2, expx2, occupation, term1, &
84 term2, tmp, tmp2, tmp3, tmp4, tmplog, &
85 u, x
86
87 nstate = SIZE(e)
88 kts = 0.0_dp
89
90 DO i = 1, nstate
91 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
92 IF (i == estate) THEN
93 occupation = festate
94 ELSE
95 occupation = maxocc
96 END IF
97 ELSE
98 occupation = maxocc
99 END IF
100
101 SELECT CASE (method)
102 CASE (smear_fermi_dirac)
103 IF (e(i) > mu) THEN
104 arg = -(e(i) - mu)/sigma
105 tmp = exp(arg)
106 tmp4 = tmp + 1.0_dp
107 tmp2 = tmp/tmp4
108 tmp3 = 1.0_dp/tmp4
109 tmplog = -log(tmp4)
110 term1 = tmp2*(arg + tmplog)
111 term2 = tmp3*tmplog
112 ELSE
113 arg = (e(i) - mu)/sigma
114 tmp = exp(arg)
115 tmp4 = tmp + 1.0_dp
116 tmp2 = 1.0_dp/tmp4
117 tmp3 = tmp/tmp4
118 tmplog = -log(tmp4)
119 term1 = tmp2*tmplog
120 term2 = tmp3*(arg + tmplog)
121 END IF
122 f(i) = occupation*tmp2
123 kts = kts + sigma*occupation*(term1 + term2)
124
125 CASE (smear_gaussian)
126 x = (e(i) - mu)/sigma
127 expx2 = exp(-x*x)
128 f(i) = occupation*0.5_dp*erfc(x)
129 kts = kts - (sigma/(2.0_dp*rootpi))*occupation*expx2
130
131 CASE (smear_mp)
132 x = (e(i) - mu)/sigma
133 expx2 = exp(-x*x)
134 f(i) = occupation*(0.5_dp*erfc(x) - x/(2.0_dp*rootpi)*expx2)
135 kts = kts + (sigma/(4.0_dp*rootpi))*occupation*(2.0_dp*x*x - 1.0_dp)*expx2
136
137 CASE (smear_mv)
138 x = (e(i) - mu)/sigma
139 u = x + sqrthalf
140 expu2 = exp(-u*u)
141 f(i) = occupation*(0.5_dp*erfc(u) + expu2/(sqrt2*rootpi))
142 kts = kts - (sigma/(sqrt2*rootpi))*occupation*u*expu2
143
144 CASE DEFAULT
145 cpabort("SmearOcc: unknown smearing method")
146 END SELECT
147 END DO
148
149 n = accurate_sum(f)
150
151 END SUBROUTINE smearocc
152
153! **************************************************************************************************
154!> \brief k-point version of SmearOcc (module-private).
155!> Computes occupations and kTS for a 2D array of eigenvalues
156!> (nmo x nkp) weighted by k-point weights.
157!> Falls back to a step function when sigma < 1e-14.
158!>
159!> \param f occupations (nmo x nkp, output)
160!> \param nel total number of electrons (output)
161!> \param kTS smearing correction (output)
162!> \param e eigenvalues (nmo x nkp, input)
163!> \param mu chemical potential (input)
164!> \param wk k-point weights (input)
165!> \param sigma smearing width (input)
166!> \param maxocc maximum occupation (input)
167!> \param method smearing method selector (input)
168! **************************************************************************************************
169 SUBROUTINE smear2(f, nel, kTS, e, mu, wk, sigma, maxocc, method)
170
171 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: f
172 REAL(kind=dp), INTENT(OUT) :: nel, kts
173 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: e
174 REAL(kind=dp), INTENT(IN) :: mu
175 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wk
176 REAL(kind=dp), INTENT(IN) :: sigma, maxocc
177 INTEGER, INTENT(IN) :: method
178
179 INTEGER :: ik, is, nkp, nmo
180 REAL(kind=dp) :: arg, expu2, expx2, term1, term2, tmp, &
181 tmp2, tmp3, tmp4, tmplog, u, x
182
183 nmo = SIZE(e, 1)
184 nkp = SIZE(e, 2)
185 kts = 0.0_dp
186
187 IF (sigma > 1.0e-14_dp) THEN
188 DO ik = 1, nkp
189 DO is = 1, nmo
190 SELECT CASE (method)
191 CASE (smear_fermi_dirac)
192 IF (e(is, ik) > mu) THEN
193 arg = -(e(is, ik) - mu)/sigma
194 tmp = exp(arg)
195 tmp4 = tmp + 1.0_dp
196 tmp2 = tmp/tmp4
197 tmp3 = 1.0_dp/tmp4
198 tmplog = -log(tmp4)
199 term1 = tmp2*(arg + tmplog)
200 term2 = tmp3*tmplog
201 ELSE
202 arg = (e(is, ik) - mu)/sigma
203 tmp = exp(arg)
204 tmp4 = tmp + 1.0_dp
205 tmp2 = 1.0_dp/tmp4
206 tmp3 = tmp/tmp4
207 tmplog = -log(tmp4)
208 term1 = tmp2*tmplog
209 term2 = tmp3*(arg + tmplog)
210 END IF
211 f(is, ik) = maxocc*tmp2
212 kts = kts + sigma*maxocc*(term1 + term2)*wk(ik)
213
214 CASE (smear_gaussian)
215 x = (e(is, ik) - mu)/sigma
216 expx2 = exp(-x*x)
217 f(is, ik) = maxocc*0.5_dp*erfc(x)
218 kts = kts - (sigma/(2.0_dp*rootpi))*maxocc*expx2*wk(ik)
219
220 CASE (smear_mp)
221 x = (e(is, ik) - mu)/sigma
222 expx2 = exp(-x*x)
223 f(is, ik) = maxocc*(0.5_dp*erfc(x) - x/(2.0_dp*rootpi)*expx2)
224 kts = kts + (sigma/(4.0_dp*rootpi))*maxocc*(2.0_dp*x*x - 1.0_dp)*expx2*wk(ik)
225
226 CASE (smear_mv)
227 x = (e(is, ik) - mu)/sigma
228 u = x + sqrthalf
229 expu2 = exp(-u*u)
230 f(is, ik) = maxocc*(0.5_dp*erfc(u) + expu2/(sqrt2*rootpi))
231 kts = kts - (sigma/(sqrt2*rootpi))*maxocc*u*expu2*wk(ik)
232
233 CASE DEFAULT
234 cpabort("Smear2: unknown smearing method")
235 END SELECT
236 END DO
237 END DO
238 ELSE
239 ! Zero-width limit: step function
240 DO ik = 1, nkp
241 DO is = 1, nmo
242 IF (e(is, ik) <= mu) THEN
243 f(is, ik) = maxocc
244 ELSE
245 f(is, ik) = 0.0_dp
246 END IF
247 END DO
248 END DO
249 END IF
250
251 nel = 0.0_dp
252 DO ik = 1, nkp
253 nel = nel + accurate_sum(f(1:nmo, ik))*wk(ik)
254 END DO
255
256 END SUBROUTINE smear2
257
258! **************************************************************************************************
259!> \brief Bisection search for the chemical potential mu such that the total
260!> electron count equals N, for a given smearing method (Gamma point).
261!> Brackets mu by expanding outward from [min(e), max(e)] in steps
262!> of sigma, then bisects to machine precision.
263!> Could fail if mu lies far outside the eigenvalue range (to be fixed).
264!>
265!> \param f occupations (output)
266!> \param mu chemical potential found by bisection (output)
267!> \param kTS smearing correction (output)
268!> \param e eigenvalues (input)
269!> \param N target number of electrons (input)
270!> \param sigma smearing width (input)
271!> \param maxocc maximum occupation (input)
272!> \param method smearing method selector (input)
273!> \param estate excited state index for core-level spectroscopy (optional)
274!> \param festate occupation of the excited state (optional)
275! **************************************************************************************************
276 SUBROUTINE smearfixed(f, mu, kTS, e, N, sigma, maxocc, method, estate, festate)
277
278 REAL(kind=dp), INTENT(OUT) :: f(:), mu, kts
279 REAL(kind=dp), INTENT(IN) :: e(:), n, sigma, maxocc
280 INTEGER, INTENT(IN) :: method
281 INTEGER, INTENT(IN), OPTIONAL :: estate
282 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
283
284 INTEGER :: iter, my_estate
285 REAL(kind=dp) :: mu_max, mu_min, mu_now, my_festate, &
286 n_now, n_tmp
287
288 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
289 my_estate = estate
290 my_festate = festate
291 ELSE
292 my_estate = nint(maxocc)
293 my_festate = my_estate
294 END IF
295
296 ! Bracket mu from below
297 mu_min = minval(e)
298 iter = 0
299 DO
300 iter = iter + 1
301 CALL smearocc(f, n_tmp, kts, e, mu_min, sigma, maxocc, method, my_estate, my_festate)
302 IF (n_tmp > n .OR. iter > 20) THEN
303 mu_min = mu_min - sigma
304 ELSE
305 EXIT
306 END IF
307 END DO
308
309 ! Bracket mu from above
310 mu_max = maxval(e)
311 iter = 0
312 DO
313 iter = iter + 1
314 CALL smearocc(f, n_tmp, kts, e, mu_max, sigma, maxocc, method, my_estate, my_festate)
315 IF (n_tmp < n .OR. iter > 20) THEN
316 mu_max = mu_max + sigma
317 ELSE
318 EXIT
319 END IF
320 END DO
321
322 ! Bisection
323 iter = 0
324 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
325 iter = iter + 1
326 mu_now = (mu_max + mu_min)/2.0_dp
327 CALL smearocc(f, n_now, kts, e, mu_now, sigma, maxocc, method, my_estate, my_festate)
328
329 IF (n_now <= n) THEN
330 mu_min = mu_now
331 ELSE
332 mu_max = mu_now
333 END IF
334
335 IF (iter > bisect_max_iter) THEN
336 cpwarn("SmearFixed: maximum bisection iterations reached")
337 EXIT
338 END IF
339 END DO
340
341 mu = (mu_max + mu_min)/2.0_dp
342 CALL smearocc(f, n_now, kts, e, mu, sigma, maxocc, method, my_estate, my_festate)
343
344 END SUBROUTINE smearfixed
345
346! **************************************************************************************************
347!> \brief Bisection search for mu given a target electron count (k-point case,
348!> single spin channel or spin-degenerate).
349!> Initial bracket width is max(10*sigma, 0.5) for Gaussian/MP/MV,
350!> or sigma*ln[(1-eps)/eps] for Fermi-Dirac, reflecting the different
351!> tail decay rates.
352!>
353!> \param f occupations (nmo x nkp, output)
354!> \param mu chemical potential (output)
355!> \param kTS smearing correction (output)
356!> \param e eigenvalues (nmo x nkp, input)
357!> \param nel target number of electrons (input)
358!> \param wk k-point weights (input)
359!> \param sigma smearing width (input)
360!> \param maxocc maximum occupation (input)
361!> \param method smearing method selector (input)
362! **************************************************************************************************
363 SUBROUTINE smearkp(f, mu, kTS, e, nel, wk, sigma, maxocc, method)
364
365 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: f
366 REAL(kind=dp), INTENT(OUT) :: mu, kts
367 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: e
368 REAL(kind=dp), INTENT(IN) :: nel
369 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wk
370 REAL(kind=dp), INTENT(IN) :: sigma, maxocc
371 INTEGER, INTENT(IN) :: method
372
373 REAL(kind=dp), PARAMETER :: epsocc = 1.0e-12_dp
374
375 INTEGER :: iter
376 REAL(kind=dp) :: de, mu_max, mu_min, n_now
377
378 SELECT CASE (method)
379 CASE (smear_fermi_dirac)
380 de = sigma*log((1.0_dp - epsocc)/epsocc)
381 CASE DEFAULT
382 de = 10.0_dp*sigma
383 END SELECT
384 de = max(de, 0.5_dp)
385
386 mu_min = minval(e) - de
387 mu_max = maxval(e) + de
388 iter = 0
389 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
390 iter = iter + 1
391 mu = (mu_max + mu_min)/2.0_dp
392 CALL smear2(f, n_now, kts, e, mu, wk, sigma, maxocc, method)
393
394 IF (abs(n_now - nel) < nel*epsocc) EXIT
395
396 IF (n_now <= nel) THEN
397 mu_min = mu
398 ELSE
399 mu_max = mu
400 END IF
401
402 IF (iter > bisect_max_iter) THEN
403 cpwarn("Smearkp: maximum bisection iterations reached")
404 EXIT
405 END IF
406 END DO
407
408 mu = (mu_max + mu_min)/2.0_dp
409 CALL smear2(f, n_now, kts, e, mu, wk, sigma, maxocc, method)
410
411 END SUBROUTINE smearkp
412
413! **************************************************************************************************
414!> \brief Bisection search for mu (k-point, spin-polarised with a shared
415!> chemical potential across both spin channels).
416!> Asserts that the third dimension of f and e is exactly 2.
417!>
418!> \param f occupations (nmo x nkp x 2, output)
419!> \param mu chemical potential (output)
420!> \param kTS smearing correction (output)
421!> \param e eigenvalues (nmo x nkp x 2, input)
422!> \param nel target total number of electrons (input)
423!> \param wk k-point weights (input)
424!> \param sigma smearing width (input)
425!> \param method smearing method selector (input)
426! **************************************************************************************************
427 SUBROUTINE smearkp2(f, mu, kTS, e, nel, wk, sigma, method)
428
429 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: f
430 REAL(kind=dp), INTENT(OUT) :: mu, kts
431 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: e
432 REAL(kind=dp), INTENT(IN) :: nel
433 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wk
434 REAL(kind=dp), INTENT(IN) :: sigma
435 INTEGER, INTENT(IN) :: method
436
437 REAL(kind=dp), PARAMETER :: epsocc = 1.0e-12_dp
438
439 INTEGER :: iter
440 REAL(kind=dp) :: de, ktsa, ktsb, mu_max, mu_min, n_now, &
441 na, nb
442
443 cpassert(SIZE(f, 3) == 2 .AND. SIZE(e, 3) == 2)
444
445 SELECT CASE (method)
446 CASE (smear_fermi_dirac)
447 de = sigma*log((1.0_dp - epsocc)/epsocc)
448 CASE DEFAULT
449 de = 10.0_dp*sigma
450 END SELECT
451 de = max(de, 0.5_dp)
452
453 mu_min = minval(e) - de
454 mu_max = maxval(e) + de
455 iter = 0
456 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
457 iter = iter + 1
458 mu = (mu_max + mu_min)/2.0_dp
459 CALL smear2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
460 CALL smear2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
461 n_now = na + nb
462
463 IF (abs(n_now - nel) < nel*epsocc) EXIT
464
465 IF (n_now <= nel) THEN
466 mu_min = mu
467 ELSE
468 mu_max = mu
469 END IF
470
471 IF (iter > bisect_max_iter) THEN
472 cpwarn("Smearkp2: maximum bisection iterations reached")
473 EXIT
474 END IF
475 END DO
476
477 mu = (mu_max + mu_min)/2.0_dp
478 CALL smear2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
479 CALL smear2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
480 kts = ktsa + ktsb
481
482 END SUBROUTINE smearkp2
483
484! **************************************************************************************************
485!> \brief Computes the smearing weight vector g_i = -df_i/de_i with mu held
486!> fixed (module-private helper for the analytical Jacobian routines).
487!>
488!> Fermi-Dirac: g_i = occ * f_norm * (1 - f_norm) / sigma
489!> where f_norm = f_i/occ_i (overflow-safe, uses
490!> pre-computed f rather than re-evaluating exp)
491!> Gaussian: g_i = occ / (sigma*sqrt(pi)) * exp(-x^2)
492!> MP-1: g_i = occ * (3 - 2*x^2) / (2*sigma*sqrt(pi)) * exp(-x^2)
493!> MV: g_i = occ * (2 + sqrt(2)*x) / (sigma*sqrt(pi)) * exp(-u^2)
494!>
495!> Note: g_i can be negative for MP-1 (|x| > sqrt(3/2)) and MV
496!> (x < -sqrt(2)). This does not break bisection (N(mu) is still
497!> monotone) but the Jacobian routines must guard against G = sum(g_i)
498!> being near zero.
499!>
500!> \param gvec weight vector (Nstate, output)
501!> \param f occupations from a prior SmearOcc/SmearFixed call (input)
502!> \param e eigenvalues (input)
503!> \param mu chemical potential (input)
504!> \param sigma smearing width (input)
505!> \param maxocc maximum occupation (input)
506!> \param Nstate number of states (input)
507!> \param method smearing method selector (input)
508!> \param estate excited state index (optional)
509!> \param festate occupation of the excited state (optional)
510! **************************************************************************************************
511 SUBROUTINE compute_gvec(gvec, f, e, mu, sigma, maxocc, Nstate, method, estate, festate)
512
513 REAL(kind=dp), INTENT(OUT) :: gvec(:)
514 REAL(kind=dp), INTENT(IN) :: f(:), e(:), mu, sigma, maxocc
515 INTEGER, INTENT(IN) :: nstate, method
516 INTEGER, INTENT(IN), OPTIONAL :: estate
517 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
518
519 INTEGER :: i
520 REAL(kind=dp) :: expu2, expx2, fi_norm, occ_i, u, x
521
522 DO i = 1, nstate
523 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
524 IF (i == estate) THEN
525 occ_i = festate
526 ELSE
527 occ_i = maxocc
528 END IF
529 ELSE
530 occ_i = maxocc
531 END IF
532
533 IF (occ_i < epsilon(occ_i)) THEN
534 gvec(i) = 0.0_dp
535 cycle
536 END IF
537
538 x = (e(i) - mu)/sigma
539
540 SELECT CASE (method)
541 CASE (smear_fermi_dirac)
542 fi_norm = f(i)/occ_i
543 gvec(i) = occ_i*fi_norm*(1.0_dp - fi_norm)/sigma
544
545 CASE (smear_gaussian)
546 expx2 = exp(-x*x)
547 gvec(i) = occ_i/(sigma*rootpi)*expx2
548
549 CASE (smear_mp)
550 expx2 = exp(-x*x)
551 gvec(i) = occ_i*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2
552
553 CASE (smear_mv)
554 u = x + sqrthalf
555 expu2 = exp(-u*u)
556 gvec(i) = occ_i*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2
557
558 END SELECT
559 END DO
560
561 END SUBROUTINE compute_gvec
562
563! **************************************************************************************************
564!> \brief Analytical Jacobian df_i/de_j for any smearing method under the
565!> electron-number constraint sum(f) = N.
566!>
567!> Differentiating f_i(e, mu(e)) where mu is implicitly defined by
568!> the constraint yields:
569!>
570!> df_i/de_j = -delta_{ij} * g_i + g_i * g_j / G
571!>
572!> where g_i = -df_i/de_i (mu fixed) and G = sum(g_i).
573!> This is a diagonal matrix plus a symmetric rank-1 update.
574!> Building it costs O(N) for g, plus O(N^2) for the outer product.
575!>
576!> Replaces the original numerical finite-difference FermiFixedDeriv
577!> which required 2N bisection solves. Exact to machine precision
578!> for all four methods.
579!>
580!> \param dfde Jacobian matrix dfde(i,j) = df_i/de_j (Nstate x Nstate, output)
581!> \param f occupations (output)
582!> \param mu chemical potential (output)
583!> \param kTS smearing correction (output)
584!> \param e eigenvalues (input)
585!> \param N target number of electrons (input)
586!> \param sigma smearing width (input)
587!> \param maxocc maximum occupation (input)
588!> \param method smearing method selector (input)
589!> \param estate excited state index (optional)
590!> \param festate occupation of the excited state (optional)
591! **************************************************************************************************
592 SUBROUTINE smearfixedderiv(dfde, f, mu, kTS, e, N, sigma, maxocc, method, estate, festate)
593
594 REAL(kind=dp), INTENT(OUT) :: dfde(:, :), f(:), mu, kts
595 REAL(kind=dp), INTENT(IN) :: e(:), n, sigma, maxocc
596 INTEGER, INTENT(IN) :: method
597 INTEGER, INTENT(IN), OPTIONAL :: estate
598 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
599
600 CHARACTER(len=*), PARAMETER :: routinen = 'SmearFixedDeriv'
601
602 INTEGER :: handle, i, j, nstate
603 REAL(kind=dp) :: gsum
604 REAL(kind=dp), ALLOCATABLE :: gvec(:)
605
606 CALL timeset(routinen, handle)
607
608 ! Step 1: find mu and f
609 CALL smearfixed(f, mu, kts, e, n, sigma, maxocc, method, estate, festate)
610
611 ! Step 2: build g vector
612 nstate = SIZE(e)
613 ALLOCATE (gvec(nstate))
614 CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, nstate, method, estate, festate)
615 gsum = accurate_sum(gvec)
616
617 ! Step 3: assemble dfde(i,j) = -delta_{ij}*g_i + g_i*g_j/G
618 IF (abs(gsum) > epsilon(gsum)) THEN
619 DO j = 1, nstate
620 DO i = 1, nstate
621 dfde(i, j) = gvec(i)*gvec(j)/gsum
622 END DO
623 dfde(j, j) = dfde(j, j) - gvec(j)
624 END DO
625 ELSE
626 dfde(:, :) = 0.0_dp
627 END IF
628
629 DEALLOCATE (gvec)
630 CALL timestop(handle)
631
632 END SUBROUTINE smearfixedderiv
633
634! **************************************************************************************************
635!> \brief Apply TRANSPOSE(df/de) to a vector WITHOUT forming the full N x N
636!> Jacobian. O(N) time and O(N) memory for all four methods.
637!>
638!> Exploiting the rank-1 structure of the constrained Jacobian:
639!>
640!> [J^T v]_j = g_j * (g . v / G - v_j)
641!>
642!> This replaces the pattern used in qs_mo_occupation:
643!> ALLOCATE(dfde(nmo,nmo))
644!> CALL SmearFixedDeriv(dfde, ...)
645!> RESULT = MATMUL(TRANSPOSE(dfde), v)
646!> DEALLOCATE(dfde)
647!> turning O(N^2) storage + O(N^2) MATMUL into O(N) throughout.
648!>
649!> Currently the sole caller (qs_ot_scf do_ener) is dead code, but
650!> this routine is ready for when it is enabled.
651!>
652!> \param RESULT output vector = TRANSPOSE(df/de) * v (Nstate, output)
653!> \param f occupations (output)
654!> \param mu chemical potential (output)
655!> \param kTS smearing correction (output)
656!> \param e eigenvalues (input)
657!> \param N_el target number of electrons (input)
658!> \param sigma smearing width (input)
659!> \param maxocc maximum occupation (input)
660!> \param method smearing method selector (input)
661!> \param v input vector to multiply (Nstate, input)
662!> \param estate excited state index (optional)
663!> \param festate occupation of the excited state (optional)
664! **************************************************************************************************
665 SUBROUTINE smearfixedderivmv(RESULT, f, mu, kTS, e, N_el, sigma, maxocc, method, v, estate, festate)
666
667 REAL(kind=dp), INTENT(OUT) :: result(:), f(:), mu, kts
668 REAL(kind=dp), INTENT(IN) :: e(:), n_el, sigma, maxocc
669 INTEGER, INTENT(IN) :: method
670 REAL(kind=dp), INTENT(IN) :: v(:)
671 INTEGER, INTENT(IN), OPTIONAL :: estate
672 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
673
674 CHARACTER(len=*), PARAMETER :: routinen = 'SmearFixedDerivMV'
675
676 INTEGER :: handle, i, nstate
677 REAL(kind=dp) :: gdotv, gsum
678 REAL(kind=dp), ALLOCATABLE :: gvec(:)
679
680 CALL timeset(routinen, handle)
681
682 ! Step 1: find mu and f
683 CALL smearfixed(f, mu, kts, e, n_el, sigma, maxocc, method, estate, festate)
684
685 ! Step 2: build g vector
686 nstate = SIZE(e)
687 ALLOCATE (gvec(nstate))
688 CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, nstate, method, estate, festate)
689 gsum = accurate_sum(gvec)
690
691 ! Step 3: RESULT_j = g_j * (g.v / G - v_j)
692 IF (abs(gsum) > epsilon(gsum)) THEN
693 gdotv = 0.0_dp
694 DO i = 1, nstate
695 gdotv = gdotv + gvec(i)*v(i)
696 END DO
697 DO i = 1, nstate
698 result(i) = gvec(i)*(gdotv/gsum - v(i))
699 END DO
700 ELSE
701 result(:) = 0.0_dp
702 END IF
703
704 DEALLOCATE (gvec)
705 CALL timestop(handle)
706
707 END SUBROUTINE smearfixedderivmv
708
709END MODULE smearing_utils
kahan_sum::accurate_sum
Definition kahan_sum.F:41
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::smear_fermi_dirac
integer, parameter, public smear_fermi_dirac
Definition input_constants.F:905
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:905
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:905
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:905
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::sqrthalf
real(kind=dp), parameter, public sqrthalf
Definition mathconstants.F:127
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::sqrt2
real(kind=dp), parameter, public sqrt2
Definition mathconstants.F:127
smearing_utils
Unified smearing module supporting four methods: smear_fermi_dirac — Fermi-Dirac distribution smear_g...
Definition smearing_utils.F:24
smearing_utils::smearkp
subroutine, public smearkp(f, mu, kts, e, nel, wk, sigma, maxocc, method)
Bisection search for mu given a target electron count (k-point case, single spin channel or spin-dege...
Definition smearing_utils.F:364
smearing_utils::smearkp2
subroutine, public smearkp2(f, mu, kts, e, nel, wk, sigma, method)
Bisection search for mu (k-point, spin-polarised with a shared chemical potential across both spin ch...
Definition smearing_utils.F:428
smearing_utils::smearfixedderiv
subroutine, public smearfixedderiv(dfde, f, mu, kts, e, n, sigma, maxocc, method, estate, festate)
Analytical Jacobian df_i/de_j for any smearing method under the electron-number constraint sum(f) = N...
Definition smearing_utils.F:593
smearing_utils::smearocc
subroutine, public smearocc(f, n, kts, e, mu, sigma, maxocc, method, estate, festate)
Returns occupations and smearing correction for a given set of energies and chemical potential,...
Definition smearing_utils.F:75
smearing_utils::smearfixed
subroutine, public smearfixed(f, mu, kts, e, n, sigma, maxocc, method, estate, festate)
Bisection search for the chemical potential mu such that the total electron count equals N,...
Definition smearing_utils.F:277
smearing_utils::smearfixedderivmv
subroutine, public smearfixedderivmv(result, f, mu, kts, e, n_el, sigma, maxocc, method, v, estate, festate)
Apply TRANSPOSE(df/de) to a vector WITHOUT forming the full N x N Jacobian. O(N) time and O(N) memory...
Definition smearing_utils.F:666