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