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_scf, ONLY: new_mixer
14 : #endif
15 : USE input_constants, ONLY: tblite_scc_mixer_auto, &
16 : tblite_scc_mixer_cp2k, &
17 : tblite_scc_mixer_none, &
18 : tblite_scc_mixer_tblite
19 : USE kinds, ONLY: dp
20 : USE mathlib, ONLY: get_pseudo_inverse_svd
21 : USE message_passing, ONLY: mp_para_env_type
22 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr, &
23 : gspace_mixing_nr, &
24 : mixing_storage_type, &
25 : multisecant_mixing_nr, &
26 : pulay_mixing_nr
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_charge_mixing'
34 :
35 : PUBLIC :: charge_mixing, charge_mixing_scc_error
36 :
37 : REAL(KIND=dp), PARAMETER, PRIVATE :: tblite_scc_pconv = 2.0E-5_dp
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief Driver for TB SCC variable mixing, calls the requested method.
43 : !> \param mixing_method ...
44 : !> \param mixing_store ...
45 : !> \param charges ...
46 : !> \param para_env ...
47 : !> \param iter_count ...
48 : !> \param scc_mixer ...
49 : !> \param tblite_mixer_damping ...
50 : !> \param tblite_mixer_memory ...
51 : !> \par History
52 : !> \author JGH
53 : ! **************************************************************************************************
54 26741 : SUBROUTINE charge_mixing(mixing_method, mixing_store, charges, para_env, iter_count, &
55 : scc_mixer, tblite_mixer_damping, tblite_mixer_memory)
56 : INTEGER, INTENT(IN) :: mixing_method
57 : TYPE(mixing_storage_type), POINTER :: mixing_store
58 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
59 : TYPE(mp_para_env_type), POINTER :: para_env
60 : INTEGER, INTENT(IN) :: iter_count
61 : INTEGER, INTENT(IN), OPTIONAL :: scc_mixer
62 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: tblite_mixer_damping
63 : INTEGER, INTENT(IN), OPTIONAL :: tblite_mixer_memory
64 :
65 : CHARACTER(len=*), PARAMETER :: routineN = 'charge_mixing'
66 :
67 : INTEGER :: effective_scc_mixer, handle, ia, ii, &
68 : imin, inow, mixer_memory, nbuffer, ns, &
69 : nvec
70 : REAL(dp) :: alpha, mixer_damping
71 :
72 26741 : CALL timeset(routineN, handle)
73 :
74 26741 : effective_scc_mixer = tblite_scc_mixer_cp2k
75 26741 : IF (PRESENT(scc_mixer)) effective_scc_mixer = scc_mixer
76 26741 : IF (ASSOCIATED(mixing_store)) mixing_store%tb_scc_mixer_error = 0.0_dp
77 :
78 495 : SELECT CASE (effective_scc_mixer)
79 : CASE (tblite_scc_mixer_auto, tblite_scc_mixer_cp2k)
80 : ! Use the regular CP2K SCC-variable mixing path below.
81 : CASE (tblite_scc_mixer_tblite)
82 495 : CPASSERT(ASSOCIATED(mixing_store))
83 495 : mixer_damping = 0.4_dp
84 495 : IF (PRESENT(tblite_mixer_damping)) mixer_damping = tblite_mixer_damping
85 495 : mixer_memory = MAX(1, mixing_store%nbuffer)
86 495 : IF (PRESENT(tblite_mixer_memory)) mixer_memory = MAX(1, tblite_mixer_memory)
87 : CALL tblite_charge_mixing(mixing_store, charges, para_env, iter_count, &
88 495 : mixer_damping, mixer_memory)
89 495 : CALL timestop(handle)
90 495 : RETURN
91 : CASE (tblite_scc_mixer_none)
92 0 : IF (ASSOCIATED(mixing_store)) mixing_store%iter_method = "NoMix"
93 0 : CALL timestop(handle)
94 0 : RETURN
95 : CASE DEFAULT
96 26741 : CPABORT("Unknown SCC mixer for TB charge mixing")
97 : END SELECT
98 :
99 26246 : IF (mixing_method >= gspace_mixing_nr) THEN
100 608 : CPASSERT(ASSOCIATED(mixing_store))
101 608 : mixing_store%ncall = mixing_store%ncall + 1
102 608 : ns = SIZE(charges, 2)
103 608 : IF (ns > mixing_store%max_shell) THEN
104 0 : CPABORT("Mixing storage too small for TB SCC variables")
105 : END IF
106 608 : alpha = mixing_store%alpha
107 608 : nbuffer = mixing_store%nbuffer
108 608 : inow = MOD(mixing_store%ncall - 1, nbuffer) + 1
109 608 : imin = inow - 1
110 608 : IF (imin == 0) imin = nbuffer
111 608 : IF (mixing_store%ncall > nbuffer) THEN
112 30 : nvec = nbuffer
113 : ELSE
114 578 : nvec = mixing_store%ncall - 1
115 : END IF
116 608 : IF (mixing_store%ncall > 1) THEN
117 : ! store in/out charge difference
118 2226 : DO ia = 1, mixing_store%nat_local
119 1702 : ii = mixing_store%atlist(ia)
120 6256 : mixing_store%dacharge(ia, 1:ns, imin) = mixing_store%acharge(ia, 1:ns, imin) - charges(ii, 1:ns)
121 : END DO
122 : END IF
123 608 : IF ((iter_count == 1) .OR. (iter_count + 1 <= mixing_store%nskip_mixing)) THEN
124 : ! skip mixing
125 84 : mixing_store%iter_method = "NoMix"
126 524 : ELSEIF (((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) .OR. (nvec == 1)) THEN
127 80 : CALL mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
128 80 : mixing_store%iter_method = "Mixing"
129 444 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
130 0 : CPABORT("Kerker method not available for Charge Mixing")
131 444 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
132 0 : CPABORT("Pulay method not available for Charge Mixing")
133 444 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
134 444 : CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
135 444 : mixing_store%iter_method = "Broy."
136 0 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
137 0 : CPABORT("Multisecant_mixing method not available for Charge Mixing")
138 : END IF
139 :
140 : ! store new 'input' charges
141 2560 : DO ia = 1, mixing_store%nat_local
142 1952 : ii = mixing_store%atlist(ia)
143 7200 : mixing_store%acharge(ia, 1:ns, inow) = charges(ii, 1:ns)
144 : END DO
145 :
146 : END IF
147 :
148 26246 : CALL timestop(handle)
149 :
150 : END SUBROUTINE charge_mixing
151 :
152 : ! **************************************************************************************************
153 : !> \brief Return the tblite SCC-mixer residual on the CP2K EPS_SCF scale.
154 : !> \param mixing_store ...
155 : !> \param eps_scf ...
156 : !> \return ...
157 : ! **************************************************************************************************
158 36950 : FUNCTION charge_mixing_scc_error(mixing_store, eps_scf) RESULT(mixer_error)
159 : TYPE(mixing_storage_type), POINTER :: mixing_store
160 : REAL(KIND=dp), INTENT(IN) :: eps_scf
161 : REAL(KIND=dp) :: mixer_error
162 :
163 36950 : mixer_error = 0.0_dp
164 36950 : IF (.NOT. ASSOCIATED(mixing_store)) RETURN
165 36950 : IF (mixing_store%tb_scc_mixer_step <= 1) RETURN
166 :
167 489 : IF (eps_scf > 0.0_dp) THEN
168 489 : mixer_error = eps_scf*mixing_store%tb_scc_mixer_error/tblite_scc_pconv
169 : ELSE
170 0 : mixer_error = mixing_store%tb_scc_mixer_error
171 : END IF
172 :
173 : END FUNCTION charge_mixing_scc_error
174 :
175 : ! **************************************************************************************************
176 : !> \brief TBLite modified-Broyden mixing for a complete TB SCC-variable vector.
177 : !> \param mixing_store ...
178 : !> \param charges ...
179 : !> \param para_env ...
180 : !> \param iter_count ...
181 : !> \param damping ...
182 : !> \param memory ...
183 : ! **************************************************************************************************
184 495 : SUBROUTINE tblite_charge_mixing(mixing_store, charges, para_env, iter_count, damping, memory)
185 : TYPE(mixing_storage_type), POINTER :: mixing_store
186 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
187 : TYPE(mp_para_env_type), POINTER :: para_env
188 : INTEGER, INTENT(IN) :: iter_count, memory
189 : REAL(KIND=dp), INTENT(IN) :: damping
190 :
191 : #if defined(__TBLITE)
192 495 : TYPE(error_type), ALLOCATABLE :: error
193 : #endif
194 : INTEGER :: natom, ndim, ns
195 : LOGICAL :: reset_mixer
196 495 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: qvec
197 :
198 : MARK_USED(para_env)
199 :
200 495 : natom = SIZE(charges, 1)
201 495 : ns = SIZE(charges, 2)
202 495 : ndim = natom*ns
203 1485 : ALLOCATE (qvec(ndim))
204 990 : qvec(:) = RESHAPE(charges, [ndim])
205 : reset_mixer = (iter_count == 1) .OR. (mixing_store%tb_scc_mixer_step == 0) .OR. &
206 : (mixing_store%tb_scc_mixer_natom /= natom) .OR. &
207 : (mixing_store%tb_scc_mixer_ns /= ns) .OR. &
208 495 : (mixing_store%tb_scc_mixer_memory /= memory)
209 495 : mixing_store%tb_scc_mixer_error = 0.0_dp
210 :
211 : #if defined(__TBLITE)
212 495 : IF (reset_mixer) THEN
213 3 : IF (ALLOCATED(mixing_store%tb_scc_mixer)) DEALLOCATE (mixing_store%tb_scc_mixer)
214 3 : CALL new_mixer(mixing_store%tb_scc_mixer, memory, ndim, damping)
215 3 : CALL mixing_store%tb_scc_mixer%set(qvec)
216 3 : mixing_store%tb_scc_mixer_natom = natom
217 3 : mixing_store%tb_scc_mixer_ns = ns
218 3 : mixing_store%tb_scc_mixer_memory = memory
219 3 : mixing_store%tb_scc_mixer_step = 1
220 3 : mixing_store%iter_method = "NoMix"
221 3 : RETURN
222 : END IF
223 :
224 492 : CPASSERT(ALLOCATED(mixing_store%tb_scc_mixer))
225 492 : CALL mixing_store%tb_scc_mixer%diff(qvec)
226 492 : mixing_store%tb_scc_mixer_error = REAL(mixing_store%tb_scc_mixer%get_error(), KIND=dp)
227 492 : CALL mixing_store%tb_scc_mixer%next(error)
228 492 : IF (ALLOCATED(error)) CPABORT("tblite SCC mixer failed")
229 492 : CALL mixing_store%tb_scc_mixer%get(qvec)
230 1476 : charges = RESHAPE(qvec, SHAPE(charges))
231 492 : mixing_store%tb_scc_mixer_step = mixing_store%tb_scc_mixer_step + 1
232 492 : mixing_store%iter_method = "TBLITE"
233 : #else
234 : MARK_USED(mixing_store)
235 : MARK_USED(charges)
236 : MARK_USED(iter_count)
237 : MARK_USED(damping)
238 : MARK_USED(memory)
239 : CPABORT("SCC_MIXER TBLITE requires CP2K to be built with tblite")
240 : #endif
241 :
242 495 : END SUBROUTINE tblite_charge_mixing
243 :
244 : ! **************************************************************************************************
245 : !> \brief Simple charge mixing
246 : !> \param mixing_store ...
247 : !> \param charges ...
248 : !> \param alpha ...
249 : !> \param imin ...
250 : !> \param ns ...
251 : !> \param para_env ...
252 : !> \author JGH
253 : ! **************************************************************************************************
254 80 : SUBROUTINE mix_charges_only(mixing_store, charges, alpha, imin, ns, para_env)
255 : TYPE(mixing_storage_type), POINTER :: mixing_store
256 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
257 : REAL(KIND=dp), INTENT(IN) :: alpha
258 : INTEGER, INTENT(IN) :: imin, ns
259 : TYPE(mp_para_env_type), POINTER :: para_env
260 :
261 : INTEGER :: ia, ii
262 :
263 1480 : charges = 0.0_dp
264 :
265 324 : DO ia = 1, mixing_store%nat_local
266 244 : ii = mixing_store%atlist(ia)
267 904 : charges(ii, 1:ns) = alpha*mixing_store%dacharge(ia, 1:ns, imin) - mixing_store%acharge(ia, 1:ns, imin)
268 : END DO
269 :
270 2880 : CALL para_env%sum(charges)
271 :
272 80 : END SUBROUTINE mix_charges_only
273 :
274 : ! **************************************************************************************************
275 : !> \brief Broyden charge mixing
276 : !> \param mixing_store ...
277 : !> \param charges ...
278 : !> \param inow ...
279 : !> \param nvec ...
280 : !> \param ns ...
281 : !> \param para_env ...
282 : !> \author JGH
283 : ! **************************************************************************************************
284 444 : SUBROUTINE broyden_mixing(mixing_store, charges, inow, nvec, ns, para_env)
285 : TYPE(mixing_storage_type), POINTER :: mixing_store
286 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: charges
287 : INTEGER, INTENT(IN) :: inow, nvec, ns
288 : TYPE(mp_para_env_type), POINTER :: para_env
289 :
290 : INTEGER :: i, ia, ii, imin, j, nbuffer, nv
291 : REAL(KIND=dp) :: alpha, maxw, minw, omega0, rskip, wdf, &
292 : wfac
293 444 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cvec, gammab
294 444 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, beta
295 444 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dq_last, dq_now, q_last, q_now
296 :
297 0 : CPASSERT(nvec > 1)
298 :
299 444 : nbuffer = mixing_store%nbuffer
300 444 : alpha = mixing_store%alpha
301 444 : imin = inow - 1
302 444 : IF (imin == 0) imin = nvec
303 444 : nv = nvec - 1
304 :
305 : ! charge vectors
306 444 : q_now => mixing_store%acharge(:, :, inow)
307 444 : q_last => mixing_store%acharge(:, :, imin)
308 444 : dq_now => mixing_store%dacharge(:, :, inow)
309 444 : dq_last => mixing_store%dacharge(:, :, imin)
310 :
311 444 : IF (nvec == nbuffer) THEN
312 : ! reshuffel Broyden storage n->n-1
313 146 : DO i = 1, nv - 1
314 116 : mixing_store%wbroy(i) = mixing_store%wbroy(i + 1)
315 6612 : mixing_store%dfbroy(:, :, i) = mixing_store%dfbroy(:, :, i + 1)
316 6642 : mixing_store%ubroy(:, :, i) = mixing_store%ubroy(:, :, i + 1)
317 : END DO
318 146 : DO i = 1, nv - 1
319 810 : DO j = 1, nv - 1
320 780 : mixing_store%abroy(i, j) = mixing_store%abroy(i + 1, j + 1)
321 : END DO
322 : END DO
323 : END IF
324 :
325 444 : omega0 = 0.01_dp
326 444 : minw = 1.0_dp
327 444 : maxw = 100000.0_dp
328 444 : wfac = 0.01_dp
329 :
330 5154 : mixing_store%wbroy(nv) = SUM(dq_now(:, 1:ns)**2)
331 444 : CALL para_env%sum(mixing_store%wbroy(nv))
332 444 : mixing_store%wbroy(nv) = SQRT(mixing_store%wbroy(nv))
333 444 : IF (mixing_store%wbroy(nv) > (wfac/maxw)) THEN
334 354 : mixing_store%wbroy(nv) = wfac/mixing_store%wbroy(nv)
335 : ELSE
336 90 : mixing_store%wbroy(nv) = maxw
337 : END IF
338 444 : IF (mixing_store%wbroy(nv) < minw) mixing_store%wbroy(nv) = minw
339 :
340 : ! dfbroy
341 11472 : mixing_store%dfbroy(:, :, nv) = 0.0_dp
342 9864 : mixing_store%dfbroy(:, 1:ns, nv) = dq_now(:, 1:ns) - dq_last(:, 1:ns)
343 5154 : wdf = SUM(mixing_store%dfbroy(:, 1:ns, nv)**2)
344 444 : CALL para_env%sum(wdf)
345 444 : wdf = 1.0_dp/SQRT(wdf)
346 5154 : mixing_store%dfbroy(:, 1:ns, nv) = wdf*mixing_store%dfbroy(:, 1:ns, nv)
347 :
348 : ! abroy matrix
349 2134 : DO i = 1, nv
350 22640 : wfac = SUM(mixing_store%dfbroy(:, 1:ns, i)*mixing_store%dfbroy(:, 1:ns, nv))
351 1690 : CALL para_env%sum(wfac)
352 1690 : mixing_store%abroy(i, nv) = wfac
353 2134 : mixing_store%abroy(nv, i) = wfac
354 : END DO
355 :
356 : ! broyden matrices
357 3996 : ALLOCATE (amat(nv, nv), beta(nv, nv), cvec(nv), gammab(nv))
358 2134 : DO i = 1, nv
359 22640 : wfac = SUM(mixing_store%dfbroy(:, 1:ns, i)*dq_now(:, 1:ns))
360 1690 : CALL para_env%sum(wfac)
361 2134 : cvec(i) = mixing_store%wbroy(i)*wfac
362 : END DO
363 :
364 2134 : DO i = 1, nv
365 11380 : DO j = 1, nv
366 11380 : beta(j, i) = mixing_store%wbroy(j)*mixing_store%wbroy(i)*mixing_store%abroy(j, i)
367 : END DO
368 2134 : beta(i, i) = beta(i, i) + omega0*omega0
369 : END DO
370 :
371 444 : rskip = 1.e-12_dp
372 444 : CALL get_pseudo_inverse_svd(beta, amat, rskip)
373 13514 : gammab(1:nv) = MATMUL(cvec(1:nv), amat(1:nv, 1:nv))
374 :
375 : ! build ubroy
376 11472 : mixing_store%ubroy(:, :, nv) = 0.0_dp
377 : mixing_store%ubroy(:, 1:ns, nv) = alpha*mixing_store%dfbroy(:, 1:ns, nv) + &
378 9864 : wdf*(q_now(:, 1:ns) - q_last(:, 1:ns))
379 :
380 8604 : charges = 0.0_dp
381 1902 : DO ia = 1, mixing_store%nat_local
382 1458 : ii = mixing_store%atlist(ia)
383 5352 : charges(ii, 1:ns) = q_now(ia, 1:ns) + alpha*dq_now(ia, 1:ns)
384 : END DO
385 2134 : DO i = 1, nv
386 7994 : DO ia = 1, mixing_store%nat_local
387 5860 : ii = mixing_store%atlist(ia)
388 23410 : charges(ii, 1:ns) = charges(ii, 1:ns) - mixing_store%wbroy(i)*gammab(i)*mixing_store%ubroy(ia, 1:ns, i)
389 : END DO
390 : END DO
391 16764 : CALL para_env%sum(charges)
392 :
393 444 : DEALLOCATE (amat, beta, cvec, gammab)
394 :
395 444 : END SUBROUTINE broyden_mixing
396 :
397 : END MODULE qs_charge_mixing
|
