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