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
44 :
45 : REAL(KIND=dp), PARAMETER, PRIVATE :: 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 34534 : 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 34534 : CALL timeset(routineN, handle)
97 :
98 34534 : effective_scc_mixer = tblite_scc_mixer_cp2k
99 34534 : IF (PRESENT(scc_mixer)) effective_scc_mixer = scc_mixer
100 34534 : 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) &
121 0 : CPABORT("tblite SCC mixer MAX_WEIGHT must not be smaller than MIN_WEIGHT")
122 32 : mixer_weight_factor = tblite_mixer_weight_factor_default
123 32 : IF (PRESENT(tblite_mixer_weight_factor)) mixer_weight_factor = tblite_mixer_weight_factor
124 32 : IF (mixer_weight_factor <= 0.0_dp) CPABORT("tblite SCC mixer WEIGHT_FACTOR must be positive")
125 32 : mixer_iterations = tblite_mixer_iterations_default
126 32 : IF (PRESENT(tblite_mixer_iterations)) mixer_iterations = tblite_mixer_iterations
127 32 : IF (mixer_iterations < 1) CPABORT("tblite SCC mixer ITERATIONS must be positive")
128 32 : IF (iter_count > mixer_iterations) CPABORT("tblite SCC mixer exceeded ITERATIONS")
129 32 : mixer_memory = MAX(1, mixing_store%nbuffer)
130 32 : IF (PRESENT(tblite_mixer_memory)) mixer_memory = tblite_mixer_memory
131 32 : IF (mixer_memory < 1) CPABORT("tblite SCC mixer MEMORY must be positive")
132 : CALL tblite_charge_mixing(mixing_store, charges, para_env, iter_count, &
133 : mixer_damping, mixer_memory, mixer_omega0, mixer_min_weight, &
134 32 : mixer_max_weight, mixer_weight_factor)
135 32 : CALL timestop(handle)
136 32 : RETURN
137 : #else
138 : MARK_USED(tblite_mixer_damping)
139 : MARK_USED(tblite_mixer_iterations)
140 : MARK_USED(tblite_mixer_max_weight)
141 : MARK_USED(tblite_mixer_memory)
142 : MARK_USED(tblite_mixer_min_weight)
143 : MARK_USED(tblite_mixer_omega0)
144 : MARK_USED(tblite_mixer_weight_factor)
145 : IF (iter_count == 1) THEN
146 : CALL cp_warn(__LOCATION__, &
147 : "SCC_MIXER TBLITE requested but CP2K was built without tblite; "// &
148 : "falling back to the CP2K SCC mixer.")
149 : END IF
150 : #endif
151 : CASE (tblite_scc_mixer_none)
152 0 : IF (ASSOCIATED(mixing_store)) mixing_store%iter_method = "NoMix"
153 0 : CALL timestop(handle)
154 0 : RETURN
155 : CASE DEFAULT
156 34534 : CPABORT("Unknown SCC mixer for TB charge mixing")
157 : END SELECT
158 :
159 34502 : IF (mixing_method >= gspace_mixing_nr) THEN
160 1086 : CPASSERT(ASSOCIATED(mixing_store))
161 1086 : mixing_store%ncall = mixing_store%ncall + 1
162 1086 : ns = SIZE(charges, 2)
163 1086 : IF (ns > mixing_store%max_shell) THEN
164 0 : CPABORT("Mixing storage too small for TB SCC variables")
165 : END IF
166 1086 : alpha = mixing_store%alpha
167 1086 : nbuffer = mixing_store%nbuffer
168 1086 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
169 1086 : imin = inow - 1
170 1086 : IF (imin == 0) imin = nbuffer
171 1086 : IF (mixing_store%ncall > nbuffer) THEN
172 300 : nvec = nbuffer
173 : ELSE
174 786 : nvec = mixing_store%ncall - 1
175 : END IF
176 1086 : IF (mixing_store%ncall > 1) THEN
177 : ! store in/out charge difference
178 4064 : DO ia = 1, mixing_store%nat_local
179 3094 : ii = mixing_store%atlist(ia)
180 13050 : mixing_store%dacharge(ia, 1:ns, imin) = mixing_store%acharge(ia, 1:ns, imin) - charges(ii, 1:ns)
181 : END DO
182 : END IF
183 1086 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
184 : ! skip mixing
185 126 : mixing_store%iter_method = "NoMix"
186 960 : ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
187 102 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
188 102 : mixing_store%iter_method = "Mixing"
189 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
190 0 : CPABORT("Kerker method not available for Charge Mixing")
191 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
192 0 : CPABORT("Pulay method not available for Charge Mixing")
193 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
194 858 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
195 858 : mixing_store%iter_method = "Broy."
196 : ELSEIF (mixing_method == modified_broyden_mixing_nr) THEN
197 0 : CPABORT("Modified Broyden mixing is only available for DFT density mixing")
198 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
199 0 : CPABORT("Multisecant_mixing method not available for Charge Mixing")
200 : ELSEIF (mixing_method == new_pulay_mixing_nr) THEN
201 0 : CPABORT("New Pulay method not available for Charge Mixing")
202 : END IF
203 :
204 : ! store new 'input' charges
205 4524 : DO ia = 1, mixing_store%nat_local
206 3438 : ii = mixing_store%atlist(ia)
207 14624 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
208 : END DO
209 :
210 : END IF
211 :
212 34502 : CALL timestop(handle)
213 :
214 : END SUBROUTINE charge_mixing
215 :
216 : ! **************************************************************************************************
217 : !> \brief Return the CP2K-side tblite SCC-mixer residual on the CP2K EPS_SCF scale.
218 : !> \param mixing_store ...
219 : !> \param eps_scf ...
220 : !> \return ...
221 : ! **************************************************************************************************
222 52768 : FUNCTION charge_mixing_scc_error(mixing_store, eps_scf) RESULT(mixer_error)
223 : TYPE(mixing_storage_type), POINTER :: mixing_store
224 : REAL(KIND=dp), INTENT(IN) :: eps_scf
225 : REAL(KIND=dp) :: mixer_error
226 :
227 52768 : mixer_error = 0.0_dp
228 52768 : IF (.NOT. ASSOCIATED(mixing_store)) RETURN
229 52768 : IF (mixing_store%tb_scc_mixer_step <= 1) RETURN
230 :
231 26 : IF (eps_scf > 0.0_dp) THEN
232 26 : mixer_error = eps_scf*mixing_store%tb_scc_mixer_error/tblite_scc_pconv
233 : ELSE
234 0 : mixer_error = mixing_store%tb_scc_mixer_error
235 : END IF
236 :
237 : END FUNCTION charge_mixing_scc_error
238 :
239 : ! **************************************************************************************************
240 : !> \brief TBLite modified-Broyden mixing for a complete TB SCC-variable vector.
241 : !> \param mixing_store ...
242 : !> \param charges ...
243 : !> \param para_env ...
244 : !> \param iter_count ...
245 : !> \param damping ...
246 : !> \param memory ...
247 : !> \param omega0 ...
248 : !> \param min_weight ...
249 : !> \param max_weight ...
250 : !> \param weight_factor ...
251 : ! **************************************************************************************************
252 32 : SUBROUTINE tblite_charge_mixing(mixing_store, charges, para_env, iter_count, damping, memory, omega0, &
253 : min_weight, max_weight, weight_factor)
254 : TYPE(mixing_storage_type), POINTER :: mixing_store
255 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
256 : TYPE(mp_para_env_type), POINTER :: para_env
257 : INTEGER, INTENT(IN) :: iter_count, memory
258 : REAL(KIND=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, &
259 : weight_factor
260 :
261 : #if defined(__TBLITE)
262 32 : TYPE(error_type), ALLOCATABLE :: error
263 : #endif
264 : INTEGER :: natom, ndim, ns
265 : LOGICAL :: on_source, reset_mixer
266 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: qvec
267 :
268 32 : natom = SIZE(charges, 1)
269 32 : ns = SIZE(charges, 2)
270 32 : ndim = natom*ns
271 96 : ALLOCATE (qvec(ndim))
272 64 : qvec(:) = RESHAPE(charges, [ndim])
273 32 : on_source = para_env%mepos == para_env%source
274 : reset_mixer = (iter_count == 1) .OR. (mixing_store%tb_scc_mixer_step == 0) .OR. &
275 : (mixing_store%tb_scc_mixer_natom /= natom) .OR. &
276 : (mixing_store%tb_scc_mixer_ns /= ns) .OR. &
277 32 : (mixing_store%tb_scc_mixer_memory /= memory)
278 32 : mixing_store%tb_scc_mixer_error = 0.0_dp
279 :
280 : #if defined(__TBLITE)
281 32 : IF (reset_mixer) THEN
282 4 : IF (ALLOCATED(mixing_store%tb_scc_mixer)) DEALLOCATE (mixing_store%tb_scc_mixer)
283 4 : IF (on_source) THEN
284 : CALL new_cp2k_tblite_mixer(mixing_store%tb_scc_mixer, memory, ndim, damping, omega0, &
285 2 : min_weight, max_weight, weight_factor)
286 2 : CALL mixing_store%tb_scc_mixer%set(qvec)
287 : END IF
288 4 : mixing_store%tb_scc_mixer_natom = natom
289 4 : mixing_store%tb_scc_mixer_ns = ns
290 4 : mixing_store%tb_scc_mixer_memory = memory
291 4 : mixing_store%tb_scc_mixer_step = 1
292 4 : mixing_store%iter_method = "NoMix"
293 4 : CALL para_env%bcast(qvec)
294 12 : charges = RESHAPE(qvec, SHAPE(charges))
295 4 : RETURN
296 : END IF
297 :
298 28 : IF (on_source) THEN
299 14 : CPASSERT(ALLOCATED(mixing_store%tb_scc_mixer))
300 14 : CALL mixing_store%tb_scc_mixer%diff(qvec)
301 14 : mixing_store%tb_scc_mixer_error = REAL(mixing_store%tb_scc_mixer%get_error(), KIND=dp)
302 14 : CALL mixing_store%tb_scc_mixer%next(error)
303 14 : IF (ALLOCATED(error)) CPABORT("tblite SCC mixer failed")
304 14 : CALL mixing_store%tb_scc_mixer%get(qvec)
305 : END IF
306 28 : CALL para_env%bcast(qvec)
307 28 : CALL para_env%bcast(mixing_store%tb_scc_mixer_error)
308 84 : charges = RESHAPE(qvec, SHAPE(charges))
309 28 : mixing_store%tb_scc_mixer_step = mixing_store%tb_scc_mixer_step + 1
310 28 : mixing_store%iter_method = "TBLITE"
311 : #else
312 : MARK_USED(mixing_store)
313 : MARK_USED(charges)
314 : MARK_USED(para_env)
315 : MARK_USED(iter_count)
316 : MARK_USED(damping)
317 : MARK_USED(memory)
318 : MARK_USED(omega0)
319 : MARK_USED(min_weight)
320 : MARK_USED(max_weight)
321 : MARK_USED(weight_factor)
322 : CPABORT("SCC_MIXER TBLITE requires CP2K to be built with tblite")
323 : #endif
324 :
325 32 : END SUBROUTINE tblite_charge_mixing
326 :
327 : ! **************************************************************************************************
328 : !> \brief Simple charge mixing
329 : !> \param mixing_store ...
330 : !> \param charges ...
331 : !> \param alpha ...
332 : !> \param imin ...
333 : !> \param ns ...
334 : !> \param para_env ...
335 : !> \author JGH
336 : ! **************************************************************************************************
337 102 : SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
338 : TYPE(mixing_storage_type), POINTER :: mixing_store
339 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
340 : REAL(KIND=dp), INTENT(IN) :: alpha
341 : INTEGER, INTENT(IN) :: imin, ns
342 : TYPE(mp_para_env_type), POINTER :: para_env
343 :
344 : INTEGER :: ia, ii
345 :
346 2372 : charges = 0.0_dp
347 :
348 422 : DO ia = 1, mixing_store%nat_local
349 320 : ii = mixing_store%atlist(ia)
350 1382 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
351 : END DO
352 :
353 4642 : CALL para_env%sum(charges)
354 :
355 102 : END SUBROUTINE mix_charges_only
356 :
357 : ! **************************************************************************************************
358 : !> \brief Broyden charge mixing
359 : !> \param mixing_store ...
360 : !> \param charges ...
361 : !> \param inow ...
362 : !> \param nvec ...
363 : !> \param ns ...
364 : !> \param para_env ...
365 : !> \author JGH
366 : ! **************************************************************************************************
367 858 : SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env)
368 : TYPE(mixing_storage_type), POINTER :: mixing_store
369 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
370 : INTEGER, INTENT(IN) :: inow, nvec, ns
371 : TYPE(mp_para_env_type), POINTER :: para_env
372 :
373 : INTEGER :: i, ia, ii, imin, j, nbuffer, nv
374 : REAL(KIND=dp) :: alpha, broy_w0, rskip, wdf, wprod
375 858 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
376 858 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
377 858 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
378 :
379 0 : CPASSERT(nvec > 1)
380 :
381 858 : nbuffer = mixing_store%nbuffer
382 858 : alpha = mixing_store%alpha
383 858 : imin = inow - 1
384 858 : IF (imin == 0) imin = nvec
385 858 : nv = nvec - 1
386 :
387 : ! charge vectors
388 858 : q_now => mixing_store%acharge(:, :, inow)
389 858 : q_last => mixing_store%acharge(:, :, imin)
390 858 : dq_now => mixing_store%dacharge(:, :, inow)
391 858 : dq_last => mixing_store%dacharge(:, :, imin)
392 :
393 858 : IF (nvec == nbuffer) THEN
394 : ! reshuffel Broyden storage n->n-1
395 1900 : DO i = 1, nv - 1
396 1600 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
397 44168 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
398 44468 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
399 : END DO
400 1900 : DO i = 1, nv - 1
401 11148 : DO j = 1, nv - 1
402 10848 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
403 : END DO
404 : END DO
405 : END IF
406 :
407 858 : broy_w0 = mixing_store%broy_w0
408 858 : mixing_store%wbroy(nv) = 1.0_dp
409 :
410 : ! dfbroy
411 46430 : mixing_store%dfbroy(:, :, nv) = 0.0_dp
412 23374 : mixing_store%dfbroy(:, 1:ns, nv) = dq_now(:, 1:ns) - dq_last(:, 1:ns)
413 12116 : wdf = SUM(mixing_store%dfbroy(:, 1:ns, nv)**2)
414 858 : CALL para_env%sum(wdf)
415 858 : wdf = 1.0_dp/SQRT(wdf)
416 12116 : mixing_store%dfbroy(:, 1:ns, nv) = wdf*mixing_store%dfbroy(:, 1:ns, nv)
417 :
418 : ! abroy matrix
419 4814 : DO i = 1, nv
420 56234 : wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*mixing_store%dfbroy(:, 1:ns, nv))
421 3956 : CALL para_env%sum(wprod)
422 3956 : mixing_store%abroy(i, nv) = wprod
423 4814 : mixing_store%abroy(nv, i) = wprod
424 : END DO
425 :
426 : ! broyden matrices
427 7722 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
428 4814 : DO i = 1, nv
429 56234 : wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*dq_now(:, 1:ns))
430 3956 : CALL para_env%sum(wprod)
431 4814 : cvec(i) = mixing_store%wbroy(i)*wprod
432 : END DO
433 :
434 4814 : DO i = 1, nv
435 27240 : DO j = 1, nv
436 27240 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
437 : END DO
438 4814 : beta(i, i) = beta(i, i) + broy_w0
439 : END DO
440 :
441 858 : rskip = 1.e-12_dp
442 858 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
443 32054 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
444 :
445 : ! build ubroy
446 46430 : mixing_store%ubroy(:, :, nv) = 0.0_dp
447 : mixing_store%ubroy(:, 1:ns, nv) = alpha*mixing_store%dfbroy(:, 1:ns, nv) + &
448 23374 : wdf*(q_now(:, 1:ns) - q_last(:, 1:ns))
449 :
450 20018 : charges = 0.0_dp
451 3614 : DO ia = 1, mixing_store%nat_local
452 2756 : ii = mixing_store%atlist(ia)
453 11516 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
454 : END DO
455 4814 : DO i = 1, nv
456 17776 : DO ia = 1, mixing_store%nat_local
457 12962 : ii = mixing_store%atlist(ia)
458 53308 : charges(ii, 1:ns) = charges(ii, 1:ns) - mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
459 : END DO
460 : END DO
461 39178 : CALL para_env%sum(charges)
462 :
463 858 : DEALLOCATE (amat, beta, cvec, gammab)
464 :
465 858 : END SUBROUTINE broyden_mixing
466 :
467 : END MODULE qs_charge_mixing
|
