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_gspace_mixing
10 :
11 : USE bibliography, ONLY: Broyden1965,&
12 : Johnson1988,&
13 : Kerker1981,&
14 : cite_reference
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE kinds, ONLY: dp
17 : USE mathconstants, ONLY: z_one,&
18 : z_zero
19 : USE mathlib, ONLY: invert_matrix
20 : USE message_passing, ONLY: mp_para_env_type
21 : USE pw_methods, ONLY: pw_integrate_function,&
22 : pw_transfer
23 : USE pw_types, ONLY: pw_c1d_gs_type,&
24 : pw_r3d_rs_type
25 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
26 : cp_1d_z_p_type,&
27 : gspace_mixing_nr,&
28 : mixing_storage_type,&
29 : modified_broyden_mixing_nr,&
30 : multisecant_mixing_nr,&
31 : new_pulay_mixing_nr,&
32 : pulay_mixing_nr
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_mixing_utils, ONLY: mixing_density
36 : USE qs_rho_atom_types, ONLY: rho_atom_type
37 : USE qs_rho_types, ONLY: qs_rho_get,&
38 : qs_rho_type
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gspace_mixing'
46 :
47 : PUBLIC :: gspace_mixing
48 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief Driver for the g-space mixing, calls the proper routine given the
53 : !>requested method
54 : !> \param qs_env ...
55 : !> \param mixing_method ...
56 : !> \param mixing_store ...
57 : !> \param rho ...
58 : !> \param para_env ...
59 : !> \param iter_count ...
60 : !> \param auxiliary Replicated density coordinates, mixed with the same history and spin channels.
61 : !> \par History
62 : !> 05.2009
63 : !> 02.2015 changed input to make g-space mixing available in linear scaling
64 : !> SCF [Patrick Seewald]
65 : !> \author MI
66 : ! **************************************************************************************************
67 11052 : SUBROUTINE gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count, auxiliary)
68 : TYPE(qs_environment_type), POINTER :: qs_env
69 : INTEGER :: mixing_method
70 : TYPE(mixing_storage_type), POINTER :: mixing_store
71 : TYPE(qs_rho_type), POINTER :: rho
72 : TYPE(mp_para_env_type), POINTER :: para_env
73 : INTEGER :: iter_count
74 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: auxiliary(:, :)
75 :
76 : CHARACTER(len=*), PARAMETER :: routineN = 'gspace_mixing'
77 :
78 : COMPLEX(KIND=dp) :: tmp
79 : INTEGER :: handle, iatom, ig, ispin, natom, ng, &
80 : nspin
81 : LOGICAL :: gapw
82 11052 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
83 11052 : TYPE(cp_1d_z_p_type), ALLOCATABLE :: density(:)
84 : TYPE(dft_control_type), POINTER :: dft_control
85 11052 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
86 11052 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
87 11052 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
88 :
89 11052 : CALL timeset(routineN, handle)
90 :
91 11052 : CALL get_qs_env(qs_env, dft_control=dft_control)
92 11052 : IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
93 : dft_control%qs_control%semi_empirical)) THEN
94 :
95 9564 : CPASSERT(ASSOCIATED(mixing_store))
96 9564 : NULLIFY (rho_atom, rho_g, rho_r, tot_rho_r)
97 9564 : CALL qs_rho_get(rho, rho_g=rho_g, rho_r=rho_r, tot_rho_r=tot_rho_r)
98 :
99 9564 : nspin = SIZE(rho_g, 1)
100 18910 : CALL mixing_density(mixing_store, rho_g, density, auxiliary)
101 9564 : ng = SIZE(density(1)%cc)
102 9564 : CPASSERT(ng == SIZE(mixing_store%rhoin(1)%cc))
103 9564 : gapw = dft_control%qs_control%gapw
104 9564 : natom = 0
105 :
106 9564 : IF (mixing_store%gmix_p .AND. gapw) THEN
107 110 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
108 110 : natom = SIZE(rho_atom)
109 : END IF
110 :
111 9564 : IF (nspin == 2) THEN
112 66939571 : DO ig = 1, ng
113 66938333 : tmp = density(1)%cc(ig)
114 66938333 : density(1)%cc(ig) = tmp + density(2)%cc(ig)
115 66939571 : density(2)%cc(ig) = tmp - density(2)%cc(ig)
116 : END DO
117 1238 : IF (mixing_store%gmix_p .AND. gapw) THEN
118 58 : CALL gapw_cpc_spin_to_charge_mag(rho_atom, mixing_store)
119 : END IF
120 : END IF
121 :
122 9564 : IF (iter_count + 1 <= mixing_store%nskip_mixing) THEN
123 : ! skip mixing
124 8 : DO ispin = 1, nspin
125 27656 : DO ig = 1, ng
126 27652 : mixing_store%rhoin(ispin)%cc(ig) = density(ispin)%cc(ig)
127 : END DO
128 : END DO
129 4 : IF (mixing_store%gmix_p .AND. gapw) THEN
130 0 : DO ispin = 1, nspin
131 0 : DO iatom = 1, natom
132 0 : IF (mixing_store%paw(iatom)) THEN
133 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
134 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
135 : END IF
136 : END DO
137 : END DO
138 : END IF
139 :
140 4 : mixing_store%iter_method = "NoMix"
141 9560 : ELSE IF ((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) THEN
142 6 : CALL cite_reference(Kerker1981)
143 6 : CALL gmix_potential_only(qs_env, mixing_store, density)
144 6 : mixing_store%iter_method = "Kerker"
145 9554 : ELSE IF (mixing_method == gspace_mixing_nr) THEN
146 120 : CALL cite_reference(Kerker1981)
147 120 : CALL gmix_potential_only(qs_env, mixing_store, density)
148 120 : mixing_store%iter_method = "Kerker"
149 : ELSE IF (mixing_method == pulay_mixing_nr) THEN
150 282 : CALL pulay_mixing(qs_env, mixing_store, density, para_env, rho_g(1)%pw_grid%vol)
151 282 : mixing_store%iter_method = "Pulay"
152 : ELSE IF (mixing_method == new_pulay_mixing_nr) THEN
153 406 : CALL pulay_mixing(qs_env, mixing_store, density, para_env, rho_g(1)%pw_grid%vol)
154 406 : mixing_store%iter_method = "NPulay"
155 : ELSE IF (mixing_method == broyden_mixing_nr) THEN
156 8738 : CALL cite_reference(Broyden1965)
157 8738 : CALL broyden_mixing(qs_env, mixing_store, density, para_env)
158 8738 : mixing_store%iter_method = "Broy."
159 : ELSE IF (mixing_method == modified_broyden_mixing_nr) THEN
160 8 : CALL cite_reference(Johnson1988)
161 8 : CALL modified_broyden_mixing(qs_env, mixing_store, density, para_env)
162 8 : mixing_store%iter_method = "MBroy"
163 : ELSE IF (mixing_method == multisecant_mixing_nr) THEN
164 0 : CPASSERT(.NOT. gapw)
165 0 : CPASSERT(rho_g(1)%pw_grid%ngpts < HUGE(ng) - mixing_store%naux)
166 : CALL multisecant_mixing(mixing_store, density, para_env, rho_g(1)%pw_grid%gsq, &
167 0 : INT(rho_g(1)%pw_grid%ngpts) + mixing_store%naux)
168 0 : mixing_store%iter_method = "MSec."
169 : END IF
170 :
171 9564 : IF (nspin == 2) THEN
172 1238 : IF (mixing_store%gmix_p .AND. gapw) THEN
173 58 : CALL gapw_cpc_charge_mag_to_spin(rho_atom, mixing_store)
174 : END IF
175 66939571 : DO ig = 1, ng
176 66938333 : tmp = density(1)%cc(ig)
177 66938333 : density(1)%cc(ig) = 0.5_dp*(tmp + density(2)%cc(ig))
178 66939571 : density(2)%cc(ig) = 0.5_dp*(tmp - density(2)%cc(ig))
179 : END DO
180 : END IF
181 :
182 9564 : IF (mixing_store%naux > 0) THEN
183 218 : IF (ALLOCATED(mixing_store%aux_buffer)) THEN
184 109 : ng = SIZE(rho_g(1)%array)
185 247 : DO ispin = 1, nspin
186 1360128 : rho_g(ispin)%array = density(ispin)%cc(:ng)
187 1075 : auxiliary(:, ispin) = REAL(density(ispin)%cc(ng + 1:), dp)*rho_g(ispin)%pw_grid%vol
188 : END DO
189 : END IF
190 4082 : CALL para_env%bcast(auxiliary, 0)
191 : END IF
192 : ! A skipped update already has a consistent real-space density.
193 9564 : IF (iter_count + 1 <= mixing_store%nskip_mixing) THEN
194 4 : CALL timestop(handle)
195 4 : RETURN
196 : END IF
197 20358 : DO ispin = 1, nspin
198 10798 : CALL pw_transfer(rho_g(ispin), rho_r(ispin))
199 20358 : tot_rho_r(ispin) = pw_integrate_function(rho_r(ispin), isign=-1)
200 : END DO
201 : END IF
202 :
203 11048 : CALL timestop(handle)
204 :
205 22108 : END SUBROUTINE gspace_mixing
206 :
207 : ! **************************************************************************************************
208 : !> \brief Convert GAPW compensation charge coefficients from spin-up/spin-down
209 : !> representation to total-charge/magnetization representation.
210 : !> \param rho_atom ...
211 : !> \param mixing_store ...
212 : ! **************************************************************************************************
213 58 : SUBROUTINE gapw_cpc_spin_to_charge_mag(rho_atom, mixing_store)
214 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
215 : TYPE(mixing_storage_type), POINTER :: mixing_store
216 :
217 : INTEGER :: iatom, icoef1, icoef2, ncoef_h1, &
218 : ncoef_h2, ncoef_s1, ncoef_s2
219 : REAL(dp) :: cpc_h_tmp, cpc_s_tmp
220 :
221 58 : CPASSERT(ASSOCIATED(rho_atom))
222 58 : CPASSERT(ASSOCIATED(mixing_store%paw))
223 :
224 522 : DO iatom = 1, SIZE(rho_atom)
225 522 : IF (mixing_store%paw(iatom)) THEN
226 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_h(1)%r_coef))
227 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_h(2)%r_coef))
228 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_s(1)%r_coef))
229 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_s(2)%r_coef))
230 464 : ncoef_h1 = SIZE(rho_atom(iatom)%cpc_h(1)%r_coef, 1)
231 464 : ncoef_h2 = SIZE(rho_atom(iatom)%cpc_h(1)%r_coef, 2)
232 464 : ncoef_s1 = SIZE(rho_atom(iatom)%cpc_s(1)%r_coef, 1)
233 464 : ncoef_s2 = SIZE(rho_atom(iatom)%cpc_s(1)%r_coef, 2)
234 464 : CPASSERT(ncoef_h1 == SIZE(rho_atom(iatom)%cpc_h(2)%r_coef, 1))
235 464 : CPASSERT(ncoef_h2 == SIZE(rho_atom(iatom)%cpc_h(2)%r_coef, 2))
236 464 : CPASSERT(ncoef_s1 == SIZE(rho_atom(iatom)%cpc_s(2)%r_coef, 1))
237 464 : CPASSERT(ncoef_s2 == SIZE(rho_atom(iatom)%cpc_s(2)%r_coef, 2))
238 :
239 12992 : DO icoef2 = 1, ncoef_h2
240 351248 : DO icoef1 = 1, ncoef_h1
241 338256 : cpc_h_tmp = rho_atom(iatom)%cpc_h(1)%r_coef(icoef1, icoef2)
242 : rho_atom(iatom)%cpc_h(1)%r_coef(icoef1, icoef2) = &
243 338256 : cpc_h_tmp + rho_atom(iatom)%cpc_h(2)%r_coef(icoef1, icoef2)
244 : rho_atom(iatom)%cpc_h(2)%r_coef(icoef1, icoef2) = &
245 350784 : cpc_h_tmp - rho_atom(iatom)%cpc_h(2)%r_coef(icoef1, icoef2)
246 : END DO
247 : END DO
248 :
249 12992 : DO icoef2 = 1, ncoef_s2
250 351248 : DO icoef1 = 1, ncoef_s1
251 338256 : cpc_s_tmp = rho_atom(iatom)%cpc_s(1)%r_coef(icoef1, icoef2)
252 : rho_atom(iatom)%cpc_s(1)%r_coef(icoef1, icoef2) = &
253 338256 : cpc_s_tmp + rho_atom(iatom)%cpc_s(2)%r_coef(icoef1, icoef2)
254 : rho_atom(iatom)%cpc_s(2)%r_coef(icoef1, icoef2) = &
255 350784 : cpc_s_tmp - rho_atom(iatom)%cpc_s(2)%r_coef(icoef1, icoef2)
256 : END DO
257 : END DO
258 : END IF
259 : END DO
260 :
261 58 : END SUBROUTINE gapw_cpc_spin_to_charge_mag
262 :
263 : ! **************************************************************************************************
264 : !> \brief Convert GAPW compensation charge coefficients from total-charge/magnetization
265 : !> representation back to spin-up/spin-down representation.
266 : !> \param rho_atom ...
267 : !> \param mixing_store ...
268 : ! **************************************************************************************************
269 58 : SUBROUTINE gapw_cpc_charge_mag_to_spin(rho_atom, mixing_store)
270 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
271 : TYPE(mixing_storage_type), POINTER :: mixing_store
272 :
273 : INTEGER :: iatom, icoef1, icoef2, ncoef_h1, &
274 : ncoef_h2, ncoef_s1, ncoef_s2
275 : REAL(dp) :: cpc_h_tmp, cpc_s_tmp
276 :
277 58 : CPASSERT(ASSOCIATED(rho_atom))
278 58 : CPASSERT(ASSOCIATED(mixing_store%paw))
279 :
280 522 : DO iatom = 1, SIZE(rho_atom)
281 522 : IF (mixing_store%paw(iatom)) THEN
282 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_h(1)%r_coef))
283 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_h(2)%r_coef))
284 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_s(1)%r_coef))
285 464 : CPASSERT(ASSOCIATED(rho_atom(iatom)%cpc_s(2)%r_coef))
286 464 : ncoef_h1 = SIZE(rho_atom(iatom)%cpc_h(1)%r_coef, 1)
287 464 : ncoef_h2 = SIZE(rho_atom(iatom)%cpc_h(1)%r_coef, 2)
288 464 : ncoef_s1 = SIZE(rho_atom(iatom)%cpc_s(1)%r_coef, 1)
289 464 : ncoef_s2 = SIZE(rho_atom(iatom)%cpc_s(1)%r_coef, 2)
290 464 : CPASSERT(ncoef_h1 == SIZE(rho_atom(iatom)%cpc_h(2)%r_coef, 1))
291 464 : CPASSERT(ncoef_h2 == SIZE(rho_atom(iatom)%cpc_h(2)%r_coef, 2))
292 464 : CPASSERT(ncoef_s1 == SIZE(rho_atom(iatom)%cpc_s(2)%r_coef, 1))
293 464 : CPASSERT(ncoef_s2 == SIZE(rho_atom(iatom)%cpc_s(2)%r_coef, 2))
294 :
295 12992 : DO icoef2 = 1, ncoef_h2
296 351248 : DO icoef1 = 1, ncoef_h1
297 338256 : cpc_h_tmp = rho_atom(iatom)%cpc_h(1)%r_coef(icoef1, icoef2)
298 : rho_atom(iatom)%cpc_h(1)%r_coef(icoef1, icoef2) = &
299 338256 : 0.5_dp*(cpc_h_tmp + rho_atom(iatom)%cpc_h(2)%r_coef(icoef1, icoef2))
300 : rho_atom(iatom)%cpc_h(2)%r_coef(icoef1, icoef2) = &
301 350784 : 0.5_dp*(cpc_h_tmp - rho_atom(iatom)%cpc_h(2)%r_coef(icoef1, icoef2))
302 : END DO
303 : END DO
304 :
305 12992 : DO icoef2 = 1, ncoef_s2
306 351248 : DO icoef1 = 1, ncoef_s1
307 338256 : cpc_s_tmp = rho_atom(iatom)%cpc_s(1)%r_coef(icoef1, icoef2)
308 : rho_atom(iatom)%cpc_s(1)%r_coef(icoef1, icoef2) = &
309 338256 : 0.5_dp*(cpc_s_tmp + rho_atom(iatom)%cpc_s(2)%r_coef(icoef1, icoef2))
310 : rho_atom(iatom)%cpc_s(2)%r_coef(icoef1, icoef2) = &
311 350784 : 0.5_dp*(cpc_s_tmp - rho_atom(iatom)%cpc_s(2)%r_coef(icoef1, icoef2))
312 : END DO
313 : END DO
314 : END IF
315 : END DO
316 :
317 58 : END SUBROUTINE gapw_cpc_charge_mag_to_spin
318 :
319 : ! **************************************************************************************************
320 : !> \brief G-space mixing performed via the Kerker damping on the density on the grid
321 : !> thus affecting only the caluclation of the potential, but not the denisity matrix
322 : !> \param qs_env ...
323 : !> \param mixing_store ...
324 : !> \param density Local mixing vectors, including any auxiliary coordinates.
325 : !> \par History
326 : !> 02.2009
327 : !> \author MI
328 : ! **************************************************************************************************
329 :
330 126 : SUBROUTINE gmix_potential_only(qs_env, mixing_store, density)
331 :
332 : TYPE(qs_environment_type), POINTER :: qs_env
333 : TYPE(mixing_storage_type), POINTER :: mixing_store
334 : TYPE(cp_1d_z_p_type), INTENT(IN) :: density(:)
335 :
336 : CHARACTER(len=*), PARAMETER :: routineN = 'gmix_potential_only'
337 :
338 126 : COMPLEX(dp), DIMENSION(:), POINTER :: cc_new
339 : INTEGER :: handle, iatom, ig, ispin, natom, ng, &
340 : nspin
341 : LOGICAL :: gapw
342 : REAL(dp) :: alpha, alpha_eff, f_mix
343 126 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
344 : TYPE(dft_control_type), POINTER :: dft_control
345 126 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
346 :
347 0 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
348 126 : CPASSERT(ASSOCIATED(mixing_store%kerker_factor))
349 :
350 126 : CALL timeset(routineN, handle)
351 126 : NULLIFY (cc_new, dft_control, rho_atom, kerker_eff)
352 :
353 126 : CALL get_qs_env(qs_env, dft_control=dft_control)
354 :
355 126 : nspin = SIZE(density)
356 126 : ng = SIZE(density(1)%cc)
357 :
358 126 : gapw = dft_control%qs_control%gapw
359 126 : alpha = mixing_store%alpha
360 :
361 252 : DO ispin = 1, nspin
362 : ! Select spin-channel-specific mixing parameters
363 126 : IF (ispin >= 2 .AND. nspin >= 2) THEN
364 0 : alpha_eff = mixing_store%alpha_mag
365 0 : kerker_eff => mixing_store%kerker_factor_mag
366 : ELSE
367 126 : alpha_eff = mixing_store%alpha
368 126 : kerker_eff => mixing_store%kerker_factor
369 : END IF
370 126 : cc_new => density(ispin)%cc
371 1334398 : DO ig = 1, mixing_store%ig_max ! ng
372 1334272 : f_mix = alpha_eff*kerker_eff(ig)
373 1334272 : cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
374 1334398 : mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
375 : END DO
376 252 : DO ig = mixing_store%ig_max + 1, ng
377 0 : f_mix = alpha_eff
378 0 : cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
379 126 : mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
380 : END DO
381 :
382 : END DO
383 :
384 126 : IF (mixing_store%gmix_p .AND. gapw) THEN
385 : CALL get_qs_env(qs_env=qs_env, &
386 8 : rho_atom_set=rho_atom)
387 8 : natom = SIZE(rho_atom)
388 16 : DO ispin = 1, nspin
389 8 : IF (ispin >= 2 .AND. nspin >= 2) THEN
390 0 : alpha_eff = mixing_store%alpha_mag
391 : ELSE
392 8 : alpha_eff = mixing_store%alpha
393 : END IF
394 80 : DO iatom = 1, natom
395 72 : IF (mixing_store%paw(iatom)) THEN
396 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
397 7408 : mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
398 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
399 7408 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
400 7408 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
401 7408 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
402 : END IF
403 : END DO
404 : END DO
405 : END IF
406 :
407 126 : CALL timestop(handle)
408 :
409 126 : END SUBROUTINE gmix_potential_only
410 :
411 : ! **************************************************************************************************
412 : !> \brief Pulay Mixing using as metrics for the residual the Kerer damping factor
413 : !> The mixing is applied directly on the density in reciprocal space,
414 : !> therefore it affects the potentials
415 : !> on the grid but not the density matrix
416 : !> \param qs_env ...
417 : !> \param mixing_store ...
418 : !> \param density Local mixing vectors, including any auxiliary coordinates.
419 : !> \param para_env ...
420 : !> \param vol Cell volume for the Pulay metric.
421 : !> \par History
422 : !> 03.2009
423 : !> \author MI
424 : ! **************************************************************************************************
425 :
426 688 : SUBROUTINE pulay_mixing(qs_env, mixing_store, density, para_env, vol)
427 :
428 : TYPE(qs_environment_type), POINTER :: qs_env
429 : TYPE(mixing_storage_type), POINTER :: mixing_store
430 : TYPE(cp_1d_z_p_type), INTENT(IN) :: density(:)
431 : TYPE(mp_para_env_type), POINTER :: para_env
432 : REAL(KIND=dp), INTENT(IN) :: vol
433 :
434 : CHARACTER(len=*), PARAMETER :: routineN = 'pulay_mixing'
435 :
436 688 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: cc_mix
437 : INTEGER :: handle, i, iatom, ib, ibb, ig, ispin, &
438 : jb, n1, n2, natom, nb, nb1, nbuffer, &
439 : ng, nspin
440 : LOGICAL :: gapw
441 : REAL(dp) :: alpha_eff, alpha_kerker, alpha_pulay, &
442 : beta, f_mix, inv_err, norm, &
443 : norm_c_inv, res_norm
444 688 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, ev
445 688 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b, c, c_inv, cpc_h_mix, cpc_s_mix
446 688 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
447 : TYPE(dft_control_type), POINTER :: dft_control
448 688 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
449 :
450 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
451 688 : CPASSERT(ASSOCIATED(mixing_store%rhoin_buffer))
452 688 : NULLIFY (dft_control, rho_atom, kerker_eff)
453 688 : CALL timeset(routineN, handle)
454 :
455 688 : CALL get_qs_env(qs_env, dft_control=dft_control)
456 688 : nspin = SIZE(density)
457 688 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
458 :
459 688 : alpha_kerker = mixing_store%alpha ! default, overridden per spin in loop
460 688 : beta = mixing_store%pulay_beta
461 688 : alpha_pulay = mixing_store%pulay_alpha
462 688 : nbuffer = mixing_store%nbuffer
463 688 : gapw = dft_control%qs_control%gapw
464 688 : IF (gapw) THEN
465 12 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
466 12 : natom = SIZE(rho_atom)
467 : END IF
468 :
469 2064 : ALLOCATE (cc_mix(ng))
470 :
471 1604 : DO ispin = 1, nspin
472 : ! Select spin-channel-specific mixing parameters
473 916 : IF (ispin >= 2 .AND. nspin >= 2) THEN
474 228 : alpha_eff = mixing_store%alpha_mag
475 228 : kerker_eff => mixing_store%kerker_factor_mag
476 : ELSE
477 688 : alpha_eff = mixing_store%alpha
478 688 : kerker_eff => mixing_store%kerker_factor
479 : END IF
480 916 : alpha_kerker = alpha_eff
481 916 : ib = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
482 :
483 916 : mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) + 1
484 916 : nb = MIN(mixing_store%ncall_p(ispin), nbuffer)
485 916 : ibb = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
486 :
487 59408980 : mixing_store%res_buffer(ib, ispin)%cc(:) = CMPLX(0._dp, 0._dp, KIND=dp)
488 916 : norm = 0.0_dp
489 19897492 : IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
490 916 : res_norm = 0.0_dp
491 59408980 : DO ig = 1, ng
492 59408064 : f_mix = kerker_eff(ig)
493 : mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(density(ispin)%cc(ig) - &
494 59408064 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
495 : res_norm = res_norm + &
496 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)*REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
497 59408980 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))*AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
498 : END DO
499 :
500 916 : CALL para_env%sum(res_norm)
501 916 : res_norm = SQRT(res_norm)
502 :
503 916 : IF (res_norm < 1.E-14_dp) THEN
504 228 : mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) - 1
505 : ELSE
506 688 : nb1 = nb + 1
507 688 : IF (ALLOCATED(b)) DEALLOCATE (b)
508 2752 : ALLOCATE (b(nb1, nb1))
509 688 : b = 0.0_dp
510 688 : IF (ALLOCATED(c)) DEALLOCATE (c)
511 2752 : ALLOCATE (c(nb, nb))
512 688 : c = 0.0_dp
513 688 : IF (ALLOCATED(c_inv)) DEALLOCATE (c_inv)
514 2064 : ALLOCATE (c_inv(nb, nb))
515 688 : c_inv = 0.0_dp
516 688 : IF (ALLOCATED(alpha_c)) DEALLOCATE (alpha_c)
517 2064 : ALLOCATE (alpha_c(nb))
518 688 : alpha_c = 0.0_dp
519 688 : norm_c_inv = 0.0_dp
520 688 : IF (ALLOCATED(ev)) DEALLOCATE (ev)
521 2064 : ALLOCATE (ev(nb1))
522 688 : ev = 0.0_dp
523 :
524 : END IF
525 :
526 4310 : DO jb = 1, nb
527 3394 : mixing_store%pulay_matrix(jb, ib, ispin) = 0.0_dp
528 219015490 : DO ig = 1, ng
529 219012096 : f_mix = mixing_store%special_metric(ig)
530 : mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin) + &
531 : f_mix*(REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp) &
532 : *REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
533 : AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
534 219015490 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig)))
535 : END DO
536 3394 : CALL para_env%sum(mixing_store%pulay_matrix(jb, ib, ispin))
537 3394 : mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)*vol*vol
538 4310 : mixing_store%pulay_matrix(ib, jb, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)
539 : END DO
540 :
541 916 : IF (nb == 1 .OR. res_norm < 1.E-14_dp) THEN
542 19896852 : DO ig = 1, ng
543 19896576 : f_mix = alpha_kerker*kerker_eff(ig)
544 : cc_mix(ig) = density(ispin)%cc(ig) - &
545 19896576 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
546 : density(ispin)%cc(ig) = f_mix*cc_mix(ig) + &
547 19896576 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
548 19896852 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = density(ispin)%cc(ig)
549 : END DO
550 276 : IF (mixing_store%gmix_p) THEN
551 2 : IF (gapw) THEN
552 18 : DO iatom = 1, natom
553 18 : IF (mixing_store%paw(iatom)) THEN
554 : mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
555 1850 : mixing_store%cpc_h_in(iatom, ispin)%r_coef
556 : mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
557 1850 : mixing_store%cpc_s_in(iatom, ispin)%r_coef
558 :
559 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
560 1850 : (1.0_dp - alpha_kerker)*mixing_store%cpc_h_in(iatom, ispin)%r_coef
561 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
562 1850 : (1.0_dp - alpha_kerker)*mixing_store%cpc_s_in(iatom, ispin)%r_coef
563 :
564 926 : mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
565 926 : mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
566 :
567 1850 : mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
568 1850 : mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
569 : END IF
570 : END DO
571 : END IF
572 : END IF
573 : ELSE
574 :
575 20148 : b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb, ispin)
576 20148 : c(1:nb, 1:nb) = b(1:nb, 1:nb)
577 3758 : b(nb1, 1:nb) = -1.0_dp
578 3758 : b(1:nb, nb1) = -1.0_dp
579 640 : b(nb1, nb1) = 0.0_dp
580 :
581 640 : CALL invert_matrix(c, c_inv, inv_err, improve=.TRUE.)
582 640 : alpha_c = 0.0_dp
583 3758 : DO i = 1, nb
584 20148 : DO jb = 1, nb
585 16390 : alpha_c(i) = alpha_c(i) + c_inv(jb, i)
586 19508 : norm_c_inv = norm_c_inv + c_inv(jb, i)
587 : END DO
588 : END DO
589 3758 : alpha_c(1:nb) = alpha_c(1:nb)/norm_c_inv
590 640 : cc_mix = CMPLX(0._dp, 0._dp, KIND=dp)
591 3758 : DO jb = 1, nb
592 199119278 : DO ig = 1, ng
593 : cc_mix(ig) = cc_mix(ig) + &
594 : alpha_c(jb)*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) + &
595 199118638 : mixing_store%pulay_beta*mixing_store%res_buffer(jb, ispin)%cc(ig))
596 : END DO
597 : END DO
598 39512128 : mixing_store%rhoin_buffer(ibb, ispin)%cc = CMPLX(0._dp, 0._dp, KIND=dp)
599 640 : IF (alpha_pulay > 0.0_dp) THEN
600 233290 : DO ig = 1, ng
601 233280 : f_mix = alpha_pulay*kerker_eff(ig)
602 : density(ispin)%cc(ig) = f_mix*density(ispin)%cc(ig) + &
603 233280 : (1.0_dp - f_mix)*cc_mix(ig)
604 233290 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = density(ispin)%cc(ig)
605 : END DO
606 : ELSE
607 39278838 : DO ig = 1, ng
608 39278208 : density(ispin)%cc(ig) = cc_mix(ig)
609 39278838 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = density(ispin)%cc(ig)
610 : END DO
611 : END IF
612 :
613 640 : IF (mixing_store%gmix_p .AND. gapw) THEN
614 10 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
615 90 : DO iatom = 1, natom
616 90 : IF (mixing_store%paw(iatom)) THEN
617 10 : n1 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 1)
618 10 : n2 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 2)
619 40 : ALLOCATE (cpc_h_mix(n1, n2))
620 30 : ALLOCATE (cpc_s_mix(n1, n2))
621 :
622 : mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
623 9250 : mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef
624 : mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
625 9250 : mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef
626 10 : cpc_h_mix = 0.0_dp
627 10 : cpc_s_mix = 0.0_dp
628 48 : DO jb = 1, nb
629 : cpc_h_mix(:, :) = cpc_h_mix(:, :) + &
630 : alpha_c(jb)*mixing_store%cpc_h_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
631 17594 : alpha_c(jb)*beta*mixing_store%cpc_h_res_buffer(jb, iatom, ispin)%r_coef(:, :)
632 : cpc_s_mix(:, :) = cpc_s_mix(:, :) + &
633 : alpha_c(jb)*mixing_store%cpc_s_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
634 17604 : alpha_c(jb)*beta*mixing_store%cpc_s_res_buffer(jb, iatom, ispin)%r_coef(:, :)
635 : END DO
636 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
637 4630 : (1.0_dp - alpha_pulay)*cpc_h_mix
638 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
639 4630 : (1.0_dp - alpha_pulay)*cpc_s_mix
640 9250 : mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
641 9250 : mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
642 10 : DEALLOCATE (cpc_h_mix)
643 10 : DEALLOCATE (cpc_s_mix)
644 : END IF
645 : END DO
646 : END IF
647 : END IF ! nb==1
648 :
649 1604 : IF (res_norm > 1.E-14_dp) THEN
650 688 : DEALLOCATE (b)
651 688 : DEALLOCATE (ev)
652 688 : DEALLOCATE (c, c_inv, alpha_c)
653 : END IF
654 :
655 : END DO ! ispin
656 :
657 688 : DEALLOCATE (cc_mix)
658 :
659 688 : CALL timestop(handle)
660 :
661 2064 : END SUBROUTINE pulay_mixing
662 :
663 : ! **************************************************************************************************
664 : !> \brief Broyden Mixing using as metrics for the residual the Kerer damping factor
665 : !> The mixing is applied directly on the density expansion in reciprocal space
666 : !> \param qs_env ...
667 : !> \param mixing_store ...
668 : !> \param density Local mixing vectors, including any auxiliary coordinates.
669 : !> \param para_env ...
670 : !> \par History
671 : !> 03.2009
672 : !> \author MI
673 : ! **************************************************************************************************
674 :
675 8738 : SUBROUTINE broyden_mixing(qs_env, mixing_store, density, para_env)
676 :
677 : TYPE(qs_environment_type), POINTER :: qs_env
678 : TYPE(mixing_storage_type), POINTER :: mixing_store
679 : TYPE(cp_1d_z_p_type), INTENT(IN) :: density(:)
680 : TYPE(mp_para_env_type), POINTER :: para_env
681 :
682 : CHARACTER(len=*), PARAMETER :: routineN = 'broyden_mixing'
683 :
684 : COMPLEX(dp) :: cc_mix
685 8738 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: res_rho
686 : INTEGER :: handle, i, iatom, ib, ig, ispin, j, jb, &
687 : n1, n2, natom, nb, nbuffer, ng, nspin
688 : LOGICAL :: gapw, only_kerker
689 : REAL(dp) :: alpha, alpha_eff, dcpc_h_res, &
690 : dcpc_s_res, delta_norm, f_mix, &
691 : inv_err, res_norm, valh, vals, w0
692 8738 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g
693 8738 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
694 : REAL(dp), DIMENSION(2) :: norms
695 8738 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
696 : TYPE(dft_control_type), POINTER :: dft_control
697 8738 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
698 :
699 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
700 8738 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
701 8738 : CPASSERT(ASSOCIATED(mixing_store%rhoin_old))
702 8738 : CPASSERT(ASSOCIATED(mixing_store%drho_buffer))
703 8738 : NULLIFY (dft_control, rho_atom, kerker_eff)
704 8738 : CALL timeset(routineN, handle)
705 :
706 8738 : CALL get_qs_env(qs_env, dft_control=dft_control)
707 :
708 8738 : nspin = SIZE(density)
709 8738 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
710 :
711 8738 : alpha = mixing_store%alpha
712 8738 : w0 = mixing_store%broy_w0
713 8738 : nbuffer = mixing_store%nbuffer
714 8738 : gapw = dft_control%qs_control%gapw
715 :
716 26214 : ALLOCATE (res_rho(ng))
717 :
718 8738 : mixing_store%ncall = mixing_store%ncall + 1
719 8738 : IF (mixing_store%ncall == 1) THEN
720 : only_kerker = .TRUE.
721 : ELSE
722 7880 : only_kerker = .FALSE.
723 7880 : nb = MIN(mixing_store%ncall - 1, nbuffer)
724 7880 : ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1
725 : END IF
726 :
727 8738 : IF (gapw) THEN
728 1864 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
729 1864 : natom = SIZE(rho_atom)
730 : ELSE
731 : natom = 0
732 : END IF
733 :
734 18486 : DO ispin = 1, nspin
735 : ! Select spin-channel-specific mixing parameters
736 9748 : IF (ispin >= 2 .AND. nspin >= 2) THEN
737 1010 : alpha_eff = mixing_store%alpha_mag
738 1010 : kerker_eff => mixing_store%kerker_factor_mag
739 : ELSE
740 8738 : alpha_eff = mixing_store%alpha
741 8738 : kerker_eff => mixing_store%kerker_factor
742 : END IF
743 18486 : IF (only_kerker) THEN
744 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig, f_mix) &
745 1040 : !$OMP SHARED(ng, ispin, mixing_store, density, res_rho, alpha_eff, kerker_eff)
746 : DO ig = 1, ng
747 : res_rho(ig) = density(ispin)%cc(ig) - mixing_store%rhoin(ispin)%cc(ig)
748 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
749 : f_mix = alpha_eff*kerker_eff(ig)
750 : density(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
751 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
752 : mixing_store%rhoin(ispin)%cc(ig) = density(ispin)%cc(ig)
753 : END DO
754 : !$OMP END PARALLEL DO
755 :
756 1040 : IF (mixing_store%gmix_p) THEN
757 24 : IF (gapw) THEN
758 216 : DO iatom = 1, natom
759 216 : IF (mixing_store%paw(iatom)) THEN
760 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
761 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef
762 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
763 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef
764 :
765 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
766 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
767 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
768 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
769 :
770 122972 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
771 122972 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
772 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
773 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
774 : END IF
775 : END DO
776 : END IF
777 : END IF
778 : ELSE
779 :
780 34832 : ALLOCATE (a(nb, nb))
781 8708 : a = 0.0_dp
782 26124 : ALLOCATE (b(nb, nb))
783 8708 : b = 0.0_dp
784 26124 : ALLOCATE (c(nb))
785 8708 : c = 0.0_dp
786 17416 : ALLOCATE (g(nb))
787 8708 : g = 0.0_dp
788 :
789 8708 : delta_norm = 0.0_dp
790 8708 : res_norm = 0.0_dp
791 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig) &
792 8708 : !$OMP SHARED(ng, ib, ispin, mixing_store, density, res_rho) REDUCTION(+:delta_norm, res_norm)
793 : DO ig = 1, ng
794 : res_rho(ig) = density(ispin)%cc(ig) - mixing_store%rhoin(ispin)%cc(ig)
795 : mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
796 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
797 : res_norm = res_norm + &
798 : REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + &
799 : AIMAG(res_rho(ig))*AIMAG(res_rho(ig))
800 : delta_norm = delta_norm + &
801 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)* &
802 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
803 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))* &
804 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
805 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
806 : mixing_store%rhoin(ispin)%cc(ig) - &
807 : mixing_store%rhoin_old(ispin)%cc(ig)
808 : END DO
809 : !$OMP END PARALLEL DO
810 26124 : norms(:) = [delta_norm, res_norm]
811 8708 : CALL para_env%sum(norms)
812 8708 : delta_norm = SQRT(norms(1))
813 8708 : res_norm = SQRT(norms(2))
814 :
815 8708 : IF (res_norm > 1.E-14_dp) THEN
816 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig, f_mix) &
817 8256 : !$OMP SHARED(ng, ib, ispin, mixing_store, delta_norm, alpha_eff, kerker_eff)
818 : DO ig = 1, ng
819 : mixing_store%res_buffer(ib, ispin)%cc(ig) = mixing_store%res_buffer(ib, ispin)%cc(ig)/delta_norm
820 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = mixing_store%drho_buffer(ib, ispin)%cc(ig)/delta_norm
821 : f_mix = alpha_eff*kerker_eff(ig)
822 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
823 : f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
824 : mixing_store%drho_buffer(ib, ispin)%cc(ig)
825 : END DO
826 : !$OMP END PARALLEL DO
827 :
828 8256 : IF (mixing_store%gmix_p .AND. gapw) THEN
829 900 : DO iatom = 1, natom
830 900 : IF (mixing_store%paw(iatom)) THEN
831 604 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
832 604 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
833 16744 : DO i = 1, n2
834 448996 : DO j = 1, n1
835 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
836 : (mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
837 432252 : mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
838 : dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
839 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
840 432252 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
841 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
842 432252 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
843 :
844 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
845 : (mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
846 432252 : mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
847 : dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
848 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
849 432252 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
850 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
851 432252 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
852 :
853 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
854 : alpha_eff*dcpc_h_res + &
855 432252 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
856 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
857 : alpha_eff*dcpc_s_res + &
858 448392 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
859 : END DO
860 : END DO
861 : END IF
862 : END DO
863 : END IF
864 :
865 : CALL broyden_history_products(mixing_store%res_buffer(1:nb, ispin), res_rho, &
866 8256 : mixing_store%p_metric, mixing_store%ig_max, a, c)
867 8256 : CALL para_env%sum(a)
868 8256 : CALL para_env%sum(c)
869 44150 : DO jb = 1, nb
870 44150 : a(jb, jb) = w0 + a(jb, jb)
871 : END DO
872 8256 : CALL invert_matrix(a, b, inv_err)
873 :
874 8256 : CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1)
875 : ELSE
876 452 : G = 0.0_dp
877 : END IF
878 :
879 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig, jb, cc_mix, f_mix) &
880 8708 : !$OMP SHARED(ng, nb, ispin, mixing_store, density, res_rho, res_norm, g, alpha_eff, kerker_eff)
881 : DO ig = 1, ng
882 : cc_mix = z_zero
883 : DO jb = 1, nb
884 : cc_mix = cc_mix - G(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
885 : END DO
886 : f_mix = alpha_eff*kerker_eff(ig)
887 :
888 : IF (res_norm > 1.E-14_dp) density(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
889 : f_mix*res_rho(ig) + cc_mix
890 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
891 : mixing_store%rhoin(ispin)%cc(ig) = density(ispin)%cc(ig)
892 : END DO
893 : !$OMP END PARALLEL DO
894 :
895 8708 : IF (mixing_store%gmix_p) THEN
896 124 : IF (gapw) THEN
897 1116 : DO iatom = 1, natom
898 1116 : IF (mixing_store%paw(iatom)) THEN
899 796 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
900 796 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
901 22120 : DO i = 1, n2
902 594340 : DO j = 1, n1
903 572220 : valh = 0.0_dp
904 572220 : vals = 0.0_dp
905 2884428 : DO jb = 1, nb
906 2312208 : valh = valh - G(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
907 2884428 : vals = vals - G(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
908 : END DO
909 593544 : IF (res_norm > 1.E-14_dp) THEN
910 : rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
911 : alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
912 432252 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha_eff) + valh
913 : rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
914 : alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
915 432252 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha_eff) + vals
916 : END IF
917 : END DO
918 : END DO
919 :
920 594340 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
921 594340 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
922 1187884 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
923 1187884 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
924 : END IF
925 : END DO
926 : END IF
927 : END IF
928 :
929 8708 : DEALLOCATE (a, b, c, g)
930 : END IF
931 :
932 : END DO ! ispin
933 :
934 8738 : DEALLOCATE (res_rho)
935 :
936 8738 : CALL timestop(handle)
937 :
938 26214 : END SUBROUTINE broyden_mixing
939 :
940 : ! **************************************************************************************************
941 : !> \brief Modified Broyden mixing with dynamic residual weights.
942 : !> The mixing is applied directly on the density expansion in reciprocal space.
943 : !> \param qs_env ...
944 : !> \param mixing_store ...
945 : !> \param density Local mixing vectors, including any auxiliary coordinates.
946 : !> \param para_env ...
947 : ! **************************************************************************************************
948 :
949 8 : SUBROUTINE modified_broyden_mixing(qs_env, mixing_store, density, para_env)
950 :
951 : TYPE(qs_environment_type), POINTER :: qs_env
952 : TYPE(mixing_storage_type), POINTER :: mixing_store
953 : TYPE(cp_1d_z_p_type), INTENT(IN) :: density(:)
954 : TYPE(mp_para_env_type), POINTER :: para_env
955 :
956 : CHARACTER(len=*), PARAMETER :: routineN = 'modified_broyden_mixing'
957 :
958 : COMPLEX(dp) :: cc_mix
959 8 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: res_rho
960 : INTEGER :: handle, i, iatom, ib, ig, ispin, j, jb, &
961 : kb, n1, n2, natom, nb, nbuffer, ng, &
962 : nspin
963 : LOGICAL :: can_update, gapw, only_kerker
964 : REAL(dp) :: alpha_eff, broy_w0, broy_wmax, broy_wmin, broy_wref, dcpc_h_res, dcpc_s_res, &
965 : delta_norm, f_mix, inv_err, res_norm, valh, vals
966 8 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g
967 8 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
968 : REAL(dp), DIMENSION(2) :: norms
969 8 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
970 : TYPE(dft_control_type), POINTER :: dft_control
971 8 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
972 :
973 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
974 8 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
975 8 : CPASSERT(ASSOCIATED(mixing_store%rhoin_old))
976 8 : CPASSERT(ASSOCIATED(mixing_store%drho_buffer))
977 8 : CPASSERT(ASSOCIATED(mixing_store%weight))
978 8 : NULLIFY (dft_control, rho_atom, kerker_eff)
979 8 : CALL timeset(routineN, handle)
980 :
981 8 : CALL get_qs_env(qs_env, dft_control=dft_control)
982 :
983 8 : nspin = SIZE(density)
984 8 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
985 :
986 8 : broy_w0 = mixing_store%broy_w0
987 8 : broy_wref = mixing_store%wc
988 8 : broy_wmax = mixing_store%wmax
989 8 : broy_wmin = 1.0_dp
990 8 : nbuffer = mixing_store%nbuffer
991 8 : gapw = dft_control%qs_control%gapw
992 :
993 24 : ALLOCATE (res_rho(ng))
994 :
995 8 : mixing_store%ncall = mixing_store%ncall + 1
996 8 : IF (mixing_store%ncall == 1) THEN
997 : only_kerker = .TRUE.
998 : ELSE
999 6 : only_kerker = .FALSE.
1000 6 : nb = MIN(mixing_store%ncall - 1, nbuffer)
1001 6 : ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1
1002 : END IF
1003 :
1004 8 : IF (gapw) THEN
1005 0 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
1006 0 : natom = SIZE(rho_atom)
1007 : ELSE
1008 : natom = 0
1009 : END IF
1010 :
1011 16 : DO ispin = 1, nspin
1012 8 : IF (ispin >= 2 .AND. nspin >= 2) THEN
1013 0 : alpha_eff = mixing_store%alpha_mag
1014 0 : kerker_eff => mixing_store%kerker_factor_mag
1015 : ELSE
1016 8 : alpha_eff = mixing_store%alpha
1017 8 : kerker_eff => mixing_store%kerker_factor
1018 : END IF
1019 16 : IF (only_kerker) THEN
1020 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig, f_mix) &
1021 2 : !$OMP SHARED(ng, ispin, mixing_store, density, res_rho, alpha_eff, kerker_eff)
1022 : DO ig = 1, ng
1023 : res_rho(ig) = density(ispin)%cc(ig) - mixing_store%rhoin(ispin)%cc(ig)
1024 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
1025 : f_mix = alpha_eff*kerker_eff(ig)
1026 : density(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
1027 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
1028 : mixing_store%rhoin(ispin)%cc(ig) = density(ispin)%cc(ig)
1029 : END DO
1030 : !$OMP END PARALLEL DO
1031 :
1032 2 : IF (mixing_store%gmix_p) THEN
1033 0 : IF (gapw) THEN
1034 0 : DO iatom = 1, natom
1035 0 : IF (mixing_store%paw(iatom)) THEN
1036 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
1037 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef
1038 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
1039 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef
1040 :
1041 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
1042 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
1043 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
1044 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
1045 :
1046 0 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
1047 0 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
1048 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
1049 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
1050 : END IF
1051 : END DO
1052 : END IF
1053 : END IF
1054 : ELSE
1055 :
1056 24 : ALLOCATE (a(nb, nb))
1057 6 : a = 0.0_dp
1058 18 : ALLOCATE (b(nb, nb))
1059 6 : b = 0.0_dp
1060 18 : ALLOCATE (c(nb))
1061 6 : c = 0.0_dp
1062 12 : ALLOCATE (g(nb))
1063 6 : g = 0.0_dp
1064 :
1065 6 : delta_norm = 0.0_dp
1066 6 : res_norm = 0.0_dp
1067 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig) &
1068 6 : !$OMP SHARED(ng, ib, ispin, mixing_store, density, res_rho) REDUCTION(+:delta_norm, res_norm)
1069 : DO ig = 1, ng
1070 : res_rho(ig) = density(ispin)%cc(ig) - mixing_store%rhoin(ispin)%cc(ig)
1071 : mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
1072 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
1073 : res_norm = res_norm + &
1074 : REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + &
1075 : AIMAG(res_rho(ig))*AIMAG(res_rho(ig))
1076 : delta_norm = delta_norm + &
1077 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)* &
1078 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
1079 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))* &
1080 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
1081 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
1082 : mixing_store%rhoin(ispin)%cc(ig) - &
1083 : mixing_store%rhoin_old(ispin)%cc(ig)
1084 : END DO
1085 : !$OMP END PARALLEL DO
1086 18 : norms(:) = [delta_norm, res_norm]
1087 6 : CALL para_env%sum(norms)
1088 6 : delta_norm = SQRT(norms(1))
1089 6 : res_norm = SQRT(norms(2))
1090 :
1091 6 : IF (res_norm > (broy_wref/broy_wmax)) THEN
1092 6 : mixing_store%weight(ib, ispin) = broy_wref/res_norm
1093 : ELSE
1094 0 : mixing_store%weight(ib, ispin) = broy_wmax
1095 : END IF
1096 6 : mixing_store%weight(ib, ispin) = MAX(broy_wmin, mixing_store%weight(ib, ispin))
1097 6 : can_update = res_norm > 1.E-14_dp .AND. delta_norm > 1.E-14_dp
1098 :
1099 6 : IF (can_update) THEN
1100 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig, f_mix) &
1101 6 : !$OMP SHARED(ng, ib, ispin, mixing_store, delta_norm, alpha_eff, kerker_eff)
1102 : DO ig = 1, ng
1103 : mixing_store%res_buffer(ib, ispin)%cc(ig) = mixing_store%res_buffer(ib, ispin)%cc(ig)/delta_norm
1104 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = mixing_store%drho_buffer(ib, ispin)%cc(ig)/delta_norm
1105 : f_mix = alpha_eff*kerker_eff(ig)
1106 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
1107 : f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
1108 : mixing_store%drho_buffer(ib, ispin)%cc(ig)
1109 : END DO
1110 : !$OMP END PARALLEL DO
1111 :
1112 6 : IF (mixing_store%gmix_p .AND. gapw) THEN
1113 0 : DO iatom = 1, natom
1114 0 : IF (mixing_store%paw(iatom)) THEN
1115 0 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
1116 0 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
1117 0 : DO i = 1, n2
1118 0 : DO j = 1, n1
1119 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
1120 : (mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
1121 0 : mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
1122 : dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
1123 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
1124 0 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
1125 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
1126 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
1127 :
1128 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
1129 : (mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
1130 0 : mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
1131 : dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
1132 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
1133 0 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
1134 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
1135 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
1136 :
1137 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
1138 : alpha_eff*dcpc_h_res + &
1139 0 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
1140 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
1141 : alpha_eff*dcpc_s_res + &
1142 0 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
1143 : END DO
1144 : END DO
1145 : END IF
1146 : END DO
1147 : END IF
1148 :
1149 : CALL broyden_history_products(mixing_store%res_buffer(1:nb, ispin), res_rho, &
1150 6 : mixing_store%p_metric, mixing_store%ig_max, a, c)
1151 6 : CALL para_env%sum(a)
1152 6 : CALL para_env%sum(c)
1153 :
1154 18 : DO jb = 1, nb
1155 12 : c(jb) = mixing_store%weight(jb, ispin)*c(jb)
1156 40 : DO kb = 1, nb
1157 40 : a(kb, jb) = mixing_store%weight(kb, ispin)*mixing_store%weight(jb, ispin)*a(kb, jb)
1158 : END DO
1159 18 : a(jb, jb) = broy_w0*broy_w0 + a(jb, jb)
1160 : END DO
1161 6 : CALL invert_matrix(a, b, inv_err)
1162 :
1163 6 : CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1)
1164 : ELSE
1165 0 : g = 0.0_dp
1166 0 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig) SHARED(ng, ib, ispin, mixing_store)
1167 : DO ig = 1, ng
1168 : mixing_store%res_buffer(ib, ispin)%cc(ig) = z_zero
1169 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = z_zero
1170 : END DO
1171 : !$OMP END PARALLEL DO
1172 : END IF
1173 :
1174 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) PRIVATE(ig, jb, cc_mix, f_mix) &
1175 6 : !$OMP SHARED(ng, nb, ispin, mixing_store, density, res_rho, res_norm, g, alpha_eff, kerker_eff)
1176 : DO ig = 1, ng
1177 : cc_mix = z_zero
1178 : DO jb = 1, nb
1179 : cc_mix = cc_mix - mixing_store%weight(jb, ispin)*G(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
1180 : END DO
1181 : f_mix = alpha_eff*kerker_eff(ig)
1182 :
1183 : IF (res_norm > 1.E-14_dp) density(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
1184 : f_mix*res_rho(ig) + cc_mix
1185 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
1186 : mixing_store%rhoin(ispin)%cc(ig) = density(ispin)%cc(ig)
1187 : END DO
1188 : !$OMP END PARALLEL DO
1189 :
1190 6 : IF (mixing_store%gmix_p) THEN
1191 0 : IF (gapw) THEN
1192 0 : DO iatom = 1, natom
1193 0 : IF (mixing_store%paw(iatom)) THEN
1194 0 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
1195 0 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
1196 0 : DO i = 1, n2
1197 0 : DO j = 1, n1
1198 0 : valh = 0.0_dp
1199 0 : vals = 0.0_dp
1200 0 : DO jb = 1, nb
1201 : valh = valh - mixing_store%weight(jb, ispin)*G(jb)* &
1202 0 : mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
1203 : vals = vals - mixing_store%weight(jb, ispin)*G(jb)* &
1204 0 : mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
1205 : END DO
1206 0 : IF (res_norm > 1.E-14_dp) THEN
1207 : rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
1208 : alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
1209 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha_eff) + valh
1210 : rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
1211 : alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
1212 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha_eff) + vals
1213 : END IF
1214 : END DO
1215 : END DO
1216 :
1217 0 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
1218 0 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
1219 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
1220 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
1221 : END IF
1222 : END DO
1223 : END IF
1224 : END IF
1225 :
1226 6 : DEALLOCATE (a, b, c, g)
1227 : END IF
1228 :
1229 : END DO
1230 :
1231 8 : DEALLOCATE (res_rho)
1232 :
1233 8 : CALL timestop(handle)
1234 :
1235 24 : END SUBROUTINE modified_broyden_mixing
1236 :
1237 : ! **************************************************************************************************
1238 : !> \brief Local metric products for Broyden history vectors and the current residual.
1239 : !> \param history Normalized residual differences.
1240 : !> \param residual Current residual.
1241 : !> \param metric Mixing metric.
1242 : !> \param ng Number of local entries included in the metric.
1243 : !> \param a Symmetric history overlap matrix.
1244 : !> \param c History projections of the current residual.
1245 : ! **************************************************************************************************
1246 8262 : SUBROUTINE broyden_history_products(history, residual, metric, ng, a, c)
1247 : TYPE(cp_1d_z_p_type), INTENT(IN) :: history(:)
1248 : COMPLEX(dp), INTENT(IN) :: residual(:)
1249 : REAL(dp), INTENT(IN) :: metric(:)
1250 : INTEGER, INTENT(IN) :: ng
1251 : REAL(dp), INTENT(OUT) :: a(:, :), c(:)
1252 :
1253 : INTEGER :: ig, jb, kb, nb
1254 : REAL(dp) :: prod
1255 :
1256 8262 : nb = SIZE(history)
1257 247322 : a(:, :) = 0.0_dp
1258 44168 : c(:) = 0.0_dp
1259 : ! Partition the long vectors; each thread accumulates private small arrays.
1260 : !$OMP PARALLEL DEFAULT(NONE) SHARED(history, residual, metric, ng, nb) &
1261 8262 : !$OMP PRIVATE(ig, jb, kb, prod) REDUCTION(+:A, c)
1262 : DO jb = 1, nb
1263 : DO kb = jb, nb
1264 : prod = 0.0_dp
1265 : !$OMP DO SCHEDULE(STATIC)
1266 : DO ig = 1, ng
1267 : prod = prod + metric(ig)*( &
1268 : REAL(history(jb)%cc(ig), dp)*REAL(history(kb)%cc(ig), dp) + &
1269 : AIMAG(history(jb)%cc(ig))*AIMAG(history(kb)%cc(ig)))
1270 : END DO
1271 : !$OMP END DO NOWAIT
1272 : a(kb, jb) = prod
1273 : END DO
1274 : prod = 0.0_dp
1275 : !$OMP DO SCHEDULE(STATIC)
1276 : DO ig = 1, ng
1277 : prod = prod + metric(ig)*( &
1278 : REAL(history(jb)%cc(ig), dp)*REAL(residual(ig), dp) + &
1279 : AIMAG(history(jb)%cc(ig))*AIMAG(residual(ig)))
1280 : END DO
1281 : !$OMP END DO NOWAIT
1282 : c(jb) = prod
1283 : END DO
1284 : !$OMP END PARALLEL
1285 44168 : DO jb = 1, nb
1286 127792 : DO kb = jb + 1, nb
1287 119530 : a(jb, kb) = a(kb, jb)
1288 : END DO
1289 : END DO
1290 8262 : END SUBROUTINE broyden_history_products
1291 :
1292 : ! **************************************************************************************************
1293 : !> \brief Multisecant scheme to perform charge density Mixing
1294 : !> as preconditioner we use the Kerker damping factor
1295 : !> The mixing is applied directly on the density in reciprocal space,
1296 : !> therefore it affects the potentials
1297 : !> on the grid but not the density matrix
1298 : !> \param mixing_store ...
1299 : !> \param density Local mixing vectors, including any auxiliary coordinates.
1300 : !> \param para_env ...
1301 : !> \param g2 Squared physical G vectors; auxiliary coordinates are not Fourier modes.
1302 : !> \param ng_global Total distributed vector length including auxiliary coordinates.
1303 : !> \par History
1304 : !> 05.2009
1305 : !> \author MI
1306 : ! **************************************************************************************************
1307 0 : SUBROUTINE multisecant_mixing(mixing_store, density, para_env, g2, ng_global)
1308 :
1309 : TYPE(mixing_storage_type), POINTER :: mixing_store
1310 : TYPE(cp_1d_z_p_type), INTENT(IN) :: density(:)
1311 : TYPE(mp_para_env_type), POINTER :: para_env
1312 : REAL(KIND=dp), INTENT(IN) :: g2(:)
1313 : INTEGER, INTENT(IN) :: ng_global
1314 :
1315 : CHARACTER(len=*), PARAMETER :: routineN = 'multisecant_mixing'
1316 : COMPLEX(KIND=dp), PARAMETER :: cmone = (-1.0_dp, 0.0_dp), &
1317 : cone = (1.0_dp, 0.0_dp), &
1318 : czero = (0.0_dp, 0.0_dp)
1319 :
1320 : COMPLEX(dp) :: saa, yaa
1321 0 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: gn_global, pgn, res_matrix_global, &
1322 0 : tmp_vec, ugn
1323 0 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, res_matrix, sa, step_matrix, ya
1324 0 : COMPLEX(dp), DIMENSION(:), POINTER :: gn
1325 : INTEGER :: handle, ib, ib_next, ib_prev, ig, &
1326 : ig_global, iig, ispin, jb, kb, nb, &
1327 : nbuffer, ng, nspin
1328 : LOGICAL :: use_zgemm, use_zgemm_rev
1329 : REAL(dp) :: alpha, alpha_eff, f_mix, gn_norm, gn_norm_old, inv_err, n_low, n_up, pgn_norm, &
1330 : prec, r_step, reg_par, sigma_max, sigma_tilde, step_size
1331 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: norm_res, norm_res_low, norm_res_up
1332 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b_matrix, binv_matrix
1333 0 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
1334 : REAL(dp), SAVE :: sigma_old = 1.0_dp
1335 :
1336 0 : CPASSERT(ASSOCIATED(mixing_store))
1337 :
1338 0 : CALL timeset(routineN, handle)
1339 :
1340 0 : NULLIFY (gn)
1341 :
1342 0 : use_zgemm = .FALSE.
1343 0 : use_zgemm_rev = .TRUE.
1344 :
1345 : ! prepare the parameters
1346 :
1347 0 : nspin = SIZE(density)
1348 0 : ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
1349 0 : alpha = mixing_store%alpha
1350 :
1351 0 : sigma_max = mixing_store%sigma_max
1352 0 : reg_par = mixing_store%reg_par
1353 0 : r_step = mixing_store%r_step
1354 0 : nbuffer = mixing_store%nbuffer
1355 :
1356 : ! determine the step number, and multisecant iteration
1357 0 : nb = MIN(mixing_store%ncall, nbuffer - 1)
1358 0 : ib = MODULO(mixing_store%ncall, nbuffer) + 1
1359 0 : IF (mixing_store%ncall > 0) THEN
1360 0 : ib_prev = MODULO(mixing_store%ncall - 1, nbuffer) + 1
1361 : ELSE
1362 : ib_prev = 0
1363 : END IF
1364 0 : mixing_store%ncall = mixing_store%ncall + 1
1365 0 : ib_next = MODULO(mixing_store%ncall, nbuffer) + 1
1366 :
1367 : ! compute the residual gn and its norm gn_norm
1368 0 : DO ispin = 1, nspin
1369 : ! Start from the last accepted state, including skipped/simple mixing steps.
1370 0 : IF (nb == 0) mixing_store%rhoin_buffer(ib, ispin)%cc = mixing_store%rhoin(ispin)%cc
1371 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1372 0 : gn_norm = 0.0_dp
1373 0 : DO ig = 1, ng
1374 0 : gn(ig) = (density(ispin)%cc(ig) - mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1375 : gn_norm = gn_norm + &
1376 0 : REAL(gn(ig), dp)*REAL(gn(ig), dp) + AIMAG(gn(ig))*AIMAG(gn(ig))
1377 : END DO
1378 0 : CALL para_env%sum(gn_norm)
1379 0 : gn_norm = SQRT(gn_norm)
1380 0 : mixing_store%norm_res_buffer(ib, ispin) = gn_norm
1381 : END DO
1382 :
1383 0 : IF (nb == 0) THEN
1384 : !simple mixing
1385 0 : DO ispin = 1, nspin
1386 0 : IF (ispin >= 2 .AND. nspin >= 2) THEN
1387 0 : alpha_eff = mixing_store%alpha_mag
1388 0 : kerker_eff => mixing_store%kerker_factor_mag
1389 : ELSE
1390 0 : alpha_eff = mixing_store%alpha
1391 0 : kerker_eff => mixing_store%kerker_factor
1392 : END IF
1393 0 : DO ig = 1, ng
1394 0 : f_mix = alpha_eff*kerker_eff(ig)
1395 : density(ispin)%cc(ig) = mixing_store%rhoin_buffer(1, ispin)%cc(ig) + &
1396 0 : f_mix*mixing_store%res_buffer(1, ispin)%cc(ig)
1397 0 : mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = density(ispin)%cc(ig)
1398 : END DO
1399 : END DO
1400 0 : CALL timestop(handle)
1401 0 : RETURN
1402 : END IF
1403 :
1404 : ! allocate temporary arrays
1405 : ! step_matrix S ngxnb
1406 0 : ALLOCATE (step_matrix(ng, nb))
1407 : ! res_matrix Y ngxnb
1408 0 : ALLOCATE (res_matrix(ng, nb))
1409 : ! matrix A nbxnb
1410 0 : ALLOCATE (a_matrix(nb, ng_global))
1411 : ! PSI nb vector of norms
1412 0 : ALLOCATE (norm_res(nb))
1413 0 : ALLOCATE (norm_res_low(nb))
1414 0 : ALLOCATE (norm_res_up(nb))
1415 :
1416 : ! matrix B nbxnb
1417 0 : ALLOCATE (b_matrix(nb, nb))
1418 : ! matrix B_inv nbxnb
1419 0 : ALLOCATE (binv_matrix(nb, nb))
1420 :
1421 0 : ALLOCATE (gn_global(ng_global))
1422 0 : ALLOCATE (res_matrix_global(ng_global))
1423 : IF (use_zgemm) THEN
1424 : ALLOCATE (sa(ng, ng_global))
1425 : ALLOCATE (ya(ng, ng_global))
1426 : END IF
1427 : IF (use_zgemm_rev) THEN
1428 0 : ALLOCATE (tmp_vec(nb))
1429 : END IF
1430 0 : ALLOCATE (pgn(ng))
1431 0 : ALLOCATE (ugn(ng))
1432 :
1433 0 : DO ispin = 1, nspin
1434 : ! Select spin-channel-specific mixing parameters
1435 0 : IF (ispin >= 2 .AND. nspin >= 2) THEN
1436 0 : alpha_eff = mixing_store%alpha_mag
1437 0 : kerker_eff => mixing_store%kerker_factor_mag
1438 : ELSE
1439 0 : alpha_eff = mixing_store%alpha
1440 0 : kerker_eff => mixing_store%kerker_factor
1441 : END IF
1442 : ! generate the global vector with the present residual
1443 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1444 0 : gn_global = z_zero
1445 0 : DO ig = 1, ng
1446 0 : ig_global = mixing_store%ig_global_index(ig)
1447 0 : gn_global(ig_global) = gn(ig)
1448 : END DO
1449 0 : CALL para_env%sum(gn_global)
1450 :
1451 : ! compute steps (matrix S) and residual differences (matrix Y)
1452 : ! with respect to the present
1453 : ! step and the present residual (use stored rho_in and res_buffer)
1454 :
1455 : ! These quantities are pre-conditioned by means of Kerker factor multipliccation
1456 :
1457 0 : DO jb = 1, nb + 1
1458 0 : IF (jb < ib) THEN
1459 : kb = jb
1460 0 : ELSE IF (jb > ib) THEN
1461 0 : kb = jb - 1
1462 : ELSE
1463 : CYCLE
1464 : END IF
1465 0 : norm_res(kb) = 0.0_dp
1466 0 : norm_res_low(kb) = 0.0_dp
1467 0 : norm_res_up(kb) = 0.0_dp
1468 0 : DO ig = 1, ng
1469 :
1470 0 : prec = 1.0_dp
1471 :
1472 : step_matrix(ig, kb) = prec*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) - &
1473 0 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1474 : res_matrix(ig, kb) = (mixing_store%res_buffer(jb, ispin)%cc(ig) - &
1475 0 : mixing_store%res_buffer(ib, ispin)%cc(ig))
1476 : norm_res(kb) = norm_res(kb) + REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1477 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1478 0 : IF (ig <= SIZE(g2)) THEN
1479 0 : IF (g2(ig) < 4.0_dp) THEN
1480 : norm_res_low(kb) = norm_res_low(kb) + &
1481 : REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1482 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1483 : ELSE
1484 : norm_res_up(kb) = norm_res_up(kb) + &
1485 : REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1486 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1487 : END IF
1488 : END IF
1489 0 : res_matrix(ig, kb) = prec*res_matrix(ig, kb)
1490 : END DO
1491 : END DO !jb
1492 :
1493 : ! normalize each column of S and Y => Snorm Ynorm
1494 0 : CALL para_env%sum(norm_res)
1495 0 : CALL para_env%sum(norm_res_up)
1496 0 : CALL para_env%sum(norm_res_low)
1497 0 : norm_res(1:nb) = 1.0_dp/SQRT(norm_res(1:nb))
1498 : n_low = 0.0_dp
1499 : n_up = 0.0_dp
1500 0 : DO jb = 1, nb
1501 0 : n_low = n_low + norm_res_low(jb)/norm_res(jb)
1502 0 : n_up = n_up + norm_res_up(jb)/norm_res(jb)
1503 : END DO
1504 : ! Auxiliary-only changes have no Fourier low/high norm to rebalance.
1505 0 : IF (n_low > 0.0_dp .AND. n_up > 0.0_dp) THEN
1506 0 : DO ig = 1, SIZE(g2)
1507 0 : IF (g2(ig) > 4.0_dp) THEN
1508 0 : step_matrix(ig, 1:nb) = step_matrix(ig, 1:nb)*SQRT(n_low/n_up)
1509 0 : res_matrix(ig, 1:nb) = res_matrix(ig, 1:nb)*SQRT(n_low/n_up)
1510 : END IF
1511 : END DO
1512 : END IF
1513 0 : DO kb = 1, nb
1514 0 : step_matrix(1:ng, kb) = step_matrix(1:ng, kb)*norm_res(kb)
1515 0 : res_matrix(1:ng, kb) = res_matrix(1:ng, kb)*norm_res(kb)
1516 : END DO
1517 :
1518 : ! compute A as [(Ynorm^t Ynorm) + (alpha I)]^(-1) Ynorm^t
1519 : ! compute B
1520 0 : DO jb = 1, nb
1521 0 : DO kb = 1, nb
1522 0 : b_matrix(kb, jb) = 0.0_dp
1523 0 : DO ig = 1, ng
1524 : ! it is assumed that summing over all G vector gives a real, because
1525 : ! y(-G,kb) = (y(G,kb))*
1526 0 : b_matrix(kb, jb) = b_matrix(kb, jb) + REAL(res_matrix(ig, kb)*res_matrix(ig, jb), dp)
1527 : END DO
1528 : END DO
1529 : END DO
1530 :
1531 0 : CALL para_env%sum(b_matrix)
1532 0 : DO jb = 1, nb
1533 0 : b_matrix(jb, jb) = b_matrix(jb, jb) + reg_par
1534 : END DO
1535 : ! invert B
1536 0 : CALL invert_matrix(b_matrix, binv_matrix, inv_err)
1537 :
1538 : ! A = Binv Ynorm^t
1539 0 : a_matrix = z_zero
1540 0 : DO kb = 1, nb
1541 0 : res_matrix_global = z_zero
1542 0 : DO ig = 1, ng
1543 0 : ig_global = mixing_store%ig_global_index(ig)
1544 0 : res_matrix_global(ig_global) = res_matrix(ig, kb)
1545 : END DO
1546 0 : CALL para_env%sum(res_matrix_global)
1547 :
1548 0 : DO jb = 1, nb
1549 0 : DO ig = 1, ng
1550 0 : ig_global = mixing_store%ig_global_index(ig)
1551 : a_matrix(jb, ig_global) = a_matrix(jb, ig_global) + &
1552 0 : binv_matrix(jb, kb)*res_matrix_global(ig_global)
1553 : END DO
1554 : END DO
1555 : END DO
1556 0 : CALL para_env%sum(a_matrix)
1557 :
1558 : ! compute the two components of gn that will be used to update rho
1559 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1560 0 : pgn_norm = 0.0_dp
1561 :
1562 : IF (use_zgemm) THEN
1563 :
1564 : CALL zgemm("N", "N", ng, ng_global, nb, cmone, step_matrix(1, 1), ng, &
1565 : a_matrix(1, 1), nb, czero, sa(1, 1), ng)
1566 : CALL zgemm("N", "N", ng, ng_global, nb, cmone, res_matrix(1, 1), ng, &
1567 : a_matrix(1, 1), nb, czero, ya(1, 1), ng)
1568 : DO ig = 1, ng
1569 : ig_global = mixing_store%ig_global_index(ig)
1570 : ya(ig, ig_global) = ya(ig, ig_global) + z_one
1571 : END DO
1572 :
1573 : CALL zgemv("N", ng, ng_global, cone, sa(1, 1), &
1574 : ng, gn_global(1), 1, czero, pgn(1), 1)
1575 : CALL zgemv("N", ng, ng_global, cone, ya(1, 1), &
1576 : ng, gn_global(1), 1, czero, ugn(1), 1)
1577 :
1578 : DO ig = 1, ng
1579 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1580 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1581 : END DO
1582 : CALL para_env%sum(pgn_norm)
1583 : ELSE IF (use_zgemm_rev) THEN
1584 :
1585 : CALL zgemv("N", nb, ng_global, cone, a_matrix(1, 1), &
1586 0 : nb, gn_global(1), 1, czero, tmp_vec(1), 1)
1587 :
1588 : CALL zgemv("N", ng, nb, cmone, step_matrix(1, 1), ng, &
1589 0 : tmp_vec(1), 1, czero, pgn(1), 1)
1590 :
1591 : CALL zgemv("N", ng, nb, cmone, res_matrix(1, 1), ng, &
1592 0 : tmp_vec(1), 1, czero, ugn(1), 1)
1593 :
1594 0 : DO ig = 1, ng
1595 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1596 0 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1597 0 : ugn(ig) = ugn(ig) + gn(ig)
1598 : END DO
1599 0 : CALL para_env%sum(pgn_norm)
1600 :
1601 : ELSE
1602 : DO ig = 1, ng
1603 : pgn(ig) = z_zero
1604 : ugn(ig) = z_zero
1605 : ig_global = mixing_store%ig_global_index(ig)
1606 : DO iig = 1, ng_global
1607 : saa = z_zero
1608 : yaa = z_zero
1609 :
1610 : IF (ig_global == iig) yaa = z_one
1611 :
1612 : DO jb = 1, nb
1613 : saa = saa - step_matrix(ig, jb)*a_matrix(jb, iig)
1614 : yaa = yaa - res_matrix(ig, jb)*a_matrix(jb, iig)
1615 : END DO
1616 : pgn(ig) = pgn(ig) + saa*gn_global(iig)
1617 : ugn(ig) = ugn(ig) + yaa*gn_global(iig)
1618 : END DO
1619 : END DO
1620 : DO ig = 1, ng
1621 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1622 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1623 : END DO
1624 : CALL para_env%sum(pgn_norm)
1625 : END IF
1626 :
1627 0 : gn_norm = mixing_store%norm_res_buffer(ib, ispin)
1628 0 : gn_norm_old = mixing_store%norm_res_buffer(ib_prev, ispin)
1629 : IF (ib_prev /= 0) THEN
1630 0 : sigma_tilde = sigma_old*MAX(0.5_dp, MIN(2.0_dp, gn_norm_old/gn_norm))
1631 : ELSE
1632 : sigma_tilde = 0.5_dp
1633 : END IF
1634 0 : sigma_tilde = 0.1_dp
1635 : ! Step size for the unpredicted component
1636 0 : step_size = MIN(sigma_tilde, r_step*pgn_norm/gn_norm, sigma_max)
1637 0 : sigma_old = step_size
1638 :
1639 : ! update the density
1640 0 : DO ig = 1, ng
1641 0 : prec = kerker_eff(ig)
1642 : density(ispin)%cc(ig) = mixing_store%rhoin_buffer(ib, ispin)%cc(ig) &
1643 0 : - prec*step_size*ugn(ig) + prec*pgn(ig) ! - 0.1_dp * prec* gn(ig)
1644 0 : mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = density(ispin)%cc(ig)
1645 : END DO
1646 :
1647 : END DO ! ispin
1648 :
1649 : ! Deallocate temporary arrays
1650 0 : DEALLOCATE (step_matrix, res_matrix)
1651 0 : DEALLOCATE (norm_res)
1652 0 : DEALLOCATE (norm_res_up)
1653 0 : DEALLOCATE (norm_res_low)
1654 0 : DEALLOCATE (a_matrix, b_matrix, binv_matrix)
1655 0 : DEALLOCATE (ugn, pgn)
1656 : IF (use_zgemm) THEN
1657 : DEALLOCATE (sa, ya)
1658 : END IF
1659 : IF (use_zgemm_rev) THEN
1660 0 : DEALLOCATE (tmp_vec)
1661 : END IF
1662 0 : DEALLOCATE (gn_global, res_matrix_global)
1663 :
1664 0 : CALL timestop(handle)
1665 :
1666 0 : END SUBROUTINE multisecant_mixing
1667 :
1668 : END MODULE qs_gspace_mixing
|
