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 39592 : 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 39592 : CALL timeset(routineN, handle)
97 :
98 39592 : effective_scc_mixer = tblite_scc_mixer_cp2k
99 39592 : IF (PRESENT(scc_mixer)) effective_scc_mixer = scc_mixer
100 39592 : 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 39592 : CPABORT("Unknown SCC mixer for TB charge mixing")
158 : END SELECT
159 :
160 39560 : IF (mixing_method >= gspace_mixing_nr) THEN
161 1074 : CPASSERT(ASSOCIATED(mixing_store))
162 1074 : mixing_store%ncall = mixing_store%ncall + 1
163 1074 : ns = SIZE(charges, 2)
164 1074 : IF (ns > mixing_store%max_shell) THEN
165 0 : CPABORT("Mixing storage too small for TB SCC variables")
166 : END IF
167 1074 : alpha = mixing_store%alpha
168 1074 : nbuffer = mixing_store%nbuffer
169 1074 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
170 1074 : imin = inow - 1
171 1074 : IF (imin == 0) imin = nbuffer
172 1074 : IF (mixing_store%ncall > nbuffer) THEN
173 286 : nvec = nbuffer
174 : ELSE
175 788 : nvec = mixing_store%ncall - 1
176 : END IF
177 1074 : IF (mixing_store%ncall > 1) THEN
178 : ! store in/out charge difference
179 4036 : DO ia = 1, mixing_store%nat_local
180 3080 : ii = mixing_store%atlist(ia)
181 12680 : 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 1074 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
185 : ! skip mixing
186 128 : mixing_store%iter_method = "NoMix"
187 946 : ELSE IF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
188 102 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
189 102 : 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 844 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
196 844 : mixing_store%iter_method = "Broy."
197 : ELSE IF (mixing_method == modified_broyden_mixing_nr) THEN
198 0 : CPABORT("Modified Broyden mixing is only available for DFT density mixing")
199 : ELSE IF (mixing_method == multisecant_mixing_nr) THEN
200 0 : CPABORT("Multisecant_mixing method not available for Charge Mixing")
201 : ELSE IF (mixing_method == new_pulay_mixing_nr) THEN
202 0 : CPABORT("New Pulay method not available for Charge Mixing")
203 : END IF
204 :
205 : ! store new 'input' charges
206 4526 : DO ia = 1, mixing_store%nat_local
207 3452 : ii = mixing_store%atlist(ia)
208 14424 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
209 : END DO
210 :
211 : END IF
212 :
213 39560 : CALL timestop(handle)
214 :
215 : END SUBROUTINE charge_mixing
216 :
217 : ! **************************************************************************************************
218 : !> \brief Return the CP2K-side tblite SCC-mixer residual on the CP2K EPS_SCF scale.
219 : !> \param mixing_store ...
220 : !> \param eps_scf ...
221 : !> \return ...
222 : ! **************************************************************************************************
223 54844 : FUNCTION charge_mixing_scc_error(mixing_store, eps_scf) RESULT(mixer_error)
224 : TYPE(mixing_storage_type), POINTER :: mixing_store
225 : REAL(KIND=dp), INTENT(IN) :: eps_scf
226 : REAL(KIND=dp) :: mixer_error
227 :
228 54844 : mixer_error = 0.0_dp
229 54844 : IF (.NOT. ASSOCIATED(mixing_store)) RETURN
230 54844 : IF (mixing_store%tb_scc_mixer_step <= 1) RETURN
231 :
232 26 : IF (eps_scf > 0.0_dp) THEN
233 26 : mixer_error = eps_scf*mixing_store%tb_scc_mixer_error/tblite_scc_pconv
234 : ELSE
235 0 : mixer_error = mixing_store%tb_scc_mixer_error
236 : END IF
237 :
238 : END FUNCTION charge_mixing_scc_error
239 :
240 : ! **************************************************************************************************
241 : !> \brief TBLite modified-Broyden mixing for a complete TB SCC-variable vector.
242 : !> \param mixing_store ...
243 : !> \param charges ...
244 : !> \param para_env ...
245 : !> \param iter_count ...
246 : !> \param damping ...
247 : !> \param memory ...
248 : !> \param omega0 ...
249 : !> \param min_weight ...
250 : !> \param max_weight ...
251 : !> \param weight_factor ...
252 : ! **************************************************************************************************
253 32 : SUBROUTINE tblite_charge_mixing(mixing_store, charges, para_env, iter_count, damping, memory, omega0, &
254 : min_weight, max_weight, weight_factor)
255 : TYPE(mixing_storage_type), POINTER :: mixing_store
256 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
257 : TYPE(mp_para_env_type), POINTER :: para_env
258 : INTEGER, INTENT(IN) :: iter_count, memory
259 : REAL(KIND=dp), INTENT(IN) :: damping, max_weight, min_weight, omega0, &
260 : weight_factor
261 :
262 : #if defined(__TBLITE)
263 32 : TYPE(error_type), ALLOCATABLE :: error
264 : #endif
265 : INTEGER :: natom, ndim, ns
266 : LOGICAL :: on_source, reset_mixer
267 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: qvec
268 :
269 32 : natom = SIZE(charges, 1)
270 32 : ns = SIZE(charges, 2)
271 32 : ndim = natom*ns
272 96 : ALLOCATE (qvec(ndim))
273 64 : qvec(:) = RESHAPE(charges, [ndim])
274 32 : on_source = para_env%mepos == para_env%source
275 : reset_mixer = (iter_count == 1) .OR. (mixing_store%tb_scc_mixer_step == 0) .OR. &
276 : (mixing_store%tb_scc_mixer_natom /= natom) .OR. &
277 : (mixing_store%tb_scc_mixer_ns /= ns) .OR. &
278 32 : (mixing_store%tb_scc_mixer_memory /= memory)
279 32 : mixing_store%tb_scc_mixer_error = 0.0_dp
280 :
281 : #if defined(__TBLITE)
282 32 : IF (reset_mixer) THEN
283 4 : IF (ALLOCATED(mixing_store%tb_scc_mixer)) DEALLOCATE (mixing_store%tb_scc_mixer)
284 4 : IF (on_source) THEN
285 : CALL new_cp2k_tblite_mixer(mixing_store%tb_scc_mixer, memory, ndim, damping, omega0, &
286 2 : min_weight, max_weight, weight_factor)
287 2 : CALL mixing_store%tb_scc_mixer%set(qvec)
288 : END IF
289 4 : mixing_store%tb_scc_mixer_natom = natom
290 4 : mixing_store%tb_scc_mixer_ns = ns
291 4 : mixing_store%tb_scc_mixer_memory = memory
292 4 : mixing_store%tb_scc_mixer_step = 1
293 4 : mixing_store%iter_method = "NoMix"
294 4 : CALL para_env%bcast(qvec)
295 12 : charges = RESHAPE(qvec, SHAPE(charges))
296 4 : RETURN
297 : END IF
298 :
299 28 : IF (on_source) THEN
300 14 : CPASSERT(ALLOCATED(mixing_store%tb_scc_mixer))
301 14 : CALL mixing_store%tb_scc_mixer%diff(qvec)
302 14 : mixing_store%tb_scc_mixer_error = REAL(mixing_store%tb_scc_mixer%get_error(), KIND=dp)
303 14 : CALL mixing_store%tb_scc_mixer%next(error)
304 14 : IF (ALLOCATED(error)) CPABORT("tblite SCC mixer failed")
305 14 : CALL mixing_store%tb_scc_mixer%get(qvec)
306 : END IF
307 28 : CALL para_env%bcast(qvec)
308 28 : CALL para_env%bcast(mixing_store%tb_scc_mixer_error)
309 84 : charges = RESHAPE(qvec, SHAPE(charges))
310 28 : mixing_store%tb_scc_mixer_step = mixing_store%tb_scc_mixer_step + 1
311 28 : mixing_store%iter_method = "TBLITE"
312 : #else
313 : MARK_USED(mixing_store)
314 : MARK_USED(charges)
315 : MARK_USED(para_env)
316 : MARK_USED(iter_count)
317 : MARK_USED(damping)
318 : MARK_USED(memory)
319 : MARK_USED(omega0)
320 : MARK_USED(min_weight)
321 : MARK_USED(max_weight)
322 : MARK_USED(weight_factor)
323 : CPABORT("SCC_MIXER TBLITE requires CP2K to be built with tblite")
324 : #endif
325 :
326 32 : END SUBROUTINE tblite_charge_mixing
327 :
328 : ! **************************************************************************************************
329 : !> \brief Simple charge mixing
330 : !> \param mixing_store ...
331 : !> \param charges ...
332 : !> \param alpha ...
333 : !> \param imin ...
334 : !> \param ns ...
335 : !> \param para_env ...
336 : !> \author JGH
337 : ! **************************************************************************************************
338 102 : SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
339 : TYPE(mixing_storage_type), POINTER :: mixing_store
340 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
341 : REAL(KIND=dp), INTENT(IN) :: alpha
342 : INTEGER, INTENT(IN) :: imin, ns
343 : TYPE(mp_para_env_type), POINTER :: para_env
344 :
345 : INTEGER :: ia, ii
346 :
347 2372 : charges = 0.0_dp
348 :
349 422 : DO ia = 1, mixing_store%nat_local
350 320 : ii = mixing_store%atlist(ia)
351 1382 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
352 : END DO
353 :
354 4642 : CALL para_env%sum(charges)
355 :
356 102 : END SUBROUTINE mix_charges_only
357 :
358 : ! **************************************************************************************************
359 : !> \brief Broyden charge mixing
360 : !> \param mixing_store ...
361 : !> \param charges ...
362 : !> \param inow ...
363 : !> \param nvec ...
364 : !> \param ns ...
365 : !> \param para_env ...
366 : !> \author JGH
367 : ! **************************************************************************************************
368 844 : SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env)
369 : TYPE(mixing_storage_type), POINTER :: mixing_store
370 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
371 : INTEGER, INTENT(IN) :: inow, nvec, ns
372 : TYPE(mp_para_env_type), POINTER :: para_env
373 :
374 : INTEGER :: i, ia, ii, imin, j, nbuffer, nv
375 : REAL(KIND=dp) :: alpha, broy_w0, rskip, wdf, wprod
376 844 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
377 844 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
378 844 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
379 :
380 0 : CPASSERT(nvec > 1)
381 :
382 844 : nbuffer = mixing_store%nbuffer
383 844 : alpha = mixing_store%alpha
384 844 : imin = inow - 1
385 844 : IF (imin == 0) imin = nvec
386 844 : nv = nvec - 1
387 :
388 : ! charge vectors
389 844 : q_now => mixing_store%acharge(:, :, inow)
390 844 : q_last => mixing_store%acharge(:, :, imin)
391 844 : dq_now => mixing_store%dacharge(:, :, inow)
392 844 : dq_last => mixing_store%dacharge(:, :, imin)
393 :
394 844 : IF (nvec == nbuffer) THEN
395 : ! reshuffel Broyden storage n->n-1
396 1802 : DO i = 1, nv - 1
397 1516 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
398 39380 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
399 39666 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
400 : END DO
401 1802 : DO i = 1, nv - 1
402 10546 : DO j = 1, nv - 1
403 10260 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
404 : END DO
405 : END DO
406 : END IF
407 :
408 844 : broy_w0 = mixing_store%broy_w0
409 844 : mixing_store%wbroy(nv) = 1.0_dp
410 :
411 : ! dfbroy
412 45632 : mixing_store%dfbroy(:, :, nv) = 0.0_dp
413 21992 : mixing_store%dfbroy(:, 1:ns, nv) = dq_now(:, 1:ns) - dq_last(:, 1:ns)
414 11418 : wdf = SUM(mixing_store%dfbroy(:, 1:ns, nv)**2)
415 844 : CALL para_env%sum(wdf)
416 844 : wdf = 1.0_dp/SQRT(wdf)
417 11418 : mixing_store%dfbroy(:, 1:ns, nv) = wdf*mixing_store%dfbroy(:, 1:ns, nv)
418 :
419 : ! abroy matrix
420 4702 : DO i = 1, nv
421 51348 : wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*mixing_store%dfbroy(:, 1:ns, nv))
422 3858 : CALL para_env%sum(wprod)
423 3858 : mixing_store%abroy(i, nv) = wprod
424 4702 : mixing_store%abroy(nv, i) = wprod
425 : END DO
426 :
427 : ! broyden matrices
428 7596 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
429 4702 : DO i = 1, nv
430 51348 : wprod = SUM(mixing_store%dfbroy(:, 1:ns, i)*dq_now(:, 1:ns))
431 3858 : CALL para_env%sum(wprod)
432 4702 : cvec(i) = mixing_store%wbroy(i)*wprod
433 : END DO
434 :
435 4702 : DO i = 1, nv
436 26456 : DO j = 1, nv
437 26456 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
438 : END DO
439 4702 : beta(i, i) = beta(i, i) + broy_w0
440 : END DO
441 :
442 844 : rskip = 1.e-12_dp
443 844 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
444 31158 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
445 :
446 : ! build ubroy
447 45632 : mixing_store%ubroy(:, :, nv) = 0.0_dp
448 : mixing_store%ubroy(:, 1:ns, nv) = alpha*mixing_store%dfbroy(:, 1:ns, nv) + &
449 21992 : wdf*(q_now(:, 1:ns) - q_last(:, 1:ns))
450 :
451 18978 : charges = 0.0_dp
452 3586 : DO ia = 1, mixing_store%nat_local
453 2742 : ii = mixing_store%atlist(ia)
454 11146 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
455 : END DO
456 4702 : DO i = 1, nv
457 17566 : DO ia = 1, mixing_store%nat_local
458 12864 : ii = mixing_store%atlist(ia)
459 50718 : charges(ii, 1:ns) = charges(ii, 1:ns) - mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
460 : END DO
461 : END DO
462 37112 : CALL para_env%sum(charges)
463 :
464 844 : DEALLOCATE (amat, beta, cvec, gammab)
465 :
466 844 : END SUBROUTINE broyden_mixing
467 :
468 : END MODULE qs_charge_mixing
|
