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