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