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