(git:495eafe)
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 bibliography, ONLY: fuho1983,&
27 marzari1999,&
28 mermin1965,&
29 methfesselpaxton1989,&
30 cite_reference,&
31 dossantos2023
32 USE input_constants, ONLY: smear_fermi_dirac,&
33 smear_gaussian,&
34 smear_mp,&
35 smear_mv
36 USE kahan_sum, ONLY: accurate_sum
37 USE kinds, ONLY: dp
38 USE mathconstants, ONLY: rootpi,&
39 sqrt2,&
40 sqrthalf
41#include "base/base_uses.f90"
42
43 IMPLICIT NONE
44
45 PRIVATE
46
47 ! Unified interface (method as parameter)
48 PUBLIC :: smearocc, smearfixed, smearfixedderiv, smearfixedderivmv
49 PUBLIC :: smearkp, smearkp2
50 PRIVATE :: cite_smearing
51
52 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'smearing_utils'
53 INTEGER, PARAMETER, PRIVATE :: BISECT_MAX_ITER = 400
54 INTEGER, PARAMETER, PRIVATE :: NEWTON_MAX_ITER = 50
55
56CONTAINS
57! **************************************************************************************************
58!> \brief Citation of Smearing methods
59!> \param method ...
60! **************************************************************************************************
61 SUBROUTINE cite_smearing(method)
62 INTEGER, INTENT(IN) :: method
63
64 SELECT CASE (method)
65 CASE (smear_fermi_dirac)
66 CALL cite_reference(mermin1965)
67 CASE (smear_gaussian)
68 CALL cite_reference(fuho1983)
69 CASE (smear_mp)
70 CALL cite_reference(fuho1983)
71 CALL cite_reference(methfesselpaxton1989)
72 CALL cite_reference(dossantos2023)
73 CASE (smear_mv)
74 CALL cite_reference(fuho1983)
75 CALL cite_reference(marzari1999)
76 CALL cite_reference(dossantos2023)
77 END SELECT
78 END SUBROUTINE cite_smearing
79
80! **************************************************************************************************
81!> \brief Returns occupations and smearing correction for a given set of
82!> energies and chemical potential, using one of four smearing methods.
83!>
84!> Fermi-Dirac: f_i = occ / [1 + exp((e_i - mu)/sigma)]
85!> Gaussian: f_i = (occ/2) * erfc[(e_i - mu)/sigma]
86!> MP-1: f_i = (occ/2) * erfc(x) - occ*x/(2*sqrt(pi)) * exp(-x^2)
87!> MV: f_i = (occ/2) * erfc(u) + occ/(sqrt(2*pi)) * exp(-u^2), u = x + 1/sqrt(2)
88!>
89!> kTS is the smearing correction to the free energy (physically -TS
90!> for Fermi-Dirac; a variational correction term for the other methods).
91!> It enters the total energy and the Gillan extrapolation E(0) = E - kTS/2.
92!>
93!> \param f occupations (output)
94!> \param N total number of electrons (output)
95!> \param kTS smearing correction to the free energy (output)
96!> \param e eigenvalues (input)
97!> \param mu chemical potential (input)
98!> \param sigma smearing width: kT for Fermi-Dirac, sigma for others (input)
99!> \param maxocc maximum occupation of an orbital (input)
100!> \param method smearing method selector from input_constants (input)
101!> \param estate excited state index for core-level spectroscopy (optional)
102!> \param festate occupation of the excited state (optional)
103! **************************************************************************************************
104 SUBROUTINE smearocc(f, N, kTS, e, mu, sigma, maxocc, method, estate, festate)
105
106 REAL(kind=dp), INTENT(OUT) :: f(:), n, kts
107 REAL(kind=dp), INTENT(IN) :: e(:), mu, sigma, maxocc
108 INTEGER, INTENT(IN) :: method
109 INTEGER, INTENT(IN), OPTIONAL :: estate
110 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
111
112 INTEGER :: i, nstate
113 REAL(kind=dp) :: arg, expu2, expx2, occupation, term1, &
114 term2, tmp, tmp2, tmp3, tmp4, tmplog, &
115 u, x
116
117 nstate = SIZE(e)
118 kts = 0.0_dp
119
120 DO i = 1, nstate
121 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
122 IF (i == estate) THEN
123 occupation = festate
124 ELSE
125 occupation = maxocc
126 END IF
127 ELSE
128 occupation = maxocc
129 END IF
130
131 SELECT CASE (method)
132 CASE (smear_fermi_dirac)
133 IF (e(i) > mu) THEN
134 arg = -(e(i) - mu)/sigma
135 tmp = exp(arg)
136 tmp4 = tmp + 1.0_dp
137 tmp2 = tmp/tmp4
138 tmp3 = 1.0_dp/tmp4
139 tmplog = -log(tmp4)
140 term1 = tmp2*(arg + tmplog)
141 term2 = tmp3*tmplog
142 ELSE
143 arg = (e(i) - mu)/sigma
144 tmp = exp(arg)
145 tmp4 = tmp + 1.0_dp
146 tmp2 = 1.0_dp/tmp4
147 tmp3 = tmp/tmp4
148 tmplog = -log(tmp4)
149 term1 = tmp2*tmplog
150 term2 = tmp3*(arg + tmplog)
151 END IF
152 f(i) = occupation*tmp2
153 kts = kts + sigma*occupation*(term1 + term2)
154
155 CASE (smear_gaussian)
156 x = (e(i) - mu)/sigma
157 expx2 = exp(-x*x)
158 f(i) = occupation*0.5_dp*erfc(x)
159 kts = kts - (sigma/(2.0_dp*rootpi))*occupation*expx2
160
161 CASE (smear_mp)
162 x = (e(i) - mu)/sigma
163 expx2 = exp(-x*x)
164 f(i) = occupation*(0.5_dp*erfc(x) - x/(2.0_dp*rootpi)*expx2)
165 kts = kts + (sigma/(4.0_dp*rootpi))*occupation*(2.0_dp*x*x - 1.0_dp)*expx2
166
167 CASE (smear_mv)
168 x = (e(i) - mu)/sigma
169 u = x + sqrthalf
170 expu2 = exp(-u*u)
171 f(i) = occupation*(0.5_dp*erfc(u) + expu2/(sqrt2*rootpi))
172 kts = kts - (sigma/(sqrt2*rootpi))*occupation*u*expu2
173
174 CASE DEFAULT
175 cpabort("SmearOcc: unknown smearing method")
176 END SELECT
177 END DO
178
179 n = accurate_sum(f)
180
181 END SUBROUTINE smearocc
182
183! **************************************************************************************************
184!> \brief k-point version of SmearOcc (module-private).
185!> Computes occupations and kTS for a 2D array of eigenvalues
186!> (nmo x nkp) weighted by k-point weights.
187!> Falls back to a step function when sigma < 1e-14.
188!>
189!> \param f occupations (nmo x nkp, output)
190!> \param nel total number of electrons (output)
191!> \param kTS smearing correction (output)
192!> \param e eigenvalues (nmo x nkp, input)
193!> \param mu chemical potential (input)
194!> \param wk k-point weights (input)
195!> \param sigma smearing width (input)
196!> \param maxocc maximum occupation (input)
197!> \param method smearing method selector (input)
198! **************************************************************************************************
199 SUBROUTINE smear2(f, nel, kTS, e, mu, wk, sigma, maxocc, method)
200
201 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: f
202 REAL(kind=dp), INTENT(OUT) :: nel, kts
203 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: e
204 REAL(kind=dp), INTENT(IN) :: mu
205 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wk
206 REAL(kind=dp), INTENT(IN) :: sigma, maxocc
207 INTEGER, INTENT(IN) :: method
208
209 INTEGER :: ik, is, nkp, nmo
210 REAL(kind=dp) :: arg, expu2, expx2, term1, term2, tmp, &
211 tmp2, tmp3, tmp4, tmplog, u, x
212
213 nmo = SIZE(e, 1)
214 nkp = SIZE(e, 2)
215 kts = 0.0_dp
216
217 IF (sigma > 1.0e-14_dp) THEN
218 DO ik = 1, nkp
219 DO is = 1, nmo
220 SELECT CASE (method)
221 CASE (smear_fermi_dirac)
222 IF (e(is, ik) > mu) THEN
223 arg = -(e(is, ik) - mu)/sigma
224 tmp = exp(arg)
225 tmp4 = tmp + 1.0_dp
226 tmp2 = tmp/tmp4
227 tmp3 = 1.0_dp/tmp4
228 tmplog = -log(tmp4)
229 term1 = tmp2*(arg + tmplog)
230 term2 = tmp3*tmplog
231 ELSE
232 arg = (e(is, ik) - mu)/sigma
233 tmp = exp(arg)
234 tmp4 = tmp + 1.0_dp
235 tmp2 = 1.0_dp/tmp4
236 tmp3 = tmp/tmp4
237 tmplog = -log(tmp4)
238 term1 = tmp2*tmplog
239 term2 = tmp3*(arg + tmplog)
240 END IF
241 f(is, ik) = maxocc*tmp2
242 kts = kts + sigma*maxocc*(term1 + term2)*wk(ik)
243
244 CASE (smear_gaussian)
245 x = (e(is, ik) - mu)/sigma
246 expx2 = exp(-x*x)
247 f(is, ik) = maxocc*0.5_dp*erfc(x)
248 kts = kts - (sigma/(2.0_dp*rootpi))*maxocc*expx2*wk(ik)
249
250 CASE (smear_mp)
251 x = (e(is, ik) - mu)/sigma
252 expx2 = exp(-x*x)
253 f(is, ik) = maxocc*(0.5_dp*erfc(x) - x/(2.0_dp*rootpi)*expx2)
254 kts = kts + (sigma/(4.0_dp*rootpi))*maxocc*(2.0_dp*x*x - 1.0_dp)*expx2*wk(ik)
255
256 CASE (smear_mv)
257 x = (e(is, ik) - mu)/sigma
258 u = x + sqrthalf
259 expu2 = exp(-u*u)
260 f(is, ik) = maxocc*(0.5_dp*erfc(u) + expu2/(sqrt2*rootpi))
261 kts = kts - (sigma/(sqrt2*rootpi))*maxocc*u*expu2*wk(ik)
262
263 CASE DEFAULT
264 cpabort("Smear2: unknown smearing method")
265 END SELECT
266 END DO
267 END DO
268 ELSE
269 ! Zero-width limit: step function
270 DO ik = 1, nkp
271 DO is = 1, nmo
272 IF (e(is, ik) <= mu) THEN
273 f(is, ik) = maxocc
274 ELSE
275 f(is, ik) = 0.0_dp
276 END IF
277 END DO
278 END DO
279 END IF
280
281 nel = 0.0_dp
282 DO ik = 1, nkp
283 nel = nel + accurate_sum(f(1:nmo, ik))*wk(ik)
284 END DO
285
286 END SUBROUTINE smear2
287
288! **************************************************************************************************
289!> \brief Bisection search for the chemical potential mu such that the total
290!> electron count equals N, for a given smearing method (Gamma point).
291!> Brackets mu by expanding outward from [min(e), max(e)] in steps
292!> of sigma, then bisects to machine precision.
293!>
294!> For MP-1 and MV: the occupation function is non-monotonic, so it's
295!> possible that pure bisection find a spurious root.
296!> We first bisect with Gaussian smearing to get a reliable initial mu,
297!> then refine with Newton's method using the actual method's dN/dmu.
298!> (dos Santos & Marzari, PRB 2023)
299!>
300!> \param f occupations (output)
301!> \param mu chemical potential found by bisection (output)
302!> \param kTS smearing correction (output)
303!> \param e eigenvalues (input)
304!> \param N target number of electrons (input)
305!> \param sigma smearing width (input)
306!> \param maxocc maximum occupation (input)
307!> \param method smearing method selector (input)
308!> \param estate excited state index for core-level spectroscopy (optional)
309!> \param festate occupation of the excited state (optional)
310! **************************************************************************************************
311 SUBROUTINE smearfixed(f, mu, kTS, e, N, sigma, maxocc, method, estate, festate)
312
313 REAL(kind=dp), INTENT(OUT) :: f(:), mu, kts
314 REAL(kind=dp), INTENT(IN) :: e(:), n, sigma, maxocc
315 INTEGER, INTENT(IN) :: method
316 INTEGER, INTENT(IN), OPTIONAL :: estate
317 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
318
319 INTEGER :: iter, my_estate, nstate
320 REAL(kind=dp) :: gsum, mu_max, mu_min, mu_now, &
321 my_festate, n_now, n_tmp
322 REAL(kind=dp), ALLOCATABLE :: gvec(:)
323
324 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
325 my_estate = estate
326 my_festate = festate
327 ELSE
328 my_estate = nint(maxocc)
329 my_festate = my_estate
330 END IF
331
332 nstate = SIZE(e)
333
334 CALL cite_smearing(method)
335
336 SELECT CASE (method)
337
338 ! Non-monotonic methods: Gaussian bisection + Newton refinement
339 CASE (smear_mp, smear_mv)
340 ! Step 1: Gaussian bisection for a reliable initial mu
341 mu_min = minval(e)
342 iter = 0
343 DO
344 iter = iter + 1
345 CALL smearocc(f, n_tmp, kts, e, mu_min, sigma, maxocc, smear_gaussian, my_estate, my_festate)
346 IF (n_tmp > n .OR. iter > 20) THEN
347 mu_min = mu_min - sigma
348 ELSE
349 EXIT
350 END IF
351 END DO
352
353 mu_max = maxval(e)
354 iter = 0
355 DO
356 iter = iter + 1
357 CALL smearocc(f, n_tmp, kts, e, mu_max, sigma, maxocc, smear_gaussian, my_estate, my_festate)
358 IF (n_tmp < n .OR. iter > 20) THEN
359 mu_max = mu_max + sigma
360 ELSE
361 EXIT
362 END IF
363 END DO
364
365 iter = 0
366 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
367 iter = iter + 1
368 mu_now = (mu_max + mu_min)/2.0_dp
369 CALL smearocc(f, n_now, kts, e, mu_now, sigma, maxocc, smear_gaussian, my_estate, my_festate)
370 IF (n_now <= n) THEN
371 mu_min = mu_now
372 ELSE
373 mu_max = mu_now
374 END IF
375 IF (iter > bisect_max_iter) EXIT
376 END DO
377 mu = (mu_max + mu_min)/2.0_dp
378
379 ! Step 2: Newton refinement with the actual method
380 ALLOCATE (gvec(nstate))
381 DO iter = 1, newton_max_iter
382 CALL smearocc(f, n_now, kts, e, mu, sigma, maxocc, method, my_estate, my_festate)
383 IF (abs(n_now - n) < n*1.0e-12_dp) EXIT
384 CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, nstate, method, my_estate, my_festate)
385 gsum = accurate_sum(gvec)
386 IF (abs(gsum) < epsilon(gsum)) EXIT
387 mu = mu + (n - n_now)/gsum
388 END DO
389 DEALLOCATE (gvec)
390
391 ! Final evaluation
392 CALL smearocc(f, n_now, kts, e, mu, sigma, maxocc, method, my_estate, my_festate)
393
394 ! Monotonic methods (FD, Gaussian): pure bisection
395 CASE DEFAULT
396 mu_min = minval(e)
397 iter = 0
398 DO
399 iter = iter + 1
400 CALL smearocc(f, n_tmp, kts, e, mu_min, sigma, maxocc, method, my_estate, my_festate)
401 IF (n_tmp > n .OR. iter > 20) THEN
402 mu_min = mu_min - sigma
403 ELSE
404 EXIT
405 END IF
406 END DO
407
408 mu_max = maxval(e)
409 iter = 0
410 DO
411 iter = iter + 1
412 CALL smearocc(f, n_tmp, kts, e, mu_max, sigma, maxocc, method, my_estate, my_festate)
413 IF (n_tmp < n .OR. iter > 20) THEN
414 mu_max = mu_max + sigma
415 ELSE
416 EXIT
417 END IF
418 END DO
419
420 iter = 0
421 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
422 iter = iter + 1
423 mu_now = (mu_max + mu_min)/2.0_dp
424 CALL smearocc(f, n_now, kts, e, mu_now, sigma, maxocc, method, my_estate, my_festate)
425 IF (n_now <= n) THEN
426 mu_min = mu_now
427 ELSE
428 mu_max = mu_now
429 END IF
430 IF (iter > bisect_max_iter) THEN
431 cpwarn("SmearFixed: maximum bisection iterations reached")
432 EXIT
433 END IF
434 END DO
435
436 mu = (mu_max + mu_min)/2.0_dp
437 CALL smearocc(f, n_now, kts, e, mu, sigma, maxocc, method, my_estate, my_festate)
438
439 END SELECT
440
441 END SUBROUTINE smearfixed
442
443! **************************************************************************************************
444!> \brief Bisection search for mu given a target electron count (k-point case,
445!> single spin channel or spin-degenerate).
446!> Initial bracket width is max(10*sigma, 0.5) for Gaussian/MP/MV,
447!> or sigma*ln[(1-eps)/eps] for Fermi-Dirac, reflecting the different
448!> tail decay rates.
449!>
450!> For MP-1 and MV: Gaussian bisection + Newton refinement
451!> (dos Santos & Marzari, PRB 2023).
452!>
453!> \param f occupations (nmo x nkp, output)
454!> \param mu chemical potential (output)
455!> \param kTS smearing correction (output)
456!> \param e eigenvalues (nmo x nkp, input)
457!> \param nel target number of electrons (input)
458!> \param wk k-point weights (input)
459!> \param sigma smearing width (input)
460!> \param maxocc maximum occupation (input)
461!> \param method smearing method selector (input)
462! **************************************************************************************************
463 SUBROUTINE smearkp(f, mu, kTS, e, nel, wk, sigma, maxocc, method)
464
465 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: f
466 REAL(kind=dp), INTENT(OUT) :: mu, kts
467 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: e
468 REAL(kind=dp), INTENT(IN) :: nel
469 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wk
470 REAL(kind=dp), INTENT(IN) :: sigma, maxocc
471 INTEGER, INTENT(IN) :: method
472
473 REAL(kind=dp), PARAMETER :: epsocc = 1.0e-12_dp
474
475 INTEGER :: bisect_method, ik, is, iter, nkp, nmo
476 REAL(kind=dp) :: de, dndmu, expu2, expx2, mu_max, mu_min, &
477 n_now, u, x
478
479 nmo = SIZE(e, 1)
480 nkp = SIZE(e, 2)
481
482 CALL cite_smearing(method)
483
484 ! Choose bisection method: Gaussian for MP/MV, actual method for FD/Gaussian
485 SELECT CASE (method)
486 CASE (smear_mp, smear_mv)
487 bisect_method = smear_gaussian
488 CASE DEFAULT
489 bisect_method = method
490 END SELECT
491
492 ! Initial bracket
493 SELECT CASE (bisect_method)
494 CASE (smear_fermi_dirac)
495 de = sigma*log((1.0_dp - epsocc)/epsocc)
496 CASE DEFAULT
497 de = 10.0_dp*sigma
498 END SELECT
499 de = max(de, 0.5_dp)
500
501 ! Bisection with bisect_method
502 mu_min = minval(e) - de
503 mu_max = maxval(e) + de
504 iter = 0
505 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
506 iter = iter + 1
507 mu = (mu_max + mu_min)/2.0_dp
508 CALL smear2(f, n_now, kts, e, mu, wk, sigma, maxocc, bisect_method)
509
510 IF (abs(n_now - nel) < nel*epsocc) EXIT
511
512 IF (n_now <= nel) THEN
513 mu_min = mu
514 ELSE
515 mu_max = mu
516 END IF
517
518 IF (iter > bisect_max_iter) THEN
519 cpwarn("Smearkp: maximum bisection iterations reached")
520 EXIT
521 END IF
522 END DO
523 mu = (mu_max + mu_min)/2.0_dp
524
525 ! Newton refinement for non-monotonic methods
526 SELECT CASE (method)
527 CASE (smear_mp, smear_mv)
528 DO iter = 1, newton_max_iter
529 CALL smear2(f, n_now, kts, e, mu, wk, sigma, maxocc, method)
530 IF (abs(n_now - nel) < nel*epsocc) EXIT
531
532 ! Compute dN/dmu = sum_{ik} wk * g_i(k) inline
533 dndmu = 0.0_dp
534 DO ik = 1, nkp
535 DO is = 1, nmo
536 x = (e(is, ik) - mu)/sigma
537 SELECT CASE (method)
538 CASE (smear_mp)
539 expx2 = exp(-x*x)
540 dndmu = dndmu + maxocc*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2*wk(ik)
541 CASE (smear_mv)
542 u = x + sqrthalf
543 expu2 = exp(-u*u)
544 dndmu = dndmu + maxocc*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2*wk(ik)
545 END SELECT
546 END DO
547 END DO
548
549 IF (abs(dndmu) < epsilon(dndmu)) EXIT
550 mu = mu + (nel - n_now)/dndmu
551 END DO
552 END SELECT
553
554 ! Final evaluation with the actual method
555 CALL smear2(f, n_now, kts, e, mu, wk, sigma, maxocc, method)
556
557 END SUBROUTINE smearkp
558
559! **************************************************************************************************
560!> \brief Bisection search for mu (k-point, spin-polarised with a shared
561!> chemical potential across both spin channels).
562!> Asserts that the third dimension of f and e is exactly 2.
563!>
564!> For MP-1 and MV: Gaussian bisection + Newton refinement.
565!>
566!> \param f occupations (nmo x nkp x 2, output)
567!> \param mu chemical potential (output)
568!> \param kTS smearing correction (output)
569!> \param e eigenvalues (nmo x nkp x 2, input)
570!> \param nel target total number of electrons (input)
571!> \param wk k-point weights (input)
572!> \param sigma smearing width (input)
573!> \param method smearing method selector (input)
574! **************************************************************************************************
575 SUBROUTINE smearkp2(f, mu, kTS, e, nel, wk, sigma, method)
576
577 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: f
578 REAL(kind=dp), INTENT(OUT) :: mu, kts
579 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: e
580 REAL(kind=dp), INTENT(IN) :: nel
581 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wk
582 REAL(kind=dp), INTENT(IN) :: sigma
583 INTEGER, INTENT(IN) :: method
584
585 REAL(kind=dp), PARAMETER :: epsocc = 1.0e-12_dp
586
587 INTEGER :: bisect_method, ik, is, ispin, iter, nkp, &
588 nmo
589 REAL(kind=dp) :: de, dndmu, expu2, expx2, ktsa, ktsb, &
590 mu_max, mu_min, n_now, na, nb, u, x
591
592 cpassert(SIZE(f, 3) == 2 .AND. SIZE(e, 3) == 2)
593
594 nmo = SIZE(e, 1)
595 nkp = SIZE(e, 2)
596
597 CALL cite_smearing(method)
598
599 SELECT CASE (method)
600 CASE (smear_mp, smear_mv)
601 bisect_method = smear_gaussian
602 CASE DEFAULT
603 bisect_method = method
604 END SELECT
605
606 SELECT CASE (bisect_method)
607 CASE (smear_fermi_dirac)
608 de = sigma*log((1.0_dp - epsocc)/epsocc)
609 CASE DEFAULT
610 de = 10.0_dp*sigma
611 END SELECT
612 de = max(de, 0.5_dp)
613
614 ! Bisection with bisect_method
615 mu_min = minval(e) - de
616 mu_max = maxval(e) + de
617 iter = 0
618 DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
619 iter = iter + 1
620 mu = (mu_max + mu_min)/2.0_dp
621 CALL smear2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, sigma, 1.0_dp, bisect_method)
622 CALL smear2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, sigma, 1.0_dp, bisect_method)
623 n_now = na + nb
624
625 IF (abs(n_now - nel) < nel*epsocc) EXIT
626
627 IF (n_now <= nel) THEN
628 mu_min = mu
629 ELSE
630 mu_max = mu
631 END IF
632
633 IF (iter > bisect_max_iter) THEN
634 cpwarn("Smearkp2: maximum bisection iterations reached")
635 EXIT
636 END IF
637 END DO
638 mu = (mu_max + mu_min)/2.0_dp
639
640 ! Newton refinement for non-monotonic methods
641 SELECT CASE (method)
642 CASE (smear_mp, smear_mv)
643 DO iter = 1, newton_max_iter
644 CALL smear2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
645 CALL smear2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
646 n_now = na + nb
647 IF (abs(n_now - nel) < nel*epsocc) EXIT
648
649 ! dN/dmu across both spin channels (maxocc=1 per spin)
650 dndmu = 0.0_dp
651 DO ispin = 1, 2
652 DO ik = 1, nkp
653 DO is = 1, nmo
654 x = (e(is, ik, ispin) - mu)/sigma
655 SELECT CASE (method)
656 CASE (smear_mp)
657 expx2 = exp(-x*x)
658 dndmu = dndmu + (3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2*wk(ik)
659 CASE (smear_mv)
660 u = x + sqrthalf
661 expu2 = exp(-u*u)
662 dndmu = dndmu + (2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2*wk(ik)
663 END SELECT
664 END DO
665 END DO
666 END DO
667
668 IF (abs(dndmu) < epsilon(dndmu)) EXIT
669 mu = mu + (nel - n_now)/dndmu
670 END DO
671 END SELECT
672
673 ! Final evaluation with the actual method
674 CALL smear2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
675 CALL smear2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
676 kts = ktsa + ktsb
677
678 END SUBROUTINE smearkp2
679
680! **************************************************************************************************
681!> \brief Computes the smearing weight vector g_i = -df_i/de_i with mu held
682!> fixed (module-private helper for the analytical Jacobian routines).
683!>
684!> Fermi-Dirac: g_i = occ * f_norm * (1 - f_norm) / sigma
685!> where f_norm = f_i/occ_i (overflow-safe, uses
686!> pre-computed f rather than re-evaluating exp)
687!> Gaussian: g_i = occ / (sigma*sqrt(pi)) * exp(-x^2)
688!> MP-1: g_i = occ * (3 - 2*x^2) / (2*sigma*sqrt(pi)) * exp(-x^2)
689!> MV: g_i = occ * (2 + sqrt(2)*x) / (sigma*sqrt(pi)) * exp(-u^2)
690!>
691!> Note: g_i can be negative for MP-1 (|x| > sqrt(3/2)) and MV
692!> (x < -sqrt(2)). This does not break bisection (N(mu) is still
693!> monotone) but the Jacobian routines must guard against G = sum(g_i)
694!> being near zero.
695!>
696!> \param gvec weight vector (Nstate, output)
697!> \param f occupations from a prior SmearOcc/SmearFixed call (input)
698!> \param e eigenvalues (input)
699!> \param mu chemical potential (input)
700!> \param sigma smearing width (input)
701!> \param maxocc maximum occupation (input)
702!> \param Nstate number of states (input)
703!> \param method smearing method selector (input)
704!> \param estate excited state index (optional)
705!> \param festate occupation of the excited state (optional)
706! **************************************************************************************************
707 SUBROUTINE compute_gvec(gvec, f, e, mu, sigma, maxocc, Nstate, method, estate, festate)
708
709 REAL(kind=dp), INTENT(OUT) :: gvec(:)
710 REAL(kind=dp), INTENT(IN) :: f(:), e(:), mu, sigma, maxocc
711 INTEGER, INTENT(IN) :: nstate, method
712 INTEGER, INTENT(IN), OPTIONAL :: estate
713 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
714
715 INTEGER :: i
716 REAL(kind=dp) :: expu2, expx2, fi_norm, occ_i, u, x
717
718 DO i = 1, nstate
719 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
720 IF (i == estate) THEN
721 occ_i = festate
722 ELSE
723 occ_i = maxocc
724 END IF
725 ELSE
726 occ_i = maxocc
727 END IF
728
729 IF (occ_i < epsilon(occ_i)) THEN
730 gvec(i) = 0.0_dp
731 cycle
732 END IF
733
734 x = (e(i) - mu)/sigma
735
736 SELECT CASE (method)
737 CASE (smear_fermi_dirac)
738 fi_norm = f(i)/occ_i
739 gvec(i) = occ_i*fi_norm*(1.0_dp - fi_norm)/sigma
740
741 CASE (smear_gaussian)
742 expx2 = exp(-x*x)
743 gvec(i) = occ_i/(sigma*rootpi)*expx2
744
745 CASE (smear_mp)
746 expx2 = exp(-x*x)
747 gvec(i) = occ_i*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2
748
749 CASE (smear_mv)
750 u = x + sqrthalf
751 expu2 = exp(-u*u)
752 gvec(i) = occ_i*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2
753
754 END SELECT
755 END DO
756
757 END SUBROUTINE compute_gvec
758
759! **************************************************************************************************
760!> \brief Analytical Jacobian df_i/de_j for any smearing method under the
761!> electron-number constraint sum(f) = N.
762!>
763!> Differentiating f_i(e, mu(e)) where mu is implicitly defined by
764!> the constraint yields:
765!>
766!> df_i/de_j = -delta_{ij} * g_i + g_i * g_j / G
767!>
768!> where g_i = -df_i/de_i (mu fixed) and G = sum(g_i).
769!> This is a diagonal matrix plus a symmetric rank-1 update.
770!> Building it costs O(N) for g, plus O(N^2) for the outer product.
771!>
772!> Replaces the original numerical finite-difference FermiFixedDeriv
773!> which required 2N bisection solves. Exact to machine precision
774!> for all four methods.
775!>
776!> \param dfde Jacobian matrix dfde(i,j) = df_i/de_j (Nstate x Nstate, output)
777!> \param f occupations (output)
778!> \param mu chemical potential (output)
779!> \param kTS smearing correction (output)
780!> \param e eigenvalues (input)
781!> \param N target number of electrons (input)
782!> \param sigma smearing width (input)
783!> \param maxocc maximum occupation (input)
784!> \param method smearing method selector (input)
785!> \param estate excited state index (optional)
786!> \param festate occupation of the excited state (optional)
787! **************************************************************************************************
788 SUBROUTINE smearfixedderiv(dfde, f, mu, kTS, e, N, sigma, maxocc, method, estate, festate)
789
790 REAL(kind=dp), INTENT(OUT) :: dfde(:, :), f(:), mu, kts
791 REAL(kind=dp), INTENT(IN) :: e(:), n, sigma, maxocc
792 INTEGER, INTENT(IN) :: method
793 INTEGER, INTENT(IN), OPTIONAL :: estate
794 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
795
796 CHARACTER(len=*), PARAMETER :: routinen = 'SmearFixedDeriv'
797
798 INTEGER :: handle, i, j, nstate
799 REAL(kind=dp) :: gsum
800 REAL(kind=dp), ALLOCATABLE :: gvec(:)
801
802 CALL timeset(routinen, handle)
803
804 ! Step 1: find mu and f
805 CALL smearfixed(f, mu, kts, e, n, sigma, maxocc, method, estate, festate)
806
807 ! Step 2: build g vector
808 nstate = SIZE(e)
809 ALLOCATE (gvec(nstate))
810 CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, nstate, method, estate, festate)
811 gsum = accurate_sum(gvec)
812
813 ! Step 3: assemble dfde(i,j) = -delta_{ij}*g_i + g_i*g_j/G
814 IF (abs(gsum) > epsilon(gsum)) THEN
815 DO j = 1, nstate
816 DO i = 1, nstate
817 dfde(i, j) = gvec(i)*gvec(j)/gsum
818 END DO
819 dfde(j, j) = dfde(j, j) - gvec(j)
820 END DO
821 ELSE
822 dfde(:, :) = 0.0_dp
823 END IF
824
825 DEALLOCATE (gvec)
826 CALL timestop(handle)
827
828 END SUBROUTINE smearfixedderiv
829
830! **************************************************************************************************
831!> \brief Apply TRANSPOSE(df/de) to a vector WITHOUT forming the full N x N
832!> Jacobian. O(N) time and O(N) memory for all four methods.
833!>
834!> Exploiting the rank-1 structure of the constrained Jacobian:
835!>
836!> [J^T v]_j = g_j * (g . v / G - v_j)
837!>
838!> This replaces the pattern used in qs_mo_occupation:
839!> ALLOCATE(dfde(nmo,nmo))
840!> CALL SmearFixedDeriv(dfde, ...)
841!> RESULT = MATMUL(TRANSPOSE(dfde), v)
842!> DEALLOCATE(dfde)
843!> turning O(N^2) storage + O(N^2) MATMUL into O(N) throughout.
844!>
845!> Currently the sole caller (qs_ot_scf do_ener) is dead code, but
846!> this routine is ready for when it is enabled.
847!>
848!> \param RESULT output vector = TRANSPOSE(df/de) * v (Nstate, output)
849!> \param f occupations (output)
850!> \param mu chemical potential (output)
851!> \param kTS smearing correction (output)
852!> \param e eigenvalues (input)
853!> \param N_el target number of electrons (input)
854!> \param sigma smearing width (input)
855!> \param maxocc maximum occupation (input)
856!> \param method smearing method selector (input)
857!> \param v input vector to multiply (Nstate, input)
858!> \param estate excited state index (optional)
859!> \param festate occupation of the excited state (optional)
860! **************************************************************************************************
861 SUBROUTINE smearfixedderivmv(RESULT, f, mu, kTS, e, N_el, sigma, maxocc, method, v, estate, festate)
862
863 REAL(kind=dp), INTENT(OUT) :: result(:), f(:), mu, kts
864 REAL(kind=dp), INTENT(IN) :: e(:), n_el, sigma, maxocc
865 INTEGER, INTENT(IN) :: method
866 REAL(kind=dp), INTENT(IN) :: v(:)
867 INTEGER, INTENT(IN), OPTIONAL :: estate
868 REAL(kind=dp), INTENT(IN), OPTIONAL :: festate
869
870 CHARACTER(len=*), PARAMETER :: routinen = 'SmearFixedDerivMV'
871
872 INTEGER :: handle, i, nstate
873 REAL(kind=dp) :: gdotv, gsum
874 REAL(kind=dp), ALLOCATABLE :: gvec(:)
875
876 CALL timeset(routinen, handle)
877
878 ! Step 1: find mu and f
879 CALL smearfixed(f, mu, kts, e, n_el, sigma, maxocc, method, estate, festate)
880
881 ! Step 2: build g vector
882 nstate = SIZE(e)
883 ALLOCATE (gvec(nstate))
884 CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, nstate, method, estate, festate)
885 gsum = accurate_sum(gvec)
886
887 ! Step 3: RESULT_j = g_j * (g.v / G - v_j)
888 IF (abs(gsum) > epsilon(gsum)) THEN
889 gdotv = 0.0_dp
890 DO i = 1, nstate
891 gdotv = gdotv + gvec(i)*v(i)
892 END DO
893 DO i = 1, nstate
894 result(i) = gvec(i)*(gdotv/gsum - v(i))
895 END DO
896 ELSE
897 result(:) = 0.0_dp
898 END IF
899
900 DEALLOCATE (gvec)
901 CALL timestop(handle)
902
903 END SUBROUTINE smearfixedderivmv
904
905END MODULE smearing_utils
kahan_sum::accurate_sum
Definition kahan_sum.F:41
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::fuho1983
integer, save, public fuho1983
Definition bibliography.F:36
bibliography::mermin1965
integer, save, public mermin1965
Definition bibliography.F:36
bibliography::marzari1999
integer, save, public marzari1999
Definition bibliography.F:36
bibliography::dossantos2023
integer, save, public dossantos2023
Definition bibliography.F:36
bibliography::methfesselpaxton1989
integer, save, public methfesselpaxton1989
Definition bibliography.F:36
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:910
input_constants::smear_gaussian
integer, parameter, public smear_gaussian
Definition input_constants.F:910
input_constants::smear_mv
integer, parameter, public smear_mv
Definition input_constants.F:910
input_constants::smear_mp
integer, parameter, public smear_mp
Definition input_constants.F:910
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:464
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:576
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:789
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:105
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:312
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:862