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 : 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 :
60 : CONTAINS
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 40606 : PURE LOGICAL FUNCTION smearing_separate_spin(fixed_mag_mom)
67 : REAL(KIND=dp), INTENT(IN) :: fixed_mag_mom
68 :
69 40606 : smearing_separate_spin = fixed_mag_mom >= 0.0_dp .OR. fixed_mag_mom == -1.0_dp
70 40606 : 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 41622 : 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 124866 : targets(:) = REAL(nelectron, dp)
86 41622 : IF (fixed_mag_mom >= 0.0_dp) THEN
87 840 : total = SUM(targets)
88 280 : targets(1) = 0.5_dp*(total + fixed_mag_mom)
89 280 : targets(2) = 0.5_dp*(total - fixed_mag_mom)
90 : END IF
91 41622 : END FUNCTION smearing_spin_targets
92 :
93 : ! **************************************************************************************************
94 : !> \brief Citation of Smearing methods
95 : !> \param method ...
96 : ! **************************************************************************************************
97 70926 : SUBROUTINE cite_smearing(method)
98 : INTEGER, INTENT(IN) :: method
99 :
100 112234 : SELECT CASE (method)
101 : CASE (smear_fermi_dirac)
102 41308 : CALL cite_reference(Mermin1965)
103 : CASE (smear_gaussian)
104 29334 : CALL cite_reference(FuHo1983)
105 : CASE (smear_mp)
106 142 : CALL cite_reference(FuHo1983)
107 142 : CALL cite_reference(MethfesselPaxton1989)
108 142 : CALL cite_reference(dosSantos2023)
109 : CASE (smear_mv)
110 142 : CALL cite_reference(FuHo1983)
111 142 : CALL cite_reference(Marzari1999)
112 71068 : CALL cite_reference(dosSantos2023)
113 : END SELECT
114 70926 : 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 1506476 : 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 1506476 : Nstate = SIZE(e)
154 1506476 : kTS = 0.0_dp
155 :
156 43752920 : DO i = 1, Nstate
157 42246444 : IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
158 42202196 : IF (i == estate) THEN
159 4834 : occupation = festate
160 : ELSE
161 42197362 : occupation = maxocc
162 : END IF
163 : ELSE
164 44248 : occupation = maxocc
165 : END IF
166 :
167 1506476 : SELECT CASE (method)
168 : CASE (smear_fermi_dirac)
169 39909538 : IF (e(i) > mu) THEN
170 24635239 : arg = -(e(i) - mu)/sigma
171 24635239 : tmp = EXP(arg)
172 24635239 : tmp4 = tmp + 1.0_dp
173 24635239 : tmp2 = tmp/tmp4
174 24635239 : tmp3 = 1.0_dp/tmp4
175 24635239 : tmplog = -LOG(tmp4)
176 24635239 : term1 = tmp2*(arg + tmplog)
177 24635239 : term2 = tmp3*tmplog
178 : ELSE
179 15274299 : arg = (e(i) - mu)/sigma
180 15274299 : tmp = EXP(arg)
181 15274299 : tmp4 = tmp + 1.0_dp
182 15274299 : tmp2 = 1.0_dp/tmp4
183 15274299 : tmp3 = tmp/tmp4
184 15274299 : tmplog = -LOG(tmp4)
185 15274299 : term1 = tmp2*tmplog
186 15274299 : term2 = tmp3*(arg + tmplog)
187 : END IF
188 39909538 : f(i) = occupation*tmp2
189 39909538 : kTS = kTS + sigma*occupation*(term1 + term2)
190 :
191 : CASE (smear_gaussian)
192 2330218 : x = (e(i) - mu)/sigma
193 2330218 : expx2 = EXP(-x*x)
194 2330218 : f(i) = occupation*0.5_dp*ERFC(x)
195 2330218 : kTS = kTS - (sigma/(2.0_dp*rootpi))*occupation*expx2
196 :
197 : CASE (smear_mp)
198 2128 : x = (e(i) - mu)/sigma
199 2128 : expx2 = EXP(-x*x)
200 2128 : f(i) = occupation*(0.5_dp*ERFC(x) - x/(2.0_dp*rootpi)*expx2)
201 2128 : kTS = kTS + (sigma/(4.0_dp*rootpi))*occupation*(2.0_dp*x*x - 1.0_dp)*expx2
202 :
203 : CASE (smear_mv)
204 4560 : x = (e(i) - mu)/sigma
205 4560 : u = x + sqrthalf
206 4560 : expu2 = EXP(-u*u)
207 4560 : f(i) = occupation*(0.5_dp*ERFC(u) + expu2/(sqrt2*rootpi))
208 4560 : kTS = kTS - (sigma/(sqrt2*rootpi))*occupation*u*expu2
209 :
210 : CASE DEFAULT
211 42246444 : CPABORT("SmearOcc: unknown smearing method")
212 : END SELECT
213 : END DO
214 :
215 1506476 : N = accurate_sum(f)
216 :
217 1506476 : 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 353742 : 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 353742 : nmo = SIZE(e, 1)
250 353742 : nkp = SIZE(e, 2)
251 353742 : kTS = 0.0_dp
252 :
253 353742 : IF (sigma > 1.0e-14_dp) THEN
254 1177366 : DO ik = 1, nkp
255 20880204 : DO is = 1, nmo
256 907722 : SELECT CASE (method)
257 : CASE (smear_fermi_dirac)
258 18593568 : IF (e(is, ik) > mu) THEN
259 10886776 : arg = -(e(is, ik) - mu)/sigma
260 10886776 : tmp = EXP(arg)
261 10886776 : tmp4 = tmp + 1.0_dp
262 10886776 : tmp2 = tmp/tmp4
263 10886776 : tmp3 = 1.0_dp/tmp4
264 10886776 : tmplog = -LOG(tmp4)
265 10886776 : term1 = tmp2*(arg + tmplog)
266 10886776 : term2 = tmp3*tmplog
267 : ELSE
268 7706792 : arg = (e(is, ik) - mu)/sigma
269 7706792 : tmp = EXP(arg)
270 7706792 : tmp4 = tmp + 1.0_dp
271 7706792 : tmp2 = 1.0_dp/tmp4
272 7706792 : tmp3 = tmp/tmp4
273 7706792 : tmplog = -LOG(tmp4)
274 7706792 : term1 = tmp2*tmplog
275 7706792 : term2 = tmp3*(arg + tmplog)
276 : END IF
277 18593568 : f(is, ik) = maxocc*tmp2
278 18593568 : kTS = kTS + sigma*maxocc*(term1 + term2)*wk(ik)
279 :
280 : CASE (smear_gaussian)
281 1001750 : x = (e(is, ik) - mu)/sigma
282 1001750 : expx2 = EXP(-x*x)
283 1001750 : f(is, ik) = maxocc*0.5_dp*ERFC(x)
284 1001750 : kTS = kTS - (sigma/(2.0_dp*rootpi))*maxocc*expx2*wk(ik)
285 :
286 : CASE (smear_mp)
287 50096 : x = (e(is, ik) - mu)/sigma
288 50096 : expx2 = EXP(-x*x)
289 50096 : f(is, ik) = maxocc*(0.5_dp*ERFC(x) - x/(2.0_dp*rootpi)*expx2)
290 50096 : 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 57424 : x = (e(is, ik) - mu)/sigma
294 57424 : u = x + sqrthalf
295 57424 : expu2 = EXP(-u*u)
296 57424 : f(is, ik) = maxocc*(0.5_dp*ERFC(u) + expu2/(sqrt2*rootpi))
297 57424 : kTS = kTS - (sigma/(sqrt2*rootpi))*maxocc*u*expu2*wk(ik)
298 :
299 : CASE DEFAULT
300 19702838 : CPABORT("Smear2: unknown smearing method")
301 : END SELECT
302 : END DO
303 : END DO
304 : ELSE
305 : ! Zero-width limit: step function
306 345408 : DO ik = 1, nkp
307 2799428 : DO is = 1, nmo
308 2715330 : IF (e(is, ik) <= mu) THEN
309 1571842 : f(is, ik) = maxocc
310 : ELSE
311 882178 : f(is, ik) = 0.0_dp
312 : END IF
313 : END DO
314 : END DO
315 : END IF
316 :
317 353742 : nel = 0.0_dp
318 1522774 : DO ik = 1, nkp
319 1522774 : nel = nel + accurate_sum(f(1:nmo, ik))*wk(ik)
320 : END DO
321 :
322 353742 : 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 26856 : 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 26856 : REAL(KIND=dp), ALLOCATABLE :: gvec(:)
359 :
360 26856 : IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
361 26800 : my_estate = estate
362 26800 : my_festate = festate
363 : ELSE
364 56 : my_estate = NINT(maxocc)
365 56 : my_festate = my_estate
366 : END IF
367 :
368 26856 : Nstate = SIZE(e)
369 :
370 26856 : CALL cite_smearing(method)
371 :
372 26856 : 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 160 : mu_min = MINVAL(e)
378 8 : iter = 0
379 6 : DO
380 14 : iter = iter + 1
381 14 : CALL SmearOcc(f, N_tmp, kTS, e, mu_min, sigma, maxocc, smear_gaussian, my_estate, my_festate)
382 14 : IF (N_tmp <= N) EXIT
383 6 : IF (iter > 20) THEN
384 0 : CPABORT("SmearFixed: failed to bracket lower chemical potential")
385 : END IF
386 6 : mu_min = mu_min - sigma
387 : END DO
388 :
389 160 : mu_max = MAXVAL(e)
390 8 : iter = 0
391 0 : DO
392 8 : iter = iter + 1
393 8 : CALL SmearOcc(f, N_tmp, kTS, e, mu_max, sigma, maxocc, smear_gaussian, my_estate, my_festate)
394 8 : IF (N_tmp >= N) EXIT
395 0 : IF (iter > 20) THEN
396 0 : CPABORT("SmearFixed: failed to bracket upper chemical potential")
397 : END IF
398 0 : mu_max = mu_max + sigma
399 : END DO
400 :
401 : iter = 0
402 378 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
403 370 : iter = iter + 1
404 370 : mu_now = (mu_max + mu_min)/2.0_dp
405 370 : CALL SmearOcc(f, N_now, kTS, e, mu_now, sigma, maxocc, smear_gaussian, my_estate, my_festate)
406 370 : IF (N_now <= N) THEN
407 192 : mu_min = mu_now
408 : ELSE
409 178 : mu_max = mu_now
410 : END IF
411 378 : IF (iter > BISECT_MAX_ITER) EXIT
412 : END DO
413 8 : 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 24 : ALLOCATE (gvec(Nstate))
421 8 : CALL SmearOcc(f, N_now, kTS, e, mu, sigma, maxocc, method, my_estate, my_festate)
422 8 : res_best = ABS(N_now - N)
423 8 : mu_best = mu
424 30 : DO iter = 1, NEWTON_MAX_ITER
425 30 : res_now = ABS(N_now - N)
426 30 : IF (res_now < N*1.0e-12_dp) EXIT
427 : CALL smearing_response_weight( &
428 22 : gvec, f, e, mu, sigma, maxocc, Nstate, method, my_estate, my_festate)
429 22 : Gsum = accurate_sum(gvec)
430 22 : IF (ABS(Gsum) < EPSILON(Gsum)) EXIT
431 :
432 22 : step = (N - N_now)/Gsum
433 22 : step = SIGN(MIN(ABS(step), MPMV_MAX_NEWTON_STEP*sigma), step)
434 22 : step_try = step
435 22 : DO iback = 1, NEWTON_MAX_BACKTRACK
436 22 : mu_trial = mu + step_try
437 : CALL SmearOcc(f, N_trial, kTS, e, mu_trial, sigma, maxocc, method, &
438 22 : my_estate, my_festate)
439 22 : res_trial = ABS(N_trial - N)
440 22 : IF (res_trial < res_now) THEN
441 22 : mu = mu_trial
442 22 : N_now = N_trial
443 22 : IF (res_trial < res_best) THEN
444 22 : res_best = res_trial
445 22 : mu_best = mu
446 : END IF
447 : EXIT
448 : END IF
449 0 : step_try = 0.5_dp*step_try
450 : END DO
451 8 : IF (iback > NEWTON_MAX_BACKTRACK) EXIT
452 : END DO
453 8 : DEALLOCATE (gvec)
454 8 : mu = mu_best
455 :
456 : ! Final evaluation
457 8 : CALL SmearOcc(f, N_now, kTS, e, mu, sigma, maxocc, method, my_estate, my_festate)
458 8 : IF (ABS(N_now - N) >= N*1.0e-12_dp) THEN
459 0 : 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 774420 : mu_min = MINVAL(e)
465 26848 : iter = 0
466 22 : DO
467 26870 : iter = iter + 1
468 26870 : CALL SmearOcc(f, N_tmp, kTS, e, mu_min, sigma, maxocc, method, my_estate, my_festate)
469 26870 : IF (N_tmp <= N) EXIT
470 22 : IF (iter > 20) THEN
471 0 : CPABORT("SmearFixed: failed to bracket lower chemical potential")
472 : END IF
473 22 : mu_min = mu_min - sigma
474 : END DO
475 :
476 774420 : mu_max = MAXVAL(e)
477 26848 : iter = 0
478 0 : DO
479 26848 : iter = iter + 1
480 26848 : CALL SmearOcc(f, N_tmp, kTS, e, mu_max, sigma, maxocc, method, my_estate, my_festate)
481 26848 : IF (N_tmp >= N) EXIT
482 0 : IF (iter > 20) THEN
483 0 : CPABORT("SmearFixed: failed to bracket upper chemical potential")
484 : END IF
485 0 : mu_max = mu_max + sigma
486 : END DO
487 :
488 : iter = 0
489 1451588 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
490 1424740 : iter = iter + 1
491 1424740 : mu_now = (mu_max + mu_min)/2.0_dp
492 1424740 : CALL SmearOcc(f, N_now, kTS, e, mu_now, sigma, maxocc, method, my_estate, my_festate)
493 1424740 : IF (N_now <= N) THEN
494 700587 : mu_min = mu_now
495 : ELSE
496 724153 : mu_max = mu_now
497 : END IF
498 1451588 : IF (iter > BISECT_MAX_ITER) THEN
499 0 : CPWARN("SmearFixed: maximum bisection iterations reached")
500 0 : EXIT
501 : END IF
502 : END DO
503 :
504 26848 : mu = (mu_max + mu_min)/2.0_dp
505 53704 : CALL SmearOcc(f, N_now, kTS, e, mu, sigma, maxocc, method, my_estate, my_festate)
506 :
507 : END SELECT
508 :
509 26856 : 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 43538 : 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 43538 : nmo = SIZE(e, 1)
551 43538 : nkp = SIZE(e, 2)
552 :
553 43538 : CALL cite_smearing(method)
554 :
555 : ! Choose bisection method: Gaussian for MP/MV, actual method for FD/Gaussian
556 43814 : SELECT CASE (method)
557 : CASE (smear_mp, smear_mv)
558 276 : bisect_method = smear_gaussian
559 : CASE DEFAULT
560 43538 : bisect_method = method
561 : END SELECT
562 :
563 : ! Initial bracket
564 15524 : SELECT CASE (bisect_method)
565 : CASE (smear_fermi_dirac)
566 15524 : de = sigma*LOG((1.0_dp - epsocc)/epsocc)
567 : CASE DEFAULT
568 43538 : de = 10.0_dp*sigma
569 : END SELECT
570 43538 : de = MAX(de, 0.5_dp)
571 :
572 : ! Bisection with bisect_method
573 2380080 : mu_min = MINVAL(e) - de
574 2380080 : mu_max = MAXVAL(e) + de
575 43538 : iter = 0
576 299670 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
577 299652 : iter = iter + 1
578 299652 : mu = (mu_max + mu_min)/2.0_dp
579 299652 : CALL Smear2(f, N_now, kTS, e, mu, wk, sigma, maxocc, bisect_method)
580 :
581 299652 : IF (ABS(N_now - nel) < nel*epsocc) EXIT
582 :
583 256132 : IF (N_now <= nel) THEN
584 137440 : mu_min = mu
585 : ELSE
586 118692 : mu_max = mu
587 : END IF
588 :
589 299670 : IF (iter > BISECT_MAX_ITER) THEN
590 0 : CPWARN("Smearkp: maximum bisection iterations reached")
591 0 : EXIT
592 : END IF
593 : END DO
594 43538 : 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 276 : SELECT CASE (method)
600 : CASE (smear_mp, smear_mv)
601 276 : CALL Smear2(f, N_now, kTS, e, mu, wk, sigma, maxocc, method)
602 276 : res_best = ABS(N_now - nel)
603 276 : mu_best = mu
604 1132 : DO iter = 1, NEWTON_MAX_ITER
605 1132 : res_now = ABS(N_now - nel)
606 1132 : IF (res_now < nel*epsocc) EXIT
607 :
608 : ! Compute dN/dmu = sum_{ik} wk * g_i(k) inline
609 : dNdmu = 0.0_dp
610 8448 : DO ik = 1, nkp
611 76544 : DO is = 1, nmo
612 68096 : x = (e(is, ik) - mu)/sigma
613 7592 : SELECT CASE (method)
614 : CASE (smear_mp)
615 30384 : expx2 = EXP(-x*x)
616 30384 : dNdmu = dNdmu + maxocc*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2*wk(ik)
617 : CASE (smear_mv)
618 37712 : u = x + sqrthalf
619 37712 : expu2 = EXP(-u*u)
620 68096 : dNdmu = dNdmu + maxocc*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2*wk(ik)
621 : END SELECT
622 : END DO
623 : END DO
624 :
625 856 : IF (ABS(dNdmu) < EPSILON(dNdmu)) EXIT
626 856 : step = (nel - N_now)/dNdmu
627 856 : step = SIGN(MIN(ABS(step), MPMV_MAX_NEWTON_STEP*sigma), step)
628 856 : step_try = step
629 856 : DO iback = 1, NEWTON_MAX_BACKTRACK
630 856 : CALL Smear2(f, N_trial, kTS, e, mu + step_try, wk, sigma, maxocc, method)
631 856 : res_trial = ABS(N_trial - nel)
632 856 : IF (res_trial < res_now) THEN
633 856 : mu = mu + step_try
634 856 : N_now = N_trial
635 856 : IF (res_trial < res_best) THEN
636 856 : res_best = res_trial
637 856 : mu_best = mu
638 : END IF
639 : EXIT
640 : END IF
641 0 : step_try = 0.5_dp*step_try
642 : END DO
643 276 : IF (iback > NEWTON_MAX_BACKTRACK) EXIT
644 : END DO
645 43814 : mu = mu_best
646 : END SELECT
647 :
648 : ! Final evaluation with the actual method
649 43538 : CALL Smear2(f, N_now, kTS, e, mu, wk, sigma, maxocc, method)
650 276 : SELECT CASE (method)
651 : CASE (smear_mp, smear_mv)
652 43538 : IF (ABS(N_now - nel) >= nel*epsocc) THEN
653 0 : CPWARN("Smearkp: MP/MV smearing did not reach the requested electron count")
654 : END IF
655 : END SELECT
656 :
657 43538 : 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 532 : 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 532 : CPASSERT(SIZE(f, 3) == 2 .AND. SIZE(e, 3) == 2)
693 :
694 532 : nmo = SIZE(e, 1)
695 532 : nkp = SIZE(e, 2)
696 :
697 532 : CALL cite_smearing(method)
698 :
699 532 : SELECT CASE (method)
700 : CASE (smear_mp, smear_mv)
701 0 : bisect_method = smear_gaussian
702 : CASE DEFAULT
703 532 : bisect_method = method
704 : END SELECT
705 :
706 500 : SELECT CASE (bisect_method)
707 : CASE (smear_fermi_dirac)
708 500 : de = sigma*LOG((1.0_dp - epsocc)/epsocc)
709 : CASE DEFAULT
710 532 : de = 10.0_dp*sigma
711 : END SELECT
712 532 : de = MAX(de, 0.5_dp)
713 :
714 : ! Bisection with bisect_method
715 32472 : mu_min = MINVAL(e) - de
716 32472 : mu_max = MAXVAL(e) + de
717 532 : iter = 0
718 4178 : DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min)))
719 4178 : iter = iter + 1
720 4178 : mu = (mu_max + mu_min)/2.0_dp
721 4178 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, sigma, 1.0_dp, bisect_method)
722 4178 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, sigma, 1.0_dp, bisect_method)
723 4178 : N_now = na + nb
724 :
725 4178 : IF (ABS(N_now - nel) < nel*epsocc) EXIT
726 :
727 3646 : IF (N_now <= nel) THEN
728 1822 : mu_min = mu
729 : ELSE
730 1824 : mu_max = mu
731 : END IF
732 :
733 4178 : IF (iter > BISECT_MAX_ITER) THEN
734 0 : CPWARN("Smearkp2: maximum bisection iterations reached")
735 0 : EXIT
736 : END IF
737 : END DO
738 532 : 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 0 : SELECT CASE (method)
744 : CASE (smear_mp, smear_mv)
745 0 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
746 0 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
747 0 : N_now = na + nb
748 0 : res_best = ABS(N_now - nel)
749 0 : mu_best = mu
750 0 : DO iter = 1, NEWTON_MAX_ITER
751 0 : res_now = ABS(N_now - nel)
752 0 : IF (res_now < nel*epsocc) EXIT
753 :
754 : ! dN/dmu across both spin channels (maxocc=1 per spin)
755 : dNdmu = 0.0_dp
756 0 : DO ispin = 1, 2
757 0 : DO ik = 1, nkp
758 0 : DO is = 1, nmo
759 0 : x = (e(is, ik, ispin) - mu)/sigma
760 0 : SELECT CASE (method)
761 : CASE (smear_mp)
762 0 : expx2 = EXP(-x*x)
763 0 : dNdmu = dNdmu + (3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2*wk(ik)
764 : CASE (smear_mv)
765 0 : u = x + sqrthalf
766 0 : expu2 = EXP(-u*u)
767 0 : 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 0 : IF (ABS(dNdmu) < EPSILON(dNdmu)) EXIT
774 0 : step = (nel - N_now)/dNdmu
775 0 : step = SIGN(MIN(ABS(step), MPMV_MAX_NEWTON_STEP*sigma), step)
776 0 : step_try = step
777 0 : DO iback = 1, NEWTON_MAX_BACKTRACK
778 0 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu + step_try, wk, sigma, 1.0_dp, method)
779 0 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu + step_try, wk, sigma, 1.0_dp, method)
780 0 : N_trial = na + nb
781 0 : res_trial = ABS(N_trial - nel)
782 0 : IF (res_trial < res_now) THEN
783 0 : mu = mu + step_try
784 0 : N_now = N_trial
785 0 : IF (res_trial < res_best) THEN
786 0 : res_best = res_trial
787 0 : mu_best = mu
788 : END IF
789 : EXIT
790 : END IF
791 0 : step_try = 0.5_dp*step_try
792 : END DO
793 0 : IF (iback > NEWTON_MAX_BACKTRACK) EXIT
794 : END DO
795 532 : mu = mu_best
796 : END SELECT
797 :
798 : ! Final evaluation with the actual method
799 532 : CALL Smear2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, sigma, 1.0_dp, method)
800 532 : CALL Smear2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, sigma, 1.0_dp, method)
801 532 : N_now = na + nb
802 532 : kTS = kTSa + kTSb
803 0 : SELECT CASE (method)
804 : CASE (smear_mp, smear_mv)
805 532 : IF (ABS(N_now - nel) >= nel*epsocc) THEN
806 0 : CPWARN("Smearkp2: MP/MV smearing did not reach the requested electron count")
807 : END IF
808 : END SELECT
809 :
810 532 : 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 4590 : SUBROUTINE smearing_response_weight( &
840 4590 : 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 35620 : DO i = 1, Nstate
852 31030 : IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
853 418 : IF (i == estate) THEN
854 22 : occ_i = festate
855 : ELSE
856 396 : occ_i = maxocc
857 : END IF
858 : ELSE
859 30612 : occ_i = maxocc
860 : END IF
861 :
862 31030 : IF (occ_i < EPSILON(occ_i)) THEN
863 0 : gvec(i) = 0.0_dp
864 0 : CYCLE
865 : END IF
866 :
867 31030 : x = (e(i) - mu)/sigma
868 :
869 4590 : SELECT CASE (method)
870 : CASE (smear_fermi_dirac)
871 27864 : fi_norm = f(i)/occ_i
872 27864 : gvec(i) = occ_i*fi_norm*(1.0_dp - fi_norm)/sigma
873 :
874 : CASE (smear_gaussian)
875 916 : expx2 = EXP(-x*x)
876 916 : gvec(i) = occ_i/(sigma*rootpi)*expx2
877 :
878 : CASE (smear_mp)
879 1030 : expx2 = EXP(-x*x)
880 1030 : gvec(i) = occ_i*(3.0_dp - 2.0_dp*x*x)/(2.0_dp*sigma*rootpi)*expx2
881 :
882 : CASE (smear_mv)
883 1220 : u = x + sqrthalf
884 1220 : expu2 = EXP(-u*u)
885 31030 : gvec(i) = occ_i*(2.0_dp + sqrt2*x)/(sigma*rootpi)*expu2
886 :
887 : END SELECT
888 : END DO
889 :
890 4590 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE :: gvec(:)
934 :
935 0 : CALL timeset(routineN, handle)
936 :
937 : ! Step 1: find mu and f
938 0 : CALL SmearFixed(f, mu, kTS, e, N, sigma, maxocc, method, estate, festate)
939 :
940 : ! Step 2: build g vector
941 0 : Nstate = SIZE(e)
942 0 : ALLOCATE (gvec(Nstate))
943 0 : CALL smearing_response_weight(gvec, f, e, mu, sigma, maxocc, Nstate, method, estate, festate)
944 0 : Gsum = accurate_sum(gvec)
945 :
946 : ! Step 3: assemble dfde(i,j) = -delta_{ij}*g_i + g_i*g_j/G
947 0 : IF (ABS(Gsum) > EPSILON(Gsum)) THEN
948 0 : DO j = 1, Nstate
949 0 : DO i = 1, Nstate
950 0 : dfde(i, j) = gvec(i)*gvec(j)/Gsum
951 : END DO
952 0 : dfde(j, j) = dfde(j, j) - gvec(j)
953 : END DO
954 : ELSE
955 0 : dfde(:, :) = 0.0_dp
956 : END IF
957 :
958 0 : DEALLOCATE (gvec)
959 0 : CALL timestop(handle)
960 :
961 0 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE :: gvec(:)
1008 :
1009 0 : CALL timeset(routineN, handle)
1010 :
1011 : ! Step 1: find mu and f
1012 0 : CALL SmearFixed(f, mu, kTS, e, N_el, sigma, maxocc, method, estate, festate)
1013 :
1014 : ! Step 2: build g vector
1015 0 : Nstate = SIZE(e)
1016 0 : ALLOCATE (gvec(Nstate))
1017 0 : CALL smearing_response_weight(gvec, f, e, mu, sigma, maxocc, Nstate, method, estate, festate)
1018 0 : Gsum = accurate_sum(gvec)
1019 :
1020 : ! Step 3: RESULT_j = g_j * (g.v / G - v_j)
1021 0 : IF (ABS(Gsum) > EPSILON(Gsum)) THEN
1022 : gdotv = 0.0_dp
1023 0 : DO i = 1, Nstate
1024 0 : gdotv = gdotv + gvec(i)*v(i)
1025 : END DO
1026 0 : DO i = 1, Nstate
1027 0 : RESULT(i) = gvec(i)*(gdotv/Gsum - v(i))
1028 : END DO
1029 : ELSE
1030 0 : RESULT(:) = 0.0_dp
1031 : END IF
1032 :
1033 0 : DEALLOCATE (gvec)
1034 0 : CALL timestop(handle)
1035 :
1036 0 : END SUBROUTINE SmearFixedDerivMV
1037 :
1038 : END MODULE smearing_utils
|
