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 : !> \brief Dense density mixing for active-space embedding.
10 : ! **************************************************************************************************
11 : MODULE qs_active_space_mixing
12 : USE cp_fm_types, ONLY: cp_fm_get_element,&
13 : cp_fm_set_element,&
14 : cp_fm_type
15 : USE input_section_types, ONLY: section_vals_get,&
16 : section_vals_get_subs_vals,&
17 : section_vals_type,&
18 : section_vals_val_get
19 : USE kinds, ONLY: dp
20 : USE mathlib, ONLY: invert_matrix
21 : USE qs_active_space_types, ONLY: active_space_type
22 : USE qs_density_mixing_types, ONLY: &
23 : broyden_mixing_nr, direct_mixing_nr, gspace_mixing_nr, mixing_storage_create, &
24 : mixing_storage_type, modified_broyden_mixing_nr, multisecant_mixing_nr, no_mixing_nr, &
25 : pulay_mixing_nr
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 : PRIVATE
30 :
31 : PUBLIC :: active_space_mixing_label
32 : PUBLIC :: initialize_active_space_mixing
33 : PUBLIC :: update_active_density
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief Initialize the density mixer used in the self-consistent active-space embedding loop.
39 : !> \param active_space_env active space environment
40 : !> \param as_input ACTIVE_SPACE input section
41 : ! **************************************************************************************************
42 328 : SUBROUTINE initialize_active_space_mixing(active_space_env, as_input)
43 : TYPE(active_space_type), POINTER :: active_space_env
44 : TYPE(section_vals_type), POINTER :: as_input
45 :
46 : LOGICAL :: do_mixing, legacy_alpha_explicit, &
47 : mixing_alpha_explicit, mixing_explicit
48 : REAL(KIND=dp) :: legacy_alpha
49 : TYPE(section_vals_type), POINTER :: mixing_section
50 :
51 82 : NULLIFY (mixing_section)
52 :
53 : CALL section_vals_val_get(as_input, "ALPHA", r_val=legacy_alpha, &
54 82 : explicit=legacy_alpha_explicit)
55 82 : IF (legacy_alpha < 0.0_dp .OR. legacy_alpha > 1.0_dp) THEN
56 0 : CPABORT("Specify an active-space damping factor between 0 and 1.")
57 : END IF
58 :
59 82 : mixing_section => section_vals_get_subs_vals(as_input, "MIXING")
60 82 : CALL section_vals_get(mixing_section, explicit=mixing_explicit)
61 82 : CALL section_vals_val_get(mixing_section, "_SECTION_PARAMETERS_", l_val=do_mixing)
62 :
63 82 : active_space_env%as_mixing_dim = 0
64 82 : active_space_env%as_mixing_iter = 0
65 :
66 82 : IF (.NOT. do_mixing) THEN
67 0 : active_space_env%as_mixing_method = no_mixing_nr
68 0 : active_space_env%alpha = 1.0_dp
69 0 : RETURN
70 : END IF
71 :
72 : CALL section_vals_val_get(mixing_section, "METHOD", &
73 82 : i_val=active_space_env%as_mixing_method)
74 :
75 82 : SELECT CASE (active_space_env%as_mixing_method)
76 : CASE (no_mixing_nr)
77 0 : active_space_env%alpha = 1.0_dp
78 0 : RETURN
79 : CASE (direct_mixing_nr, pulay_mixing_nr, broyden_mixing_nr, modified_broyden_mixing_nr)
80 0 : CONTINUE
81 : CASE (gspace_mixing_nr)
82 : CALL cp_abort(__LOCATION__, &
83 : "ACTIVE_SPACE%MIXING%METHOD KERKER_MIXING is not supported. "// &
84 : "The active-space density lives in the active MO subspace, "// &
85 0 : "not in G-space.")
86 : CASE (multisecant_mixing_nr)
87 0 : CPABORT("ACTIVE_SPACE%MIXING%METHOD MULTISECANT_MIXING is not yet supported.")
88 : CASE DEFAULT
89 82 : CPABORT("Unknown ACTIVE_SPACE%MIXING%METHOD.")
90 : END SELECT
91 :
92 82 : IF (ASSOCIATED(active_space_env%as_mixing_store)) THEN
93 0 : CPABORT("Active-space mixing storage already initialized.")
94 : END IF
95 246 : ALLOCATE (active_space_env%as_mixing_store)
96 : CALL mixing_storage_create(active_space_env%as_mixing_store, mixing_section, &
97 82 : active_space_env%as_mixing_method, ecut=0.0_dp)
98 :
99 82 : CALL section_vals_val_get(mixing_section, "ALPHA", explicit=mixing_alpha_explicit)
100 82 : IF ((legacy_alpha_explicit .AND. (.NOT. mixing_alpha_explicit)) .OR. &
101 : ((.NOT. mixing_explicit) .AND. (.NOT. mixing_alpha_explicit))) THEN
102 80 : active_space_env%as_mixing_store%alpha = legacy_alpha
103 : END IF
104 82 : IF (active_space_env%as_mixing_store%alpha < 0.0_dp .OR. &
105 : active_space_env%as_mixing_store%alpha > 1.0_dp) THEN
106 0 : CPABORT("Specify an active-space mixing ALPHA between 0 and 1.")
107 : END IF
108 82 : IF (active_space_env%as_mixing_store%nbuffer < 1 .AND. &
109 : active_space_env%as_mixing_method /= direct_mixing_nr) THEN
110 0 : CPABORT("ACTIVE_SPACE%MIXING%NBUFFER has to be positive.")
111 : END IF
112 :
113 82 : active_space_env%alpha = active_space_env%as_mixing_store%alpha
114 :
115 : END SUBROUTINE initialize_active_space_mixing
116 :
117 : ! **************************************************************************************************
118 : !> \brief Return the current active-space mixer label for iteration output.
119 : !> \param active_space_env active space environment
120 : !> \return short mixer label
121 : ! **************************************************************************************************
122 7 : FUNCTION active_space_mixing_label(active_space_env) RESULT(label)
123 : TYPE(active_space_type), POINTER :: active_space_env
124 : CHARACTER(len=15) :: label
125 :
126 14 : SELECT CASE (active_space_env%as_mixing_method)
127 : CASE (direct_mixing_nr)
128 7 : label = "P_Mix"
129 : CASE (pulay_mixing_nr)
130 0 : label = "Pulay"
131 : CASE (broyden_mixing_nr)
132 0 : label = "Broy."
133 : CASE (modified_broyden_mixing_nr)
134 0 : label = "MBroy"
135 : CASE DEFAULT
136 7 : label = "NoMix"
137 : END SELECT
138 7 : IF (ASSOCIATED(active_space_env%as_mixing_store)) THEN
139 7 : IF (active_space_env%as_mixing_iter > 0 .AND. &
140 : LEN_TRIM(active_space_env%as_mixing_store%iter_method) > 0) THEN
141 4 : label = active_space_env%as_mixing_store%iter_method
142 : END IF
143 : END IF
144 :
145 7 : END FUNCTION active_space_mixing_label
146 :
147 : ! **************************************************************************************************
148 : !> \brief Release dense active-space mixing history buffers.
149 : !> \param active_space_env active space environment
150 : ! **************************************************************************************************
151 0 : SUBROUTINE release_active_mixing_history(active_space_env)
152 : TYPE(active_space_type), POINTER :: active_space_env
153 :
154 0 : IF (ASSOCIATED(active_space_env%as_mix_r_old)) THEN
155 0 : DEALLOCATE (active_space_env%as_mix_r_old)
156 : END IF
157 0 : IF (ASSOCIATED(active_space_env%as_mix_weight)) THEN
158 0 : DEALLOCATE (active_space_env%as_mix_weight)
159 : END IF
160 0 : IF (ASSOCIATED(active_space_env%as_mix_x_old)) THEN
161 0 : DEALLOCATE (active_space_env%as_mix_x_old)
162 : END IF
163 0 : IF (ASSOCIATED(active_space_env%as_mix_r_buffer)) THEN
164 0 : DEALLOCATE (active_space_env%as_mix_r_buffer)
165 : END IF
166 0 : IF (ASSOCIATED(active_space_env%as_mix_x_buffer)) THEN
167 0 : DEALLOCATE (active_space_env%as_mix_x_buffer)
168 : END IF
169 0 : active_space_env%as_mixing_dim = 0
170 :
171 0 : END SUBROUTINE release_active_mixing_history
172 :
173 : ! **************************************************************************************************
174 : !> \brief Ensure dense active-space mixing history buffers are allocated.
175 : !> \param active_space_env active space environment
176 : !> \param ndim length of the flattened density vector
177 : ! **************************************************************************************************
178 0 : SUBROUTINE ensure_active_mixing_history(active_space_env, ndim)
179 : TYPE(active_space_type), POINTER :: active_space_env
180 : INTEGER, INTENT(IN) :: ndim
181 :
182 : INTEGER :: nbuffer
183 : TYPE(mixing_storage_type), POINTER :: mixing_store
184 :
185 0 : IF (.NOT. ASSOCIATED(active_space_env%as_mixing_store)) RETURN
186 0 : IF (active_space_env%as_mixing_method == direct_mixing_nr) RETURN
187 :
188 0 : mixing_store => active_space_env%as_mixing_store
189 0 : nbuffer = mixing_store%nbuffer
190 0 : IF (nbuffer < 1) CPABORT("ACTIVE_SPACE%MIXING%NBUFFER has to be positive.")
191 :
192 0 : IF (active_space_env%as_mixing_dim == ndim .AND. &
193 : ASSOCIATED(active_space_env%as_mix_r_buffer)) RETURN
194 :
195 0 : CALL release_active_mixing_history(active_space_env)
196 :
197 0 : active_space_env%as_mixing_dim = ndim
198 0 : ALLOCATE (active_space_env%as_mix_r_old(ndim))
199 0 : ALLOCATE (active_space_env%as_mix_weight(nbuffer))
200 0 : ALLOCATE (active_space_env%as_mix_x_old(ndim))
201 0 : ALLOCATE (active_space_env%as_mix_r_buffer(nbuffer, ndim))
202 0 : ALLOCATE (active_space_env%as_mix_x_buffer(nbuffer, ndim))
203 :
204 0 : active_space_env%as_mix_r_old = 0.0_dp
205 0 : active_space_env%as_mix_weight = 1.0_dp
206 0 : active_space_env%as_mix_x_old = 0.0_dp
207 0 : active_space_env%as_mix_r_buffer = 0.0_dp
208 0 : active_space_env%as_mix_x_buffer = 0.0_dp
209 0 : mixing_store%ncall = 0
210 :
211 : END SUBROUTINE ensure_active_mixing_history
212 :
213 : ! **************************************************************************************************
214 : !> \brief Apply a direct dense active-space density mix.
215 : !> \param p_old current active-space density vector
216 : !> \param p_solver active-space density vector returned by the external solver
217 : !> \param alpha mixing damping
218 : !> \param p_mixed mixed active-space density vector
219 : ! **************************************************************************************************
220 8 : SUBROUTINE active_space_direct_mix(p_old, p_solver, alpha, p_mixed)
221 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_old, p_solver
222 : REAL(KIND=dp), INTENT(IN) :: alpha
223 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: p_mixed
224 :
225 64 : p_mixed = p_old + alpha*(p_solver - p_old)
226 :
227 8 : END SUBROUTINE active_space_direct_mix
228 :
229 : ! **************************************************************************************************
230 : !> \brief Apply Pulay mixing to the dense active-space density vector.
231 : !> \param active_space_env active space environment
232 : !> \param p_old current active-space density vector
233 : !> \param p_solver active-space density vector returned by the external solver
234 : !> \param p_mixed mixed active-space density vector
235 : ! **************************************************************************************************
236 0 : SUBROUTINE active_space_pulay_mixing(active_space_env, p_old, p_solver, p_mixed)
237 : TYPE(active_space_type), POINTER :: active_space_env
238 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_old, p_solver
239 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: p_mixed
240 :
241 : INTEGER :: i, ib, ibb, j, nb, nbuffer, ndim
242 : REAL(KIND=dp) :: inv_err, norm_c_inv, res_norm
243 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, p_diis
244 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: c, c_inv
245 : TYPE(mixing_storage_type), POINTER :: mixing_store
246 :
247 0 : ndim = SIZE(p_old)
248 0 : CALL ensure_active_mixing_history(active_space_env, ndim)
249 0 : mixing_store => active_space_env%as_mixing_store
250 0 : nbuffer = mixing_store%nbuffer
251 :
252 0 : ib = MODULO(mixing_store%ncall, nbuffer) + 1
253 0 : mixing_store%ncall = mixing_store%ncall + 1
254 0 : nb = MIN(mixing_store%ncall, nbuffer)
255 0 : ibb = MODULO(mixing_store%ncall, nbuffer) + 1
256 :
257 0 : active_space_env%as_mix_x_buffer(ib, :) = p_old
258 0 : active_space_env%as_mix_r_buffer(ib, :) = p_solver - p_old
259 : res_norm = SQRT(DOT_PRODUCT(active_space_env%as_mix_r_buffer(ib, :), &
260 0 : active_space_env%as_mix_r_buffer(ib, :)))
261 :
262 0 : IF (nb == 1 .OR. res_norm < 1.E-14_dp) THEN
263 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
264 : ELSE
265 0 : ALLOCATE (c(nb, nb))
266 0 : ALLOCATE (c_inv(nb, nb))
267 0 : ALLOCATE (alpha_c(nb))
268 0 : ALLOCATE (p_diis(ndim))
269 :
270 0 : c(:, :) = 0.0_dp
271 0 : DO i = 1, nb
272 0 : DO j = i, nb
273 : c(j, i) = DOT_PRODUCT(active_space_env%as_mix_r_buffer(i, :), &
274 0 : active_space_env%as_mix_r_buffer(j, :))
275 0 : c(i, j) = c(j, i)
276 : END DO
277 : END DO
278 :
279 0 : CALL invert_matrix(c, c_inv, inv_err, improve=.TRUE.)
280 0 : norm_c_inv = SUM(c_inv)
281 0 : IF (ABS(norm_c_inv) < 1.E-14_dp) THEN
282 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
283 : ELSE
284 0 : DO i = 1, nb
285 0 : alpha_c(i) = SUM(c_inv(:, i))/norm_c_inv
286 : END DO
287 :
288 0 : p_diis(:) = 0.0_dp
289 0 : DO i = 1, nb
290 : p_diis(:) = p_diis(:) + alpha_c(i)*(active_space_env%as_mix_x_buffer(i, :) + &
291 0 : mixing_store%pulay_beta*active_space_env%as_mix_r_buffer(i, :))
292 : END DO
293 0 : IF (mixing_store%pulay_alpha > 0.0_dp) THEN
294 : p_mixed = mixing_store%pulay_alpha*p_solver + &
295 0 : (1.0_dp - mixing_store%pulay_alpha)*p_diis
296 : ELSE
297 0 : p_mixed = p_diis
298 : END IF
299 : END IF
300 :
301 0 : DEALLOCATE (alpha_c)
302 0 : DEALLOCATE (p_diis)
303 0 : DEALLOCATE (c)
304 0 : DEALLOCATE (c_inv)
305 : END IF
306 :
307 0 : active_space_env%as_mix_x_buffer(ibb, :) = p_mixed
308 0 : mixing_store%iter_method = "Pulay"
309 :
310 0 : END SUBROUTINE active_space_pulay_mixing
311 :
312 : ! **************************************************************************************************
313 : !> \brief Apply original or modified Broyden mixing to the dense active-space density vector.
314 : !> \param active_space_env active space environment
315 : !> \param p_old current active-space density vector
316 : !> \param p_solver active-space density vector returned by the external solver
317 : !> \param p_mixed mixed active-space density vector
318 : !> \param modified use dynamic residual weights of modified Broyden
319 : ! **************************************************************************************************
320 0 : SUBROUTINE active_space_broyden_mixing(active_space_env, p_old, p_solver, p_mixed, modified)
321 : TYPE(active_space_type), POINTER :: active_space_env
322 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_old, p_solver
323 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: p_mixed
324 : LOGICAL, INTENT(IN) :: modified
325 :
326 : INTEGER :: ib, j, k, nb, nbuffer, ndim
327 : LOGICAL :: can_update
328 : REAL(KIND=dp) :: delta_norm, inv_err, res_norm, weight
329 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c, g, p_res
330 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
331 : TYPE(mixing_storage_type), POINTER :: mixing_store
332 :
333 0 : ndim = SIZE(p_old)
334 0 : CALL ensure_active_mixing_history(active_space_env, ndim)
335 0 : mixing_store => active_space_env%as_mixing_store
336 0 : nbuffer = mixing_store%nbuffer
337 :
338 0 : ALLOCATE (p_res(ndim))
339 0 : p_res(:) = p_solver(:) - p_old(:)
340 0 : res_norm = SQRT(DOT_PRODUCT(p_res, p_res))
341 :
342 0 : mixing_store%ncall = mixing_store%ncall + 1
343 0 : IF (mixing_store%ncall == 1) THEN
344 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
345 0 : active_space_env%as_mix_x_old = p_old
346 0 : active_space_env%as_mix_r_old = p_res
347 0 : IF (modified) THEN
348 0 : mixing_store%iter_method = "MBroy"
349 : ELSE
350 0 : mixing_store%iter_method = "Broy."
351 : END IF
352 0 : DEALLOCATE (p_res)
353 0 : RETURN
354 : END IF
355 :
356 0 : nb = MIN(mixing_store%ncall - 1, nbuffer)
357 0 : ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1
358 :
359 0 : active_space_env%as_mix_r_buffer(ib, :) = p_res - active_space_env%as_mix_r_old
360 0 : active_space_env%as_mix_x_buffer(ib, :) = p_old - active_space_env%as_mix_x_old
361 : delta_norm = SQRT(DOT_PRODUCT(active_space_env%as_mix_r_buffer(ib, :), &
362 0 : active_space_env%as_mix_r_buffer(ib, :)))
363 0 : can_update = res_norm > 1.E-14_dp .AND. delta_norm > 1.E-14_dp
364 :
365 : IF (can_update) THEN
366 0 : active_space_env%as_mix_r_buffer(ib, :) = active_space_env%as_mix_r_buffer(ib, :)/delta_norm
367 : active_space_env%as_mix_x_buffer(ib, :) = active_space_env%as_mix_x_buffer(ib, :)/delta_norm + &
368 0 : mixing_store%alpha*active_space_env%as_mix_r_buffer(ib, :)
369 :
370 0 : IF (modified) THEN
371 0 : IF (res_norm > (mixing_store%wc/mixing_store%wmax)) THEN
372 0 : active_space_env%as_mix_weight(ib) = mixing_store%wc/res_norm
373 : ELSE
374 0 : active_space_env%as_mix_weight(ib) = mixing_store%wmax
375 : END IF
376 0 : active_space_env%as_mix_weight(ib) = MAX(1.0_dp, active_space_env%as_mix_weight(ib))
377 : END IF
378 :
379 0 : ALLOCATE (a(nb, nb))
380 0 : ALLOCATE (b(nb, nb))
381 0 : ALLOCATE (c(nb))
382 0 : ALLOCATE (g(nb))
383 :
384 0 : a(:, :) = 0.0_dp
385 0 : c(:) = 0.0_dp
386 0 : DO j = 1, nb
387 0 : DO k = j, nb
388 : a(k, j) = DOT_PRODUCT(active_space_env%as_mix_r_buffer(j, :), &
389 0 : active_space_env%as_mix_r_buffer(k, :))
390 0 : a(j, k) = a(k, j)
391 : END DO
392 : END DO
393 :
394 0 : DO j = 1, nb
395 0 : c(j) = DOT_PRODUCT(active_space_env%as_mix_r_buffer(j, :), p_res)
396 0 : IF (modified) THEN
397 0 : c(j) = active_space_env%as_mix_weight(j)*c(j)
398 0 : DO k = 1, nb
399 : a(k, j) = active_space_env%as_mix_weight(k)* &
400 0 : active_space_env%as_mix_weight(j)*a(k, j)
401 : END DO
402 0 : a(j, j) = mixing_store%broy_w0*mixing_store%broy_w0 + a(j, j)
403 : ELSE
404 0 : a(j, j) = mixing_store%broy_w0 + a(j, j)
405 : END IF
406 : END DO
407 :
408 0 : CALL invert_matrix(a, b, inv_err)
409 0 : g(:) = 0.0_dp
410 0 : DO j = 1, nb
411 0 : DO k = 1, nb
412 0 : g(j) = g(j) + b(k, j)*c(k)
413 : END DO
414 : END DO
415 :
416 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
417 0 : DO j = 1, nb
418 0 : weight = 1.0_dp
419 0 : IF (modified) weight = active_space_env%as_mix_weight(j)
420 0 : p_mixed = p_mixed - weight*g(j)*active_space_env%as_mix_x_buffer(j, :)
421 : END DO
422 :
423 0 : DEALLOCATE (a)
424 0 : DEALLOCATE (b)
425 0 : DEALLOCATE (c)
426 0 : DEALLOCATE (g)
427 : ELSE
428 0 : active_space_env%as_mix_r_buffer(ib, :) = 0.0_dp
429 0 : active_space_env%as_mix_x_buffer(ib, :) = 0.0_dp
430 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
431 : END IF
432 :
433 0 : active_space_env%as_mix_x_old = p_old
434 0 : active_space_env%as_mix_r_old = p_res
435 0 : IF (modified) THEN
436 0 : mixing_store%iter_method = "MBroy"
437 : ELSE
438 0 : mixing_store%iter_method = "Broy."
439 : END IF
440 :
441 0 : DEALLOCATE (p_res)
442 :
443 0 : END SUBROUTINE active_space_broyden_mixing
444 :
445 : ! **************************************************************************************************
446 : !> \brief Mix the dense active-space density vector.
447 : !> \param active_space_env active space environment
448 : !> \param p_old current active-space density vector
449 : !> \param p_solver active-space density vector returned by the external solver
450 : !> \param p_mixed mixed active-space density vector
451 : ! **************************************************************************************************
452 8 : SUBROUTINE mix_active_density_vector(active_space_env, p_old, p_solver, p_mixed)
453 : TYPE(active_space_type), POINTER :: active_space_env
454 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_old, p_solver
455 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: p_mixed
456 :
457 : INTEGER :: active_mix_iter
458 : TYPE(mixing_storage_type), POINTER :: mixing_store
459 :
460 8 : IF (active_space_env%as_mixing_method == no_mixing_nr .OR. &
461 : .NOT. ASSOCIATED(active_space_env%as_mixing_store)) THEN
462 0 : p_mixed = p_solver
463 : RETURN
464 : END IF
465 :
466 8 : mixing_store => active_space_env%as_mixing_store
467 8 : active_space_env%as_mixing_iter = active_space_env%as_mixing_iter + 1
468 :
469 8 : IF (active_space_env%as_mixing_iter <= mixing_store%nskip_mixing) THEN
470 0 : p_mixed = p_solver
471 0 : mixing_store%iter_method = "NoMix"
472 0 : RETURN
473 : END IF
474 :
475 8 : active_mix_iter = active_space_env%as_mixing_iter - mixing_store%nskip_mixing
476 8 : IF (active_space_env%as_mixing_method == direct_mixing_nr .OR. &
477 : active_mix_iter <= mixing_store%n_simple_mix) THEN
478 8 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
479 8 : mixing_store%iter_method = "P_Mix"
480 8 : RETURN
481 : END IF
482 :
483 0 : SELECT CASE (active_space_env%as_mixing_method)
484 : CASE (pulay_mixing_nr)
485 0 : CALL active_space_pulay_mixing(active_space_env, p_old, p_solver, p_mixed)
486 : CASE (broyden_mixing_nr)
487 0 : CALL active_space_broyden_mixing(active_space_env, p_old, p_solver, p_mixed, .FALSE.)
488 : CASE (modified_broyden_mixing_nr)
489 0 : CALL active_space_broyden_mixing(active_space_env, p_old, p_solver, p_mixed, .TRUE.)
490 : CASE DEFAULT
491 0 : CPABORT("Unsupported ACTIVE_SPACE%MIXING%METHOD.")
492 : END SELECT
493 :
494 : END SUBROUTINE mix_active_density_vector
495 :
496 : ! **************************************************************************************************
497 : !> \brief Update active space density matrix from Fortran arrays
498 : !> \param p_act_mo_a alpha density matrix in active space MO basis
499 : !> \param active_space_env active space environment
500 : !> \param p_act_mo_b beta density matrix in active space MO basis
501 : !> \author Vladimir Rybkin
502 : ! **************************************************************************************************
503 8 : SUBROUTINE update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
504 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_act_mo_a
505 : TYPE(active_space_type), POINTER :: active_space_env
506 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: p_act_mo_b
507 :
508 : INTEGER :: i1, i2, idx, ispin, m1, m2, nact2, &
509 : nmo_active, nspins
510 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: p_mixed, p_old, p_solver
511 : TYPE(cp_fm_type), POINTER :: p_active
512 :
513 8 : nmo_active = active_space_env%nmo_active
514 8 : nact2 = nmo_active*nmo_active
515 8 : nspins = active_space_env%nspins
516 8 : CPASSERT(SIZE(p_act_mo_a) == nact2)
517 8 : IF (nspins == 2) THEN
518 6 : IF (.NOT. PRESENT(p_act_mo_b)) CPABORT("Missing beta active-space density.")
519 6 : CPASSERT(SIZE(p_act_mo_b) == nact2)
520 : END IF
521 :
522 24 : ALLOCATE (p_mixed(nspins*nact2))
523 16 : ALLOCATE (p_old(nspins*nact2))
524 16 : ALLOCATE (p_solver(nspins*nact2))
525 :
526 40 : p_solver(1:nact2) = p_act_mo_a
527 8 : IF (nspins == 2) THEN
528 30 : p_solver(nact2 + 1:2*nact2) = p_act_mo_b
529 : END IF
530 :
531 : idx = 0
532 22 : DO ispin = 1, nspins
533 14 : p_active => active_space_env%p_active(ispin)
534 50 : DO i1 = 1, nmo_active
535 28 : m1 = active_space_env%active_orbitals(i1, ispin)
536 98 : DO i2 = 1, nmo_active
537 56 : idx = idx + 1
538 56 : m2 = active_space_env%active_orbitals(i2, ispin)
539 84 : CALL cp_fm_get_element(p_active, m1, m2, p_old(idx))
540 : END DO
541 : END DO
542 : END DO
543 :
544 8 : CALL mix_active_density_vector(active_space_env, p_old, p_solver, p_mixed)
545 :
546 8 : idx = 0
547 22 : DO ispin = 1, nspins
548 14 : p_active => active_space_env%p_active(ispin)
549 50 : DO i1 = 1, nmo_active
550 28 : m1 = active_space_env%active_orbitals(i1, ispin)
551 98 : DO i2 = 1, nmo_active
552 56 : idx = idx + 1
553 56 : m2 = active_space_env%active_orbitals(i2, ispin)
554 84 : CALL cp_fm_set_element(p_active, m1, m2, p_mixed(idx))
555 : END DO
556 : END DO
557 : END DO
558 :
559 8 : DEALLOCATE (p_mixed)
560 8 : DEALLOCATE (p_old)
561 8 : DEALLOCATE (p_solver)
562 :
563 8 : END SUBROUTINE update_active_density
564 :
565 : END MODULE qs_active_space_mixing
|
