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