Line data Source code
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 : ! **************************************************************************************************
24 : MODULE 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 :
56 : CONTAINS
57 : ! **************************************************************************************************
58 : !> \brief Citation of Smearing methods
59 : !> \param method ...
60 : ! **************************************************************************************************
61 25462 : SUBROUTINE cite_smearing(method)
62 : INTEGER, INTENT(IN) :: method
63 :
64 44440 : SELECT CASE (method)
65 : CASE (smear_fermi_dirac)
66 18978 : CALL cite_reference(Mermin1965)
67 : CASE (smear_gaussian)
68 6428 : CALL cite_reference(FuHo1983)
69 : CASE (smear_mp)
70 28 : CALL cite_reference(FuHo1983)
71 28 : CALL cite_reference(MethfesselPaxton1989)
72 28 : CALL cite_reference(dosSantos2023)
73 : CASE (smear_mv)
74 28 : CALL cite_reference(FuHo1983)
75 28 : CALL cite_reference(Marzari1999)
76 25490 : CALL cite_reference(dosSantos2023)
77 : END SELECT
78 25462 : 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 1112896 : 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 1112896 : Nstate = SIZE(e)
118 1112896 : kTS = 0.0_dp
119 :
120 22861324 : DO i = 1, Nstate
121 21748428 : IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
122 21714364 : IF (i == estate) THEN
123 4012 : occupation = festate
124 : ELSE
125 21710352 : occupation = maxocc
126 : END IF
127 : ELSE
128 34064 : occupation = maxocc
129 : END IF
130 :
131 1112896 : SELECT CASE (method)
132 : CASE (smear_fermi_dirac)
133 19506100 : IF (e(i) > mu) THEN
134 11236829 : arg = -(e(i) - mu)/sigma
135 11236829 : tmp = EXP(arg)
136 11236829 : tmp4 = tmp + 1.0_dp
137 11236829 : tmp2 = tmp/tmp4
138 11236829 : tmp3 = 1.0_dp/tmp4
139 11236829 : tmplog = -LOG(tmp4)
140 11236829 : term1 = tmp2*(arg + tmplog)
141 11236829 : term2 = tmp3*tmplog
142 : ELSE
143 8269271 : arg = (e(i) - mu)/sigma
144 8269271 : tmp = EXP(arg)
145 8269271 : tmp4 = tmp + 1.0_dp
146 8269271 : tmp2 = 1.0_dp/tmp4
147 8269271 : tmp3 = tmp/tmp4
148 8269271 : tmplog = -LOG(tmp4)
149 8269271 : term1 = tmp2*tmplog
150 8269271 : term2 = tmp3*(arg + tmplog)
151 : END IF
152 19506100 : f(i) = occupation*tmp2
153 19506100 : kTS = kTS + sigma*occupation*(term1 + term2)
154 :
155 : CASE (smear_gaussian)
156 2242328 : x = (e(i) - mu)/sigma
157 2242328 : expx2 = EXP(-x*x)
158 2242328 : f(i) = occupation*0.5_dp*ERFC(x)
159 2242328 : kTS = kTS - (sigma/(2.0_dp*rootpi))*occupation*expx2
160 :
161 : CASE (smear_mp)
162 0 : x = (e(i) - mu)/sigma
163 0 : expx2 = EXP(-x*x)
164 0 : f(i) = occupation*(0.5_dp*ERFC(x) - x/(2.0_dp*rootpi)*expx2)
165 0 : kTS = kTS + (sigma/(4.0_dp*rootpi))*occupation*(2.0_dp*x*x - 1.0_dp)*expx2
166 :
167 : CASE (smear_mv)
168 0 : x = (e(i) - mu)/sigma
169 0 : u = x + sqrthalf
170 0 : expu2 = EXP(-u*u)
171 0 : f(i) = occupation*(0.5_dp*ERFC(u) + expu2/(sqrt2*rootpi))
172 0 : kTS = kTS - (sigma/(sqrt2*rootpi))*occupation*u*expu2
173 :
174 : CASE DEFAULT
175 21748428 : CPABORT("SmearOcc: unknown smearing method")
176 : END SELECT
177 : END DO
178 :
179 1112896 : N = accurate_sum(f)
180 :
181 1112896 : 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 44342 : 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 44342 : nmo = SIZE(e, 1)
214 44342 : nkp = SIZE(e, 2)
215 44342 : kTS = 0.0_dp
216 :
217 44342 : IF (sigma > 1.0e-14_dp) THEN
218 295888 : DO ik = 1, nkp
219 10189684 : DO is = 1, nmo
220 267112 : SELECT CASE (method)
221 : CASE (smear_fermi_dirac)
222 8939920 : IF (e(is, ik) > mu) THEN
223 4912626 : arg = -(e(is, ik) - mu)/sigma
224 4912626 : tmp = EXP(arg)
225 4912626 : tmp4 = tmp + 1.0_dp
226 4912626 : tmp2 = tmp/tmp4
227 4912626 : tmp3 = 1.0_dp/tmp4
228 4912626 : tmplog = -LOG(tmp4)
229 4912626 : term1 = tmp2*(arg + tmplog)
230 4912626 : term2 = tmp3*tmplog
231 : ELSE
232 4027294 : arg = (e(is, ik) - mu)/sigma
233 4027294 : tmp = EXP(arg)
234 4027294 : tmp4 = tmp + 1.0_dp
235 4027294 : tmp2 = 1.0_dp/tmp4
236 4027294 : tmp3 = tmp/tmp4
237 4027294 : tmplog = -LOG(tmp4)
238 4027294 : term1 = tmp2*tmplog
239 4027294 : term2 = tmp3*(arg + tmplog)
240 : END IF
241 8939920 : f(is, ik) = maxocc*tmp2
242 8939920 : kTS = kTS + sigma*maxocc*(term1 + term2)*wk(ik)
243 :
244 : CASE (smear_gaussian)
245 855316 : x = (e(is, ik) - mu)/sigma
246 855316 : expx2 = EXP(-x*x)
247 855316 : f(is, ik) = maxocc*0.5_dp*ERFC(x)
248 855316 : kTS = kTS - (sigma/(2.0_dp*rootpi))*maxocc*expx2*wk(ik)
249 :
250 : CASE (smear_mp)
251 44800 : x = (e(is, ik) - mu)/sigma
252 44800 : expx2 = EXP(-x*x)
253 44800 : f(is, ik) = maxocc*(0.5_dp*ERFC(x) - x/(2.0_dp*rootpi)*expx2)
254 44800 : 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 53760 : x = (e(is, ik) - mu)/sigma
258 53760 : u = x + sqrthalf
259 53760 : expu2 = EXP(-u*u)
260 53760 : f(is, ik) = maxocc*(0.5_dp*ERFC(u) + expu2/(sqrt2*rootpi))
261 53760 : kTS = kTS - (sigma/(sqrt2*rootpi))*maxocc*u*expu2*wk(ik)
262 :
263 : CASE DEFAULT
264 9893796 : CPABORT("Smear2: unknown smearing method")
265 : END SELECT
266 : END DO
267 : END DO
268 : ELSE
269 : ! Zero-width limit: step function
270 80530 : DO ik = 1, nkp
271 689150 : DO is = 1, nmo
272 673584 : IF (e(is, ik) <= mu) THEN
273 351784 : f(is, ik) = maxocc
274 : ELSE
275 256836 : f(is, ik) = 0.0_dp
276 : END IF
277 : END DO
278 : END DO
279 : END IF
280 :
281 44342 : nel = 0.0_dp
282 376418 : DO ik = 1, nkp
283 376418 : nel = nel + accurate_sum(f(1:nmo, ik))*wk(ik)
284 : END DO
285 :
286 44342 : 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 19808 : 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 19808 : REAL(KIND=dp), ALLOCATABLE :: gvec(:)
323 :
324 19808 : IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
325 19768 : my_estate = estate
326 19768 : my_festate = festate
327 : ELSE
328 40 : my_estate = NINT(maxocc)
329 40 : my_festate = my_estate
330 : END IF
331 :
332 19808 : Nstate = SIZE(e)
333 :
334 19808 : CALL cite_smearing(method)
335 :
336 19808 : 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 0 : mu_min = MINVAL(e)
342 0 : iter = 0
343 0 : DO
344 0 : iter = iter + 1
345 0 : CALL SmearOcc(f, N_tmp, kTS, e, mu_min, sigma, maxocc, smear_gaussian, my_estate, my_festate)
346 0 : IF (N_tmp > N .OR. iter > 20) THEN
347 0 : mu_min = mu_min - sigma
348 : ELSE
349 : EXIT
350 : END IF
351 : END DO
352 :
353 0 : mu_max = MAXVAL(e)
354 0 : iter = 0
355 0 : DO
356 0 : iter = iter + 1
357 0 : CALL SmearOcc(f, N_tmp, kTS, e, mu_max, sigma, maxocc, smear_gaussian, my_estate, my_festate)
358 0 : IF (N_tmp < N .OR. iter > 20) THEN
359 0 : mu_max = mu_max + sigma
360 : ELSE
361 : EXIT
362 : END IF
363 : END DO
364 :
365 : iter = 0
366 0 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
367 0 : iter = iter + 1
368 0 : mu_now = (mu_max + mu_min)/2.0_dp
369 0 : CALL SmearOcc(f, N_now, kTS, e, mu_now, sigma, maxocc, smear_gaussian, my_estate, my_festate)
370 0 : IF (N_now <= N) THEN
371 0 : mu_min = mu_now
372 : ELSE
373 0 : mu_max = mu_now
374 : END IF
375 0 : IF (iter > BISECT_MAX_ITER) EXIT
376 : END DO
377 0 : mu = (mu_max + mu_min)/2.0_dp
378 :
379 : ! Step 2: Newton refinement with the actual method
380 0 : ALLOCATE (gvec(Nstate))
381 0 : DO iter = 1, NEWTON_MAX_ITER
382 0 : CALL SmearOcc(f, N_now, kTS, e, mu, sigma, maxocc, method, my_estate, my_festate)
383 0 : IF (ABS(N_now - N) < N*1.0e-12_dp) EXIT
384 0 : CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, Nstate, method, my_estate, my_festate)
385 0 : Gsum = accurate_sum(gvec)
386 0 : IF (ABS(Gsum) < EPSILON(Gsum)) EXIT
387 0 : mu = mu + (N - N_now)/Gsum
388 : END DO
389 0 : DEALLOCATE (gvec)
390 :
391 : ! Final evaluation
392 0 : 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 424918 : mu_min = MINVAL(e)
397 19808 : iter = 0
398 0 : DO
399 19808 : iter = iter + 1
400 19808 : CALL SmearOcc(f, N_tmp, kTS, e, mu_min, sigma, maxocc, method, my_estate, my_festate)
401 19808 : IF (N_tmp > N .OR. iter > 20) THEN
402 0 : mu_min = mu_min - sigma
403 : ELSE
404 : EXIT
405 : END IF
406 : END DO
407 :
408 424918 : mu_max = MAXVAL(e)
409 19808 : iter = 0
410 0 : DO
411 19808 : iter = iter + 1
412 19808 : CALL SmearOcc(f, N_tmp, kTS, e, mu_max, sigma, maxocc, method, my_estate, my_festate)
413 19808 : IF (N_tmp < N .OR. iter > 20) THEN
414 0 : mu_max = mu_max + sigma
415 : ELSE
416 : EXIT
417 : END IF
418 : END DO
419 :
420 : iter = 0
421 1072904 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
422 1053096 : iter = iter + 1
423 1053096 : mu_now = (mu_max + mu_min)/2.0_dp
424 1053096 : CALL SmearOcc(f, N_now, kTS, e, mu_now, sigma, maxocc, method, my_estate, my_festate)
425 1053096 : IF (N_now <= N) THEN
426 516751 : mu_min = mu_now
427 : ELSE
428 536345 : mu_max = mu_now
429 : END IF
430 1072904 : IF (iter > BISECT_MAX_ITER) THEN
431 0 : CPWARN("SmearFixed: maximum bisection iterations reached")
432 0 : EXIT
433 : END IF
434 : END DO
435 :
436 19808 : mu = (mu_max + mu_min)/2.0_dp
437 39616 : CALL SmearOcc(f, N_now, kTS, e, mu, sigma, maxocc, method, my_estate, my_festate)
438 :
439 : END SELECT
440 :
441 19808 : 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 5626 : 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 5626 : nmo = SIZE(e, 1)
480 5626 : nkp = SIZE(e, 2)
481 :
482 5626 : CALL cite_smearing(method)
483 :
484 : ! Choose bisection method: Gaussian for MP/MV, actual method for FD/Gaussian
485 5682 : SELECT CASE (method)
486 : CASE (smear_mp, smear_mv)
487 56 : bisect_method = smear_gaussian
488 : CASE DEFAULT
489 5626 : bisect_method = method
490 : END SELECT
491 :
492 : ! Initial bracket
493 696 : SELECT CASE (bisect_method)
494 : CASE (smear_fermi_dirac)
495 696 : de = sigma*LOG((1.0_dp - epsocc)/epsocc)
496 : CASE DEFAULT
497 5626 : de = 10.0_dp*sigma
498 : END SELECT
499 5626 : de = MAX(de, 0.5_dp)
500 :
501 : ! Bisection with bisect_method
502 493066 : mu_min = MINVAL(e) - de
503 493066 : mu_max = MAXVAL(e) + de
504 5626 : iter = 0
505 36668 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
506 36668 : iter = iter + 1
507 36668 : mu = (mu_max + mu_min)/2.0_dp
508 36668 : CALL Smear2(f, N_now, kTS, e, mu, wk, sigma, maxocc, bisect_method)
509 :
510 36668 : IF (ABS(N_now - nel) < nel*epsocc) EXIT
511 :
512 31042 : IF (N_now <= nel) THEN
513 16914 : mu_min = mu
514 : ELSE
515 14128 : mu_max = mu
516 : END IF
517 :
518 36668 : IF (iter > BISECT_MAX_ITER) THEN
519 0 : CPWARN("Smearkp: maximum bisection iterations reached")
520 0 : EXIT
521 : END IF
522 : END DO
523 5626 : 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 5822 : DO iter = 1, NEWTON_MAX_ITER
529 252 : CALL Smear2(f, N_now, kTS, e, mu, wk, sigma, maxocc, method)
530 252 : 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 6468 : DO ik = 1, nkp
535 69188 : DO is = 1, nmo
536 62720 : x = (e(is, ik) - mu)/sigma
537 6272 : SELECT CASE (method)
538 : CASE (smear_mp)
539 26880 : expx2 = EXP(-x*x)
540 26880 : dNdmu = dNdmu + maxocc*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2*wk(ik)
541 : CASE (smear_mv)
542 35840 : u = x + sqrthalf
543 35840 : expu2 = EXP(-u*u)
544 62720 : dNdmu = dNdmu + maxocc*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2*wk(ik)
545 : END SELECT
546 : END DO
547 : END DO
548 :
549 196 : IF (ABS(dNdmu) < EPSILON(dNdmu)) EXIT
550 252 : mu = mu + (nel - N_now)/dNdmu
551 : END DO
552 : END SELECT
553 :
554 : ! Final evaluation with the actual method
555 5626 : CALL Smear2(f, N_now, kTS, e, mu, wk, sigma, maxocc, method)
556 :
557 5626 : 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 28 : 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 28 : CPASSERT(SIZE(f, 3) == 2 .AND. SIZE(e, 3) == 2)
593 :
594 28 : nmo = SIZE(e, 1)
595 28 : nkp = SIZE(e, 2)
596 :
597 28 : CALL cite_smearing(method)
598 :
599 28 : SELECT CASE (method)
600 : CASE (smear_mp, smear_mv)
601 0 : bisect_method = smear_gaussian
602 : CASE DEFAULT
603 28 : bisect_method = method
604 : END SELECT
605 :
606 0 : SELECT CASE (bisect_method)
607 : CASE (smear_fermi_dirac)
608 0 : de = sigma*LOG((1.0_dp - epsocc)/epsocc)
609 : CASE DEFAULT
610 28 : de = 10.0_dp*sigma
611 : END SELECT
612 28 : de = MAX(de, 0.5_dp)
613 :
614 : ! Bisection with bisect_method
615 1740 : mu_min = MINVAL(e) - de
616 1740 : mu_max = MAXVAL(e) + de
617 28 : iter = 0
618 870 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
619 870 : iter = iter + 1
620 870 : mu = (mu_max + mu_min)/2.0_dp
621 870 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, sigma, 1.0_dp, bisect_method)
622 870 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, sigma, 1.0_dp, bisect_method)
623 870 : N_now = na + nb
624 :
625 870 : IF (ABS(N_now - nel) < nel*epsocc) EXIT
626 :
627 842 : IF (N_now <= nel) THEN
628 380 : mu_min = mu
629 : ELSE
630 462 : mu_max = mu
631 : END IF
632 :
633 870 : IF (iter > BISECT_MAX_ITER) THEN
634 0 : CPWARN("Smearkp2: maximum bisection iterations reached")
635 0 : EXIT
636 : END IF
637 : END DO
638 28 : 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 28 : DO iter = 1, NEWTON_MAX_ITER
644 0 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
645 0 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
646 0 : N_now = na + nb
647 0 : 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 0 : DO ispin = 1, 2
652 0 : DO ik = 1, nkp
653 0 : DO is = 1, nmo
654 0 : x = (e(is, ik, ispin) - mu)/sigma
655 0 : SELECT CASE (method)
656 : CASE (smear_mp)
657 0 : expx2 = EXP(-x*x)
658 0 : dNdmu = dNdmu + (3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2*wk(ik)
659 : CASE (smear_mv)
660 0 : u = x + sqrthalf
661 0 : expu2 = EXP(-u*u)
662 0 : 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 0 : IF (ABS(dNdmu) < EPSILON(dNdmu)) EXIT
669 0 : mu = mu + (nel - N_now)/dNdmu
670 : END DO
671 : END SELECT
672 :
673 : ! Final evaluation with the actual method
674 28 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
675 28 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
676 28 : kTS = kTSa + kTSb
677 :
678 28 : 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 0 : 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 0 : DO i = 1, Nstate
719 0 : IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
720 0 : IF (i == estate) THEN
721 0 : occ_i = festate
722 : ELSE
723 0 : occ_i = maxocc
724 : END IF
725 : ELSE
726 0 : occ_i = maxocc
727 : END IF
728 :
729 0 : IF (occ_i < EPSILON(occ_i)) THEN
730 0 : gvec(i) = 0.0_dp
731 0 : CYCLE
732 : END IF
733 :
734 0 : x = (e(i) - mu)/sigma
735 :
736 0 : SELECT CASE (method)
737 : CASE (smear_fermi_dirac)
738 0 : fi_norm = f(i)/occ_i
739 0 : gvec(i) = occ_i*fi_norm*(1.0_dp - fi_norm)/sigma
740 :
741 : CASE (smear_gaussian)
742 0 : expx2 = EXP(-x*x)
743 0 : gvec(i) = occ_i/(sigma*rootpi)*expx2
744 :
745 : CASE (smear_mp)
746 0 : expx2 = EXP(-x*x)
747 0 : gvec(i) = occ_i*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2
748 :
749 : CASE (smear_mv)
750 0 : u = x + sqrthalf
751 0 : expu2 = EXP(-u*u)
752 0 : gvec(i) = occ_i*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2
753 :
754 : END SELECT
755 : END DO
756 :
757 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE :: gvec(:)
801 :
802 0 : CALL timeset(routineN, handle)
803 :
804 : ! Step 1: find mu and f
805 0 : CALL SmearFixed(f, mu, kTS, e, N, sigma, maxocc, method, estate, festate)
806 :
807 : ! Step 2: build g vector
808 0 : Nstate = SIZE(e)
809 0 : ALLOCATE (gvec(Nstate))
810 0 : CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, Nstate, method, estate, festate)
811 0 : Gsum = accurate_sum(gvec)
812 :
813 : ! Step 3: assemble dfde(i,j) = -delta_{ij}*g_i + g_i*g_j/G
814 0 : IF (ABS(Gsum) > EPSILON(Gsum)) THEN
815 0 : DO j = 1, Nstate
816 0 : DO i = 1, Nstate
817 0 : dfde(i, j) = gvec(i)*gvec(j)/Gsum
818 : END DO
819 0 : dfde(j, j) = dfde(j, j) - gvec(j)
820 : END DO
821 : ELSE
822 0 : dfde(:, :) = 0.0_dp
823 : END IF
824 :
825 0 : DEALLOCATE (gvec)
826 0 : CALL timestop(handle)
827 :
828 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE :: gvec(:)
875 :
876 0 : CALL timeset(routineN, handle)
877 :
878 : ! Step 1: find mu and f
879 0 : CALL SmearFixed(f, mu, kTS, e, N_el, sigma, maxocc, method, estate, festate)
880 :
881 : ! Step 2: build g vector
882 0 : Nstate = SIZE(e)
883 0 : ALLOCATE (gvec(Nstate))
884 0 : CALL compute_gvec(gvec, f, e, mu, sigma, maxocc, Nstate, method, estate, festate)
885 0 : Gsum = accurate_sum(gvec)
886 :
887 : ! Step 3: RESULT_j = g_j * (g.v / G - v_j)
888 0 : IF (ABS(Gsum) > EPSILON(Gsum)) THEN
889 : gdotv = 0.0_dp
890 0 : DO i = 1, Nstate
891 0 : gdotv = gdotv + gvec(i)*v(i)
892 : END DO
893 0 : DO i = 1, Nstate
894 0 : RESULT(i) = gvec(i)*(gdotv/Gsum - v(i))
895 : END DO
896 : ELSE
897 0 : RESULT(:) = 0.0_dp
898 : END IF
899 :
900 0 : DEALLOCATE (gvec)
901 0 : CALL timestop(handle)
902 :
903 0 : END SUBROUTINE SmearFixedDerivMV
904 :
905 : END MODULE smearing_utils
|
