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