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 0 : res_norm = NORM2(active_space_env%as_mix_r_buffer(ib, :))
260 :
261 0 : IF (nb == 1 .OR. res_norm < 1.E-14_dp) THEN
262 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
263 : ELSE
264 0 : ALLOCATE (c(nb, nb))
265 0 : ALLOCATE (c_inv(nb, nb))
266 0 : ALLOCATE (alpha_c(nb))
267 0 : ALLOCATE (p_diis(ndim))
268 :
269 0 : c(:, :) = 0.0_dp
270 0 : DO i = 1, nb
271 0 : DO j = i, nb
272 : c(j, i) = DOT_PRODUCT(active_space_env%as_mix_r_buffer(i, :), &
273 0 : active_space_env%as_mix_r_buffer(j, :))
274 0 : c(i, j) = c(j, i)
275 : END DO
276 : END DO
277 :
278 0 : CALL invert_matrix(c, c_inv, inv_err, improve=.TRUE.)
279 0 : norm_c_inv = SUM(c_inv)
280 0 : IF (ABS(norm_c_inv) < 1.E-14_dp) THEN
281 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
282 : ELSE
283 0 : DO i = 1, nb
284 0 : alpha_c(i) = SUM(c_inv(:, i))/norm_c_inv
285 : END DO
286 :
287 0 : p_diis(:) = 0.0_dp
288 0 : DO i = 1, nb
289 : p_diis(:) = p_diis(:) + alpha_c(i)*(active_space_env%as_mix_x_buffer(i, :) + &
290 0 : mixing_store%pulay_beta*active_space_env%as_mix_r_buffer(i, :))
291 : END DO
292 0 : IF (mixing_store%pulay_alpha > 0.0_dp) THEN
293 : p_mixed = mixing_store%pulay_alpha*p_solver + &
294 0 : (1.0_dp - mixing_store%pulay_alpha)*p_diis
295 : ELSE
296 0 : p_mixed = p_diis
297 : END IF
298 : END IF
299 :
300 0 : DEALLOCATE (alpha_c)
301 0 : DEALLOCATE (p_diis)
302 0 : DEALLOCATE (c)
303 0 : DEALLOCATE (c_inv)
304 : END IF
305 :
306 0 : active_space_env%as_mix_x_buffer(ibb, :) = p_mixed
307 0 : mixing_store%iter_method = "Pulay"
308 :
309 0 : END SUBROUTINE active_space_pulay_mixing
310 :
311 : ! **************************************************************************************************
312 : !> \brief Apply original or modified Broyden mixing to the dense active-space density vector.
313 : !> \param active_space_env active space environment
314 : !> \param p_old current active-space density vector
315 : !> \param p_solver active-space density vector returned by the external solver
316 : !> \param p_mixed mixed active-space density vector
317 : !> \param modified use dynamic residual weights of modified Broyden
318 : ! **************************************************************************************************
319 0 : SUBROUTINE active_space_broyden_mixing(active_space_env, p_old, p_solver, p_mixed, modified)
320 : TYPE(active_space_type), POINTER :: active_space_env
321 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_old, p_solver
322 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: p_mixed
323 : LOGICAL, INTENT(IN) :: modified
324 :
325 : INTEGER :: ib, j, k, nb, nbuffer, ndim
326 : LOGICAL :: can_update
327 : REAL(KIND=dp) :: delta_norm, inv_err, res_norm, weight
328 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c, g, p_res
329 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
330 : TYPE(mixing_storage_type), POINTER :: mixing_store
331 :
332 0 : ndim = SIZE(p_old)
333 0 : CALL ensure_active_mixing_history(active_space_env, ndim)
334 0 : mixing_store => active_space_env%as_mixing_store
335 0 : nbuffer = mixing_store%nbuffer
336 :
337 0 : ALLOCATE (p_res(ndim))
338 0 : p_res(:) = p_solver(:) - p_old(:)
339 0 : res_norm = NORM2(p_res)
340 :
341 0 : mixing_store%ncall = mixing_store%ncall + 1
342 0 : IF (mixing_store%ncall == 1) THEN
343 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
344 0 : active_space_env%as_mix_x_old = p_old
345 0 : active_space_env%as_mix_r_old = p_res
346 0 : IF (modified) THEN
347 0 : mixing_store%iter_method = "MBroy"
348 : ELSE
349 0 : mixing_store%iter_method = "Broy."
350 : END IF
351 0 : DEALLOCATE (p_res)
352 0 : RETURN
353 : END IF
354 :
355 0 : nb = MIN(mixing_store%ncall - 1, nbuffer)
356 0 : ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1
357 :
358 0 : active_space_env%as_mix_r_buffer(ib, :) = p_res - active_space_env%as_mix_r_old
359 0 : active_space_env%as_mix_x_buffer(ib, :) = p_old - active_space_env%as_mix_x_old
360 0 : delta_norm = NORM2(active_space_env%as_mix_r_buffer(ib, :))
361 0 : can_update = res_norm > 1.E-14_dp .AND. delta_norm > 1.E-14_dp
362 :
363 : IF (can_update) THEN
364 0 : active_space_env%as_mix_r_buffer(ib, :) = active_space_env%as_mix_r_buffer(ib, :)/delta_norm
365 : active_space_env%as_mix_x_buffer(ib, :) = active_space_env%as_mix_x_buffer(ib, :)/delta_norm + &
366 0 : mixing_store%alpha*active_space_env%as_mix_r_buffer(ib, :)
367 :
368 0 : IF (modified) THEN
369 0 : IF (res_norm > (mixing_store%wc/mixing_store%wmax)) THEN
370 0 : active_space_env%as_mix_weight(ib) = mixing_store%wc/res_norm
371 : ELSE
372 0 : active_space_env%as_mix_weight(ib) = mixing_store%wmax
373 : END IF
374 0 : active_space_env%as_mix_weight(ib) = MAX(1.0_dp, active_space_env%as_mix_weight(ib))
375 : END IF
376 :
377 0 : ALLOCATE (a(nb, nb))
378 0 : ALLOCATE (b(nb, nb))
379 0 : ALLOCATE (c(nb))
380 0 : ALLOCATE (g(nb))
381 :
382 0 : a(:, :) = 0.0_dp
383 0 : c(:) = 0.0_dp
384 0 : DO j = 1, nb
385 0 : DO k = j, nb
386 : a(k, j) = DOT_PRODUCT(active_space_env%as_mix_r_buffer(j, :), &
387 0 : active_space_env%as_mix_r_buffer(k, :))
388 0 : a(j, k) = a(k, j)
389 : END DO
390 : END DO
391 :
392 0 : DO j = 1, nb
393 0 : c(j) = DOT_PRODUCT(active_space_env%as_mix_r_buffer(j, :), p_res)
394 0 : IF (modified) THEN
395 0 : c(j) = active_space_env%as_mix_weight(j)*c(j)
396 0 : DO k = 1, nb
397 : a(k, j) = active_space_env%as_mix_weight(k)* &
398 0 : active_space_env%as_mix_weight(j)*a(k, j)
399 : END DO
400 0 : a(j, j) = mixing_store%broy_w0*mixing_store%broy_w0 + a(j, j)
401 : ELSE
402 0 : a(j, j) = mixing_store%broy_w0 + a(j, j)
403 : END IF
404 : END DO
405 :
406 0 : CALL invert_matrix(a, b, inv_err)
407 0 : g(:) = 0.0_dp
408 0 : DO j = 1, nb
409 0 : DO k = 1, nb
410 0 : g(j) = g(j) + b(k, j)*c(k)
411 : END DO
412 : END DO
413 :
414 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
415 0 : DO j = 1, nb
416 0 : weight = 1.0_dp
417 0 : IF (modified) weight = active_space_env%as_mix_weight(j)
418 0 : p_mixed = p_mixed - weight*g(j)*active_space_env%as_mix_x_buffer(j, :)
419 : END DO
420 :
421 0 : DEALLOCATE (a)
422 0 : DEALLOCATE (b)
423 0 : DEALLOCATE (c)
424 0 : DEALLOCATE (g)
425 : ELSE
426 0 : active_space_env%as_mix_r_buffer(ib, :) = 0.0_dp
427 0 : active_space_env%as_mix_x_buffer(ib, :) = 0.0_dp
428 0 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
429 : END IF
430 :
431 0 : active_space_env%as_mix_x_old = p_old
432 0 : active_space_env%as_mix_r_old = p_res
433 0 : IF (modified) THEN
434 0 : mixing_store%iter_method = "MBroy"
435 : ELSE
436 0 : mixing_store%iter_method = "Broy."
437 : END IF
438 :
439 0 : DEALLOCATE (p_res)
440 :
441 0 : END SUBROUTINE active_space_broyden_mixing
442 :
443 : ! **************************************************************************************************
444 : !> \brief Mix the dense active-space density vector.
445 : !> \param active_space_env active space environment
446 : !> \param p_old current active-space density vector
447 : !> \param p_solver active-space density vector returned by the external solver
448 : !> \param p_mixed mixed active-space density vector
449 : ! **************************************************************************************************
450 8 : SUBROUTINE mix_active_density_vector(active_space_env, p_old, p_solver, p_mixed)
451 : TYPE(active_space_type), POINTER :: active_space_env
452 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_old, p_solver
453 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: p_mixed
454 :
455 : INTEGER :: active_mix_iter
456 : TYPE(mixing_storage_type), POINTER :: mixing_store
457 :
458 8 : IF (active_space_env%as_mixing_method == no_mixing_nr .OR. &
459 : .NOT. ASSOCIATED(active_space_env%as_mixing_store)) THEN
460 0 : p_mixed = p_solver
461 : RETURN
462 : END IF
463 :
464 8 : mixing_store => active_space_env%as_mixing_store
465 8 : active_space_env%as_mixing_iter = active_space_env%as_mixing_iter + 1
466 :
467 8 : IF (active_space_env%as_mixing_iter <= mixing_store%nskip_mixing) THEN
468 0 : p_mixed = p_solver
469 0 : mixing_store%iter_method = "NoMix"
470 0 : RETURN
471 : END IF
472 :
473 8 : active_mix_iter = active_space_env%as_mixing_iter - mixing_store%nskip_mixing
474 8 : IF (active_space_env%as_mixing_method == direct_mixing_nr .OR. &
475 : active_mix_iter <= mixing_store%n_simple_mix) THEN
476 8 : CALL active_space_direct_mix(p_old, p_solver, mixing_store%alpha, p_mixed)
477 8 : mixing_store%iter_method = "P_Mix"
478 8 : RETURN
479 : END IF
480 :
481 0 : SELECT CASE (active_space_env%as_mixing_method)
482 : CASE (pulay_mixing_nr)
483 0 : CALL active_space_pulay_mixing(active_space_env, p_old, p_solver, p_mixed)
484 : CASE (broyden_mixing_nr)
485 0 : CALL active_space_broyden_mixing(active_space_env, p_old, p_solver, p_mixed, .FALSE.)
486 : CASE (modified_broyden_mixing_nr)
487 0 : CALL active_space_broyden_mixing(active_space_env, p_old, p_solver, p_mixed, .TRUE.)
488 : CASE DEFAULT
489 0 : CPABORT("Unsupported ACTIVE_SPACE%MIXING%METHOD.")
490 : END SELECT
491 :
492 : END SUBROUTINE mix_active_density_vector
493 :
494 : ! **************************************************************************************************
495 : !> \brief Update active space density matrix from Fortran arrays
496 : !> \param p_act_mo_a alpha density matrix in active space MO basis
497 : !> \param active_space_env active space environment
498 : !> \param p_act_mo_b beta density matrix in active space MO basis
499 : !> \author Vladimir Rybkin
500 : ! **************************************************************************************************
501 8 : SUBROUTINE update_active_density(p_act_mo_a, active_space_env, p_act_mo_b)
502 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: p_act_mo_a
503 : TYPE(active_space_type), POINTER :: active_space_env
504 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: p_act_mo_b
505 :
506 : INTEGER :: i1, i2, idx, ispin, m1, m2, nact2, &
507 : nmo_active, nspins
508 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: p_mixed, p_old, p_solver
509 : TYPE(cp_fm_type), POINTER :: p_active
510 :
511 8 : nmo_active = active_space_env%nmo_active
512 8 : nact2 = nmo_active*nmo_active
513 8 : nspins = active_space_env%nspins
514 8 : CPASSERT(SIZE(p_act_mo_a) == nact2)
515 8 : IF (nspins == 2) THEN
516 6 : IF (.NOT. PRESENT(p_act_mo_b)) CPABORT("Missing beta active-space density.")
517 6 : CPASSERT(SIZE(p_act_mo_b) == nact2)
518 : END IF
519 :
520 24 : ALLOCATE (p_mixed(nspins*nact2))
521 16 : ALLOCATE (p_old(nspins*nact2))
522 16 : ALLOCATE (p_solver(nspins*nact2))
523 :
524 40 : p_solver(1:nact2) = p_act_mo_a
525 8 : IF (nspins == 2) THEN
526 30 : p_solver(nact2 + 1:2*nact2) = p_act_mo_b
527 : END IF
528 :
529 : idx = 0
530 22 : DO ispin = 1, nspins
531 14 : p_active => active_space_env%p_active(ispin)
532 50 : DO i1 = 1, nmo_active
533 28 : m1 = active_space_env%active_orbitals(i1, ispin)
534 98 : DO i2 = 1, nmo_active
535 56 : idx = idx + 1
536 56 : m2 = active_space_env%active_orbitals(i2, ispin)
537 84 : CALL cp_fm_get_element(p_active, m1, m2, p_old(idx))
538 : END DO
539 : END DO
540 : END DO
541 :
542 8 : CALL mix_active_density_vector(active_space_env, p_old, p_solver, p_mixed)
543 :
544 8 : idx = 0
545 22 : DO ispin = 1, nspins
546 14 : p_active => active_space_env%p_active(ispin)
547 50 : DO i1 = 1, nmo_active
548 28 : m1 = active_space_env%active_orbitals(i1, ispin)
549 98 : DO i2 = 1, nmo_active
550 56 : idx = idx + 1
551 56 : m2 = active_space_env%active_orbitals(i2, ispin)
552 84 : CALL cp_fm_set_element(p_active, m1, m2, p_mixed(idx))
553 : END DO
554 : END DO
555 : END DO
556 :
557 8 : DEALLOCATE (p_mixed)
558 8 : DEALLOCATE (p_old)
559 8 : DEALLOCATE (p_solver)
560 :
561 8 : END SUBROUTINE update_active_density
562 :
563 : END MODULE qs_active_space_mixing
|
