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 : MODULE qs_charge_mixing
10 :
11 : #if defined(__TBLITE)
12 : USE mctc_env, ONLY: error_type
13 : USE tblite_scc_mixer, ONLY: new_cp2k_tblite_mixer
14 : #endif
15 : USE input_constants, ONLY: tblite_mixer_damping_default, &
16 : tblite_mixer_iterations_default, &
17 : tblite_mixer_max_weight_default, &
18 : tblite_mixer_min_weight_default, &
19 : tblite_mixer_omega0_default, &
20 : tblite_mixer_weight_factor_default, &
21 : tblite_scc_mixer_auto, &
22 : tblite_scc_mixer_cp2k, &
23 : tblite_scc_mixer_none, &
24 : tblite_scc_mixer_tblite
25 : USE kinds, ONLY: dp
26 : USE mathlib, ONLY: get_pseudo_inverse_svd
27 : USE message_passing, ONLY: mp_para_env_type
28 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr, &
29 : gspace_mixing_nr, &
30 : mixing_storage_type, &
31 : modified_broyden_mixing_nr, &
32 : multisecant_mixing_nr, &
33 : new_pulay_mixing_nr, &
34 : pulay_mixing_nr
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_charge_mixing'
42 :
43 : PUBLIC :: charge_mixing, charge_mixing_scc_error, tblite_scc_error_on_cp2k_scale
44 :
45 : REAL(KIND=dp), PARAMETER, PUBLIC :: tblite_scc_pconv = 2.0E-5_dp
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Driver for TB SCC variable mixing, calls the requested method.
51 : !> \param mixing_method ...
52 : !> \param mixing_store ...
53 : !> \param charges ...
54 : !> \param para_env ...
55 : !> \param iter_count ...
56 : !> \param scc_mixer ...
57 : !> \param tblite_mixer_iterations ...
58 : !> \param tblite_mixer_damping ...
59 : !> \param tblite_mixer_memory ...
60 : !> \param tblite_mixer_omega0 ...
61 : !> \param tblite_mixer_min_weight ...
62 : !> \param tblite_mixer_max_weight ...
63 : !> \param tblite_mixer_weight_factor ...
64 : !> \par History
65 : !> \author JGH
66 : ! **************************************************************************************************
67 38814 : SUBROUTINE charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count, &
68 : scc_mixer, tblite_mixer_iterations, tblite_mixer_damping, &
69 : tblite_mixer_memory, tblite_mixer_omega0, tblite_mixer_min_weight, &
70 : tblite_mixer_max_weight, tblite_mixer_weight_factor)
71 : INTEGER, INTENT(IN) :: mixing_method
72 : TYPE(mixing_storage_type), POINTER :: mixing_store
73 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 : INTEGER, INTENT(IN) :: iter_count
76 : INTEGER, INTENT(IN), OPTIONAL :: scc_mixer, tblite_mixer_iterations
77 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: tblite_mixer_damping
78 : INTEGER, INTENT(IN), OPTIONAL :: tblite_mixer_memory
79 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: tblite_mixer_omega0, &
80 : tblite_mixer_min_weight, &
81 : tblite_mixer_max_weight, &
82 : tblite_mixer_weight_factor
83 :
84 : CHARACTER(len=*), PARAMETER :: routineN = 'charge_mixing'
85 :
86 : INTEGER :: effective_scc_mixer, handle, ia, ii, &
87 : imin, inow, nbuffer, ns, nvec
88 : REAL(dp) :: alpha
89 : #if defined(__TBLITE)
90 : INTEGER :: mixer_iterations, mixer_memory
91 : REAL(dp) :: mixer_damping, mixer_max_weight, &
92 : mixer_min_weight, mixer_omega0, &
93 : mixer_weight_factor
94 : #endif
95 :
96 38814 : CALL timeset(routineN, handle)
97 :
98 38814 : effective_scc_mixer = tblite_scc_mixer_cp2k
99 38814 : IF (PRESENT(scc_mixer)) effective_scc_mixer = scc_mixer
100 38814 : IF (ASSOCIATED(mixing_store)) mixing_store%tb_scc_mixer_error = 0.0_dp
101 :
102 32 : SELECT CASE (effective_scc_mixer)
103 : CASE (tblite_scc_mixer_auto, tblite_scc_mixer_cp2k)
104 : ! Use the regular CP2K SCC-variable mixing path below.
105 : CASE (tblite_scc_mixer_tblite)
106 32 : CPASSERT(ASSOCIATED(mixing_store))
107 : #if defined(__TBLITE)
108 32 : mixer_damping = tblite_mixer_damping_default
109 32 : IF (PRESENT(tblite_mixer_damping)) mixer_damping = tblite_mixer_damping
110 32 : IF (mixer_damping <= 0.0_dp) CPABORT("tblite SCC mixer DAMPING must be positive")
111 32 : mixer_omega0 = tblite_mixer_omega0_default
112 32 : IF (PRESENT(tblite_mixer_omega0)) mixer_omega0 = tblite_mixer_omega0
113 32 : IF (mixer_omega0 <= 0.0_dp) CPABORT("tblite SCC mixer OMEGA0 must be positive")
114 32 : mixer_min_weight = tblite_mixer_min_weight_default
115 32 : IF (PRESENT(tblite_mixer_min_weight)) mixer_min_weight = tblite_mixer_min_weight
116 32 : IF (mixer_min_weight <= 0.0_dp) CPABORT("tblite SCC mixer MIN_WEIGHT must be positive")
117 32 : mixer_max_weight = tblite_mixer_max_weight_default
118 32 : IF (PRESENT(tblite_mixer_max_weight)) mixer_max_weight = tblite_mixer_max_weight
119 32 : IF (mixer_max_weight <= 0.0_dp) CPABORT("tblite SCC mixer MAX_WEIGHT must be positive")
120 32 : IF (mixer_max_weight < mixer_min_weight) THEN
121 0 : CPABORT("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT")
122 : END IF
123 32 : mixer_weight_factor = tblite_mixer_weight_factor_default
124 32 : IF (PRESENT(tblite_mixer_weight_factor)) mixer_weight_factor = tblite_mixer_weight_factor
125 32 : IF (mixer_weight_factor <= 0.0_dp) CPABORT("tblite SCC mixer WEIGHT_FACTOR must be positive")
126 32 : mixer_iterations = tblite_mixer_iterations_default
127 32 : IF (PRESENT(tblite_mixer_iterations)) mixer_iterations = tblite_mixer_iterations
128 32 : IF (mixer_iterations < 1) CPABORT("tblite SCC mixer ITERATIONS must be positive")
129 32 : IF (iter_count > mixer_iterations) CPABORT("tblite SCC mixer exceeded ITERATIONS")
130 32 : mixer_memory = MAX(1, mixing_store%nbuffer)
131 32 : IF (PRESENT(tblite_mixer_memory)) mixer_memory = tblite_mixer_memory
132 32 : IF (mixer_memory < 1) CPABORT("tblite SCC mixer MEMORY must be positive")
133 : CALL tblite_charge_mixing(mixing_store, charges, para_env, iter_count, &
134 : mixer_damping, mixer_memory, mixer_omega0, mixer_min_weight, &
135 32 : mixer_max_weight, mixer_weight_factor)
136 32 : CALL timestop(handle)
137 32 : RETURN
138 : #else
139 : MARK_USED(tblite_mixer_damping)
140 : MARK_USED(tblite_mixer_iterations)
141 : MARK_USED(tblite_mixer_max_weight)
142 : MARK_USED(tblite_mixer_memory)
143 : MARK_USED(tblite_mixer_min_weight)
144 : MARK_USED(tblite_mixer_omega0)
145 : MARK_USED(tblite_mixer_weight_factor)
146 : IF (iter_count == 1) THEN
147 : CALL cp_warn(__LOCATION__, &
148 : "SCC_MIXER TBLITE requested but CP2K was built without tblite; "// &
149 : "falling back to the CP2K SCC mixer.")
150 : END IF
151 : #endif
152 : CASE (tblite_scc_mixer_none)
153 0 : IF (ASSOCIATED(mixing_store)) mixing_store%iter_method = "NoMix"
154 0 : CALL timestop(handle)
155 0 : RETURN
156 : CASE DEFAULT
157 38814 : CPABORT("Unknown SCC mixer for TB charge mixing")
158 : END SELECT
159 :
160 38782 : IF (mixing_method >= gspace_mixing_nr) THEN
161 1800 : CPASSERT(ASSOCIATED(mixing_store))
162 1800 : mixing_store%ncall = mixing_store%ncall + 1
163 1800 : ns = SIZE(charges, 2)
164 1800 : IF (ns > mixing_store%max_shell) THEN
165 0 : CPABORT("Mixing storage too small for TB SCC variables")
166 : END IF
167 1800 : alpha = mixing_store%alpha
168 1800 : nbuffer = mixing_store%nbuffer
169 1800 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
170 1800 : imin = inow - 1
171 1800 : IF (imin == 0) imin = nbuffer
172 1800 : IF (mixing_store%ncall > nbuffer) THEN
173 734 : nvec = nbuffer
174 : ELSE
175 1066 : nvec = mixing_store%ncall - 1
176 : END IF
177 1800 : IF (mixing_store%ncall > 1) THEN
178 : ! store in/out charge difference
179 5434 : DO ia = 1, mixing_store%nat_local
180 3788 : ii = mixing_store%atlist(ia)
181 17902 : mixing_store%dacharge(ia, 1:ns, imin) = mixing_store%acharge(ia, 1:ns, imin) - charges(ii, 1:ns)
182 : END DO
183 : END IF
184 1800 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
185 : ! skip mixing
186 166 : mixing_store%iter_method = "NoMix"
187 1634 : ELSE IF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
188 136 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
189 136 : mixing_store%iter_method = "Mixing"
190 : ELSE IF (mixing_method == gspace_mixing_nr) THEN
191 0 : CPABORT("Kerker method not available for Charge Mixing")
192 : ELSE IF (mixing_method == pulay_mixing_nr) THEN
193 0 : CPABORT("Pulay method not available for Charge Mixing")
194 : ELSE IF (mixing_method == broyden_mixing_nr) THEN
195 1410 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env, modified=.FALSE.)
196 1410 : mixing_store%iter_method = "Broy."
197 : ELSE IF (mixing_method == modified_broyden_mixing_nr) THEN
198 88 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env, modified=.TRUE.)
199 88 : mixing_store%iter_method = "MBroy"
200 : ELSE IF (mixing_method == multisecant_mixing_nr) THEN
201 0 : CPABORT("Multisecant_mixing method not available for Charge Mixing")
202 : ELSE IF (mixing_method == new_pulay_mixing_nr) THEN
203 0 : CPABORT("New Pulay method not available for Charge Mixing")
204 : END IF
205 :
206 : ! store new 'input' charges
207 5998 : DO ia = 1, mixing_store%nat_local
208 4198 : ii = mixing_store%atlist(ia)
209 19898 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
210 : END DO
211 :
212 : END IF
213 :
214 38782 : CALL timestop(handle)
215 :
216 : END SUBROUTINE charge_mixing
217 :
218 : ! **************************************************************************************************
219 : !> \brief Map a raw tblite SCC residual to CP2K's EPS_SCF reporting scale.
220 : !> \param raw_error raw tblite SCC residual
221 : !> \param eps_scf CP2K SCF convergence threshold
222 : !> \param pconv tblite SCC convergence reference
223 : !> \return residual on the CP2K convergence scale
224 : ! **************************************************************************************************
225 22384 : PURE FUNCTION tblite_scc_error_on_cp2k_scale(raw_error, eps_scf, pconv) RESULT(scaled_error)
226 : REAL(KIND=dp), INTENT(IN) :: raw_error, eps_scf, pconv
227 : REAL(KIND=dp) :: scaled_error
228 :
229 22384 : IF (eps_scf > 0.0_dp .AND. pconv > 0.0_dp) THEN
230 22384 : scaled_error = eps_scf*raw_error/pconv
231 : ELSE
232 0 : scaled_error = raw_error
233 : END IF
234 :
235 22384 : END FUNCTION tblite_scc_error_on_cp2k_scale
236 :
237 : ! **************************************************************************************************
238 : !> \brief Return the CP2K-side tblite SCC-mixer residual on the CP2K EPS_SCF scale.
239 : !> \param mixing_store ...
240 : !> \param eps_scf ...
241 : !> \return ...
242 : ! **************************************************************************************************
243 56572 : FUNCTION charge_mixing_scc_error(mixing_store, eps_scf) RESULT(mixer_error)
244 : TYPE(mixing_storage_type), POINTER :: mixing_store
245 : REAL(KIND=dp), INTENT(IN) :: eps_scf
246 : REAL(KIND=dp) :: mixer_error
247 :
248 56572 : mixer_error = 0.0_dp
249 56572 : IF (.NOT. ASSOCIATED(mixing_store)) RETURN
250 56572 : IF (mixing_store%tb_scc_mixer_step <= 1) RETURN
251 :
252 : mixer_error = tblite_scc_error_on_cp2k_scale(mixing_store%tb_scc_mixer_error, &
253 26 : eps_scf, tblite_scc_pconv)
254 :
255 26 : END FUNCTION charge_mixing_scc_error
256 :
257 : ! **************************************************************************************************
258 : !> \brief TBLite modified-Broyden mixing for a complete TB SCC-variable vector.
259 : !> \param mixing_store ...
260 : !> \param charges ...
261 : !> \param para_env ...
262 : !> \param iter_count ...
263 : !> \param damping ...
264 : !> \param memory ...
265 : !> \param omega0 ...
266 : !> \param min_weight ...
267 : !> \param max_weight ...
268 : !> \param weight_factor ...
269 : ! **************************************************************************************************
270 32 : SUBROUTINE tblite_charge_mixing(mixing_store, charges, para_env, iter_count, damping, memory, omega0, &
271 : min_weight, max_weight, weight_factor)
272 : TYPE(mixing_storage_type), POINTER :: mixing_store
273 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
274 : TYPE(mp_para_env_type), POINTER :: para_env
275 : INTEGER, INTENT(IN) :: iter_count, memory
276 : REAL(KIND=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, &
277 : weight_factor
278 :
279 : #if defined(__TBLITE)
280 32 : TYPE(error_type), ALLOCATABLE :: error
281 : #endif
282 : INTEGER :: natom, ndim, ns
283 : LOGICAL :: on_source, reset_mixer
284 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: qvec
285 :
286 32 : natom = SIZE(charges, 1)
287 32 : ns = SIZE(charges, 2)
288 32 : ndim = natom*ns
289 96 : ALLOCATE (qvec(ndim))
290 64 : qvec(:) = RESHAPE(charges, [ndim])
291 32 : on_source = para_env%mepos == para_env%source
292 : reset_mixer = (iter_count == 1) .OR. (mixing_store%tb_scc_mixer_step == 0) .OR. &
293 : (mixing_store%tb_scc_mixer_natom /= natom) .OR. &
294 : (mixing_store%tb_scc_mixer_ns /= ns) .OR. &
295 32 : (mixing_store%tb_scc_mixer_memory /= memory)
296 32 : mixing_store%tb_scc_mixer_error = 0.0_dp
297 :
298 : #if defined(__TBLITE)
299 32 : IF (reset_mixer) THEN
300 4 : IF (ALLOCATED(mixing_store%tb_scc_mixer)) DEALLOCATE (mixing_store%tb_scc_mixer)
301 4 : IF (on_source) THEN
302 : CALL new_cp2k_tblite_mixer(mixing_store%tb_scc_mixer, memory, ndim, damping, omega0, &
303 2 : min_weight, max_weight, weight_factor)
304 2 : CALL mixing_store%tb_scc_mixer%set(qvec)
305 : END IF
306 4 : mixing_store%tb_scc_mixer_natom = natom
307 4 : mixing_store%tb_scc_mixer_ns = ns
308 4 : mixing_store%tb_scc_mixer_memory = memory
309 4 : mixing_store%tb_scc_mixer_step = 1
310 4 : mixing_store%iter_method = "NoMix"
311 4 : CALL para_env%bcast(qvec)
312 12 : charges = RESHAPE(qvec, SHAPE(charges))
313 4 : RETURN
314 : END IF
315 :
316 28 : IF (on_source) THEN
317 14 : CPASSERT(ALLOCATED(mixing_store%tb_scc_mixer))
318 14 : CALL mixing_store%tb_scc_mixer%diff(qvec)
319 14 : mixing_store%tb_scc_mixer_error = REAL(mixing_store%tb_scc_mixer%get_error(), KIND=dp)
320 14 : CALL mixing_store%tb_scc_mixer%next(error)
321 14 : IF (ALLOCATED(error)) CPABORT("tblite SCC mixer failed")
322 14 : CALL mixing_store%tb_scc_mixer%get(qvec)
323 : END IF
324 28 : CALL para_env%bcast(qvec)
325 28 : CALL para_env%bcast(mixing_store%tb_scc_mixer_error)
326 84 : charges = RESHAPE(qvec, SHAPE(charges))
327 28 : mixing_store%tb_scc_mixer_step = mixing_store%tb_scc_mixer_step + 1
328 28 : mixing_store%iter_method = "TBLITE"
329 : #else
330 : MARK_USED(mixing_store)
331 : MARK_USED(charges)
332 : MARK_USED(para_env)
333 : MARK_USED(iter_count)
334 : MARK_USED(damping)
335 : MARK_USED(memory)
336 : MARK_USED(omega0)
337 : MARK_USED(min_weight)
338 : MARK_USED(max_weight)
339 : MARK_USED(weight_factor)
340 : CPABORT("SCC_MIXER TBLITE requires CP2K to be built with tblite")
341 : #endif
342 :
343 32 : END SUBROUTINE tblite_charge_mixing
344 :
345 : ! **************************************************************************************************
346 : !> \brief Simple charge mixing
347 : !> \param mixing_store ...
348 : !> \param charges ...
349 : !> \param alpha ...
350 : !> \param imin ...
351 : !> \param ns ...
352 : !> \param para_env ...
353 : !> \author JGH
354 : ! **************************************************************************************************
355 136 : SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
356 : TYPE(mixing_storage_type), POINTER :: mixing_store
357 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
358 : REAL(KIND=dp), INTENT(IN) :: alpha
359 : INTEGER, INTENT(IN) :: imin, ns
360 : TYPE(mp_para_env_type), POINTER :: para_env
361 :
362 : INTEGER :: ia, ii
363 :
364 2870 : charges = 0.0_dp
365 :
366 492 : DO ia = 1, mixing_store%nat_local
367 356 : ii = mixing_store%atlist(ia)
368 1608 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
369 : END DO
370 :
371 5604 : CALL para_env%sum(charges)
372 :
373 136 : END SUBROUTINE mix_charges_only
374 :
375 : ! **************************************************************************************************
376 : !> \brief Broyden charge mixing
377 : !> \param mixing_store ...
378 : !> \param charges ...
379 : !> \param inow ...
380 : !> \param nvec ...
381 : !> \param ns ...
382 : !> \param para_env ...
383 : !> \param modified use dynamic residual weights of modified Broyden
384 : !> \author JGH
385 : ! **************************************************************************************************
386 1498 : SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env, modified)
387 : TYPE(mixing_storage_type), POINTER :: mixing_store
388 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
389 : INTEGER, INTENT(IN) :: inow, nvec, ns
390 : TYPE(mp_para_env_type), POINTER :: para_env
391 : LOGICAL, INTENT(IN) :: modified
392 :
393 : INTEGER :: i, ia, ii, imin, j, nbuffer, nv
394 : REAL(KIND=dp) :: alpha, broy_w0, res_norm, rskip, wdf, &
395 : wprod
396 1498 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
397 1498 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
398 1498 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
399 :
400 0 : CPASSERT(nvec > 1)
401 :
402 1498 : nbuffer = mixing_store%nbuffer
403 1498 : alpha = mixing_store%alpha
404 1498 : imin = inow - 1
405 1498 : IF (imin == 0) imin = nvec
406 1498 : nv = nvec - 1
407 :
408 : ! charge vectors
409 1498 : q_now => mixing_store%acharge(:, :, inow)
410 1498 : q_last => mixing_store%acharge(:, :, imin)
411 1498 : dq_now => mixing_store%dacharge(:, :, inow)
412 1498 : dq_last => mixing_store%dacharge(:, :, imin)
413 :
414 1498 : IF (nvec == nbuffer) THEN
415 : ! reshuffel Broyden storage n->n-1
416 4938 : DO i = 1, nv - 1
417 4204 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
418 191324 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
419 192058 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
420 : END DO
421 4938 : DO i = 1, nv - 1
422 29810 : DO j = 1, nv - 1
423 29076 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
424 : END DO
425 : END DO
426 : END IF
427 :
428 1498 : broy_w0 = mixing_store%broy_w0
429 1498 : IF (modified) THEN
430 1720 : res_norm = SUM(dq_now(:, 1:ns)**2)
431 88 : CALL para_env%sum(res_norm)
432 88 : res_norm = SQRT(res_norm)
433 88 : IF (res_norm > mixing_store%wc/mixing_store%wmax) THEN
434 66 : mixing_store%wbroy(nv) = mixing_store%wc/res_norm
435 : ELSE
436 22 : mixing_store%wbroy(nv) = mixing_store%wmax
437 : END IF
438 88 : mixing_store%wbroy(nv) = MAX(1.0_dp, mixing_store%wbroy(nv))
439 : ELSE
440 1410 : mixing_store%wbroy(nv) = 1.0_dp
441 : END IF
442 :
443 : ! dfbroy
444 81954 : mixing_store%dfbroy(:, :, nv) = 0.0_dp
445 37198 : mixing_store%dfbroy(:, 1:ns, nv) = dq_now(:, 1:ns) - dq_last(:, 1:ns)
446 19348 : wdf = SUM(mixing_store%dfbroy(:, 1:ns, nv)**2)
447 1498 : CALL para_env%sum(wdf)
448 1498 : IF (wdf > TINY(1.0_dp) .AND. wdf < HUGE(1.0_dp)) THEN
449 1496 : wdf = 1.0_dp/SQRT(wdf)
450 19342 : mixing_store%dfbroy(:, 1:ns, nv) = wdf*mixing_store%dfbroy(:, 1:ns, nv)
451 : ELSE
452 : ! Identical consecutive residuals do not define a Broyden direction.
453 : ! Keep a zero history vector so it does not enter the Broyden update.
454 2 : wdf = 0.0_dp
455 6 : mixing_store%dfbroy(:, 1:ns, nv) = 0.0_dp
456 : END IF
457 :
458 : ! abroy matrix
459 9208 : DO i = 1, nv
460 98758 : wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*mixing_store%dfbroy(:, 1:ns, nv))
461 7710 : CALL para_env%sum(wprod)
462 7710 : mixing_store%abroy(i, nv) = wprod
463 9208 : mixing_store%abroy(nv, i) = wprod
464 : END DO
465 :
466 : ! broyden matrices
467 13482 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
468 9208 : DO i = 1, nv
469 98758 : wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*dq_now(:, 1:ns))
470 7710 : CALL para_env%sum(wprod)
471 9208 : cvec(i) = mixing_store%wbroy(i)*wprod
472 : END DO
473 :
474 9208 : DO i = 1, nv
475 55356 : DO j = 1, nv
476 55356 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
477 : END DO
478 9208 : IF (modified) THEN
479 532 : beta(i, i) = beta(i, i) + broy_w0*broy_w0
480 : ELSE
481 7178 : beta(i, i) = beta(i, i) + broy_w0
482 : END IF
483 : END DO
484 :
485 1498 : rskip = 1.e-12_dp
486 1498 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
487 64564 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
488 :
489 : ! build ubroy
490 81954 : mixing_store%ubroy(:, :, nv) = 0.0_dp
491 : mixing_store%ubroy(:, 1:ns, nv) = alpha*mixing_store%dfbroy(:, 1:ns, nv) + &
492 37198 : wdf*(q_now(:, 1:ns) - q_last(:, 1:ns))
493 :
494 30552 : charges = 0.0_dp
495 4910 : DO ia = 1, mixing_store%nat_local
496 3412 : ii = mixing_store%atlist(ia)
497 16116 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
498 : END DO
499 9208 : DO i = 1, nv
500 25994 : DO ia = 1, mixing_store%nat_local
501 16786 : ii = mixing_store%atlist(ia)
502 80306 : charges(ii, 1:ns) = charges(ii, 1:ns) - mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
503 : END DO
504 : END DO
505 59606 : CALL para_env%sum(charges)
506 :
507 1498 : DEALLOCATE (amat, beta, cvec, gammab)
508 :
509 1498 : END SUBROUTINE broyden_mixing
510 :
511 : END MODULE qs_charge_mixing
|
