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 cp_control_types, ONLY: dft_control_type
12 : USE kinds, ONLY: dp
13 : USE mathconstants, ONLY: z_one,&
14 : z_zero
15 : USE mathlib, ONLY: invert_matrix
16 : USE message_passing, ONLY: mp_para_env_type
17 : USE pw_env_types, ONLY: pw_env_get,&
18 : pw_env_type
19 : USE pw_methods, ONLY: pw_axpy,&
20 : pw_copy,&
21 : pw_integrate_function,&
22 : pw_scale,&
23 : pw_transfer,&
24 : pw_zero
25 : USE pw_pool_types, ONLY: pw_pool_type
26 : USE pw_types, ONLY: pw_c1d_gs_type,&
27 : pw_r3d_rs_type
28 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
29 : gspace_mixing_nr,&
30 : mixing_storage_type,&
31 : multisecant_mixing_nr,&
32 : pulay_mixing_nr
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_rho_atom_types, ONLY: rho_atom_type
36 : USE qs_rho_types, ONLY: qs_rho_get,&
37 : qs_rho_type
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gspace_mixing'
45 :
46 : PUBLIC :: gspace_mixing
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Driver for the g-space mixing, calls the proper routine given the
52 : !>requested method
53 : !> \param qs_env ...
54 : !> \param mixing_method ...
55 : !> \param mixing_store ...
56 : !> \param rho ...
57 : !> \param para_env ...
58 : !> \param iter_count ...
59 : !> \par History
60 : !> 05.2009
61 : !> 02.2015 changed input to make g-space mixing available in linear scaling
62 : !> SCF [Patrick Seewald]
63 : !> \author MI
64 : ! **************************************************************************************************
65 2518 : SUBROUTINE gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
66 : TYPE(qs_environment_type), POINTER :: qs_env
67 : INTEGER :: mixing_method
68 : TYPE(mixing_storage_type), POINTER :: mixing_store
69 : TYPE(qs_rho_type), POINTER :: rho
70 : TYPE(mp_para_env_type), POINTER :: para_env
71 : INTEGER :: iter_count
72 :
73 : CHARACTER(len=*), PARAMETER :: routineN = 'gspace_mixing'
74 :
75 : INTEGER :: handle, iatom, ig, ispin, natom, ng, &
76 : nimg, nspin
77 : LOGICAL :: gapw
78 : REAL(dp) :: alpha
79 2518 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
80 : TYPE(dft_control_type), POINTER :: dft_control
81 : TYPE(pw_c1d_gs_type) :: rho_tmp
82 2518 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
83 : TYPE(pw_env_type), POINTER :: pw_env
84 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
85 2518 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
86 2518 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
87 :
88 2518 : CALL timeset(routineN, handle)
89 :
90 2518 : CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env)
91 2518 : IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
92 : dft_control%qs_control%semi_empirical)) THEN
93 :
94 2120 : CPASSERT(ASSOCIATED(mixing_store))
95 2120 : NULLIFY (auxbas_pw_pool, rho_atom, rho_g, rho_r, tot_rho_r)
96 2120 : CALL qs_rho_get(rho, rho_g=rho_g, rho_r=rho_r, tot_rho_r=tot_rho_r)
97 :
98 2120 : nspin = SIZE(rho_g, 1)
99 2120 : nimg = dft_control%nimages
100 2120 : ng = SIZE(rho_g(1)%pw_grid%gsq)
101 2120 : CPASSERT(ng == SIZE(mixing_store%rhoin(1)%cc))
102 :
103 2120 : alpha = mixing_store%alpha
104 2120 : gapw = dft_control%qs_control%gapw
105 :
106 2120 : IF (nspin == 2) THEN
107 302 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
108 302 : CALL auxbas_pw_pool%create_pw(rho_tmp)
109 302 : CALL pw_zero(rho_tmp)
110 302 : CALL pw_copy(rho_g(1), rho_tmp)
111 302 : CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
112 302 : CALL pw_axpy(rho_tmp, rho_g(2), -1.0_dp)
113 302 : CALL pw_scale(rho_g(2), -1.0_dp)
114 : END IF
115 :
116 2120 : IF (iter_count + 1 <= mixing_store%nskip_mixing) THEN
117 : ! skip mixing
118 8 : DO ispin = 1, nspin
119 27656 : DO ig = 1, ng
120 27652 : mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
121 : END DO
122 : END DO
123 4 : IF (mixing_store%gmix_p .AND. gapw) THEN
124 0 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
125 0 : natom = SIZE(rho_atom)
126 0 : DO ispin = 1, nspin
127 0 : DO iatom = 1, natom
128 0 : IF (mixing_store%paw(iatom)) THEN
129 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
130 0 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
131 : END IF
132 : END DO
133 : END DO
134 : END IF
135 :
136 4 : mixing_store%iter_method = "NoMix"
137 4 : IF (nspin == 2) THEN
138 0 : CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
139 0 : CALL pw_scale(rho_g(1), 0.5_dp)
140 0 : CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
141 0 : CALL pw_scale(rho_g(2), -1.0_dp)
142 0 : CALL auxbas_pw_pool%give_back_pw(rho_tmp)
143 : END IF
144 4 : CALL timestop(handle)
145 4 : RETURN
146 : END IF
147 :
148 2116 : IF ((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) THEN
149 6 : CALL gmix_potential_only(qs_env, mixing_store, rho)
150 6 : mixing_store%iter_method = "Kerker"
151 2110 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
152 112 : CALL gmix_potential_only(qs_env, mixing_store, rho)
153 112 : mixing_store%iter_method = "Kerker"
154 1998 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
155 246 : CALL pulay_mixing(qs_env, mixing_store, rho, para_env)
156 246 : mixing_store%iter_method = "Pulay"
157 1752 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
158 1752 : CALL broyden_mixing(qs_env, mixing_store, rho, para_env)
159 1752 : mixing_store%iter_method = "Broy."
160 0 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
161 0 : CPASSERT(.NOT. gapw)
162 0 : CALL multisecant_mixing(mixing_store, rho, para_env)
163 0 : mixing_store%iter_method = "MSec."
164 : END IF
165 :
166 2116 : IF (nspin == 2) THEN
167 302 : CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
168 302 : CALL pw_scale(rho_g(1), 0.5_dp)
169 302 : CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
170 302 : CALL pw_scale(rho_g(2), -1.0_dp)
171 302 : CALL auxbas_pw_pool%give_back_pw(rho_tmp)
172 : END IF
173 :
174 4534 : DO ispin = 1, nspin
175 2418 : CALL pw_transfer(rho_g(ispin), rho_r(ispin))
176 4534 : tot_rho_r(ispin) = pw_integrate_function(rho_r(ispin), isign=-1)
177 : END DO
178 : END IF
179 :
180 2514 : CALL timestop(handle)
181 :
182 2518 : END SUBROUTINE gspace_mixing
183 :
184 : ! **************************************************************************************************
185 : !> \brief G-space mixing performed via the Kerker damping on the density on the grid
186 : !> thus affecting only the caluclation of the potential, but not the denisity matrix
187 : !> \param qs_env ...
188 : !> \param mixing_store ...
189 : !> \param rho ...
190 : !> \par History
191 : !> 02.2009
192 : !> \author MI
193 : ! **************************************************************************************************
194 :
195 236 : SUBROUTINE gmix_potential_only(qs_env, mixing_store, rho)
196 :
197 : TYPE(qs_environment_type), POINTER :: qs_env
198 : TYPE(mixing_storage_type), POINTER :: mixing_store
199 : TYPE(qs_rho_type), POINTER :: rho
200 :
201 : CHARACTER(len=*), PARAMETER :: routineN = 'gmix_potential_only'
202 :
203 118 : COMPLEX(dp), DIMENSION(:), POINTER :: cc_new
204 : INTEGER :: handle, iatom, ig, ispin, natom, ng, &
205 : nspin
206 : LOGICAL :: gapw
207 : REAL(dp) :: alpha, alpha_eff, f_mix
208 118 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
209 : TYPE(dft_control_type), POINTER :: dft_control
210 118 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
211 118 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
212 :
213 0 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
214 118 : CPASSERT(ASSOCIATED(mixing_store%kerker_factor))
215 :
216 118 : CALL timeset(routineN, handle)
217 118 : NULLIFY (cc_new, dft_control, rho_atom, rho_g, kerker_eff)
218 :
219 118 : CALL get_qs_env(qs_env, dft_control=dft_control)
220 118 : CALL qs_rho_get(rho, rho_g=rho_g)
221 :
222 118 : nspin = SIZE(rho_g, 1)
223 118 : ng = SIZE(rho_g(1)%pw_grid%gsq)
224 :
225 118 : gapw = dft_control%qs_control%gapw
226 118 : alpha = mixing_store%alpha
227 :
228 236 : DO ispin = 1, nspin
229 : ! Select spin-channel-specific mixing parameters
230 118 : IF (ispin >= 2 .AND. nspin >= 2) THEN
231 0 : alpha_eff = mixing_store%alpha_mag
232 0 : kerker_eff => mixing_store%kerker_factor_mag
233 : ELSE
234 118 : alpha_eff = mixing_store%alpha
235 118 : kerker_eff => mixing_store%kerker_factor
236 : END IF
237 118 : cc_new => rho_g(ispin)%array
238 1078390 : DO ig = 1, mixing_store%ig_max ! ng
239 1078272 : f_mix = alpha_eff*kerker_eff(ig)
240 1078272 : cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
241 1078390 : mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
242 : END DO
243 236 : DO ig = mixing_store%ig_max + 1, ng
244 0 : f_mix = alpha_eff
245 0 : cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
246 118 : mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
247 : END DO
248 :
249 : END DO
250 :
251 118 : IF (mixing_store%gmix_p .AND. gapw) THEN
252 : CALL get_qs_env(qs_env=qs_env, &
253 8 : rho_atom_set=rho_atom)
254 8 : natom = SIZE(rho_atom)
255 16 : DO ispin = 1, nspin
256 8 : IF (ispin >= 2 .AND. nspin >= 2) THEN
257 0 : alpha_eff = mixing_store%alpha_mag
258 : ELSE
259 8 : alpha_eff = mixing_store%alpha
260 : END IF
261 80 : DO iatom = 1, natom
262 72 : IF (mixing_store%paw(iatom)) THEN
263 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
264 7408 : mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
265 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
266 7408 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
267 7408 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
268 7408 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
269 : END IF
270 : END DO
271 : END DO
272 : END IF
273 :
274 118 : CALL timestop(handle)
275 :
276 118 : END SUBROUTINE gmix_potential_only
277 :
278 : ! **************************************************************************************************
279 : !> \brief Pulay Mixing using as metrics for the residual the Kerer damping factor
280 : !> The mixing is applied directly on the density in reciprocal space,
281 : !> therefore it affects the potentials
282 : !> on the grid but not the density matrix
283 : !> \param qs_env ...
284 : !> \param mixing_store ...
285 : !> \param rho ...
286 : !> \param para_env ...
287 : !> \par History
288 : !> 03.2009
289 : !> \author MI
290 : ! **************************************************************************************************
291 :
292 246 : SUBROUTINE pulay_mixing(qs_env, mixing_store, rho, para_env)
293 :
294 : TYPE(qs_environment_type), POINTER :: qs_env
295 : TYPE(mixing_storage_type), POINTER :: mixing_store
296 : TYPE(qs_rho_type), POINTER :: rho
297 : TYPE(mp_para_env_type), POINTER :: para_env
298 :
299 : CHARACTER(len=*), PARAMETER :: routineN = 'pulay_mixing'
300 :
301 246 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: cc_mix
302 : INTEGER :: handle, i, iatom, ib, ibb, ig, ispin, &
303 : jb, n1, n2, natom, nb, nb1, nbuffer, &
304 : ng, nspin
305 : LOGICAL :: gapw
306 : REAL(dp) :: alpha_eff, alpha_kerker, alpha_pulay, &
307 : beta, f_mix, inv_err, norm, &
308 : norm_c_inv, res_norm, vol
309 246 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, ev
310 246 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b, c, c_inv, cpc_h_mix, cpc_s_mix
311 246 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
312 : TYPE(dft_control_type), POINTER :: dft_control
313 246 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
314 246 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
315 :
316 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
317 246 : CPASSERT(ASSOCIATED(mixing_store%rhoin_buffer))
318 246 : NULLIFY (dft_control, rho_atom, rho_g, kerker_eff)
319 246 : CALL timeset(routineN, handle)
320 :
321 246 : CALL get_qs_env(qs_env, dft_control=dft_control)
322 246 : CALL qs_rho_get(rho, rho_g=rho_g)
323 246 : nspin = SIZE(rho_g, 1)
324 246 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
325 246 : vol = rho_g(1)%pw_grid%vol
326 :
327 246 : alpha_kerker = mixing_store%alpha ! default, overridden per spin in loop
328 246 : beta = mixing_store%pulay_beta
329 246 : alpha_pulay = mixing_store%pulay_alpha
330 246 : nbuffer = mixing_store%nbuffer
331 246 : gapw = dft_control%qs_control%gapw
332 246 : IF (gapw) THEN
333 12 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
334 12 : natom = SIZE(rho_atom)
335 : END IF
336 :
337 738 : ALLOCATE (cc_mix(ng))
338 :
339 246 : IF (nspin == 2) THEN
340 14 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
341 96782 : DO ig = 1, ng
342 96782 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
343 : END DO
344 14 : IF (gapw .AND. mixing_store%gmix_p) THEN
345 0 : DO iatom = 1, natom
346 0 : IF (mixing_store%paw(iatom)) THEN
347 0 : rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
348 0 : rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
349 : mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
350 0 : mixing_store%cpc_h_in(iatom, 2)%r_coef
351 : mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
352 0 : mixing_store%cpc_s_in(iatom, 2)%r_coef
353 : END IF
354 : END DO
355 : END IF
356 : END IF
357 :
358 506 : DO ispin = 1, nspin
359 : ! Select spin-channel-specific mixing parameters
360 260 : IF (ispin >= 2 .AND. nspin >= 2) THEN
361 14 : alpha_eff = mixing_store%alpha_mag
362 14 : kerker_eff => mixing_store%kerker_factor_mag
363 : ELSE
364 246 : alpha_eff = mixing_store%alpha
365 246 : kerker_eff => mixing_store%kerker_factor
366 : END IF
367 260 : alpha_kerker = alpha_eff
368 260 : ib = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
369 :
370 260 : mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) + 1
371 260 : nb = MIN(mixing_store%ncall_p(ispin), nbuffer)
372 260 : ibb = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
373 :
374 9070532 : mixing_store%res_buffer(ib, ispin)%cc(:) = CMPLX(0._dp, 0._dp, KIND=dp)
375 260 : norm = 0.0_dp
376 1306628 : IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
377 260 : res_norm = 0.0_dp
378 9070532 : DO ig = 1, ng
379 9070272 : f_mix = kerker_eff(ig)
380 : mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(rho_g(ispin)%array(ig) - &
381 9070272 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
382 : res_norm = res_norm + &
383 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)*REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
384 9070532 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))*AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
385 : END DO
386 :
387 260 : CALL para_env%sum(res_norm)
388 260 : res_norm = SQRT(res_norm)
389 :
390 260 : IF (res_norm < 1.E-14_dp) THEN
391 0 : mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) - 1
392 : ELSE
393 260 : nb1 = nb + 1
394 260 : IF (ALLOCATED(b)) DEALLOCATE (b)
395 1040 : ALLOCATE (b(nb1, nb1))
396 8436 : b = 0.0_dp
397 260 : IF (ALLOCATED(c)) DEALLOCATE (c)
398 1040 : ALLOCATE (c(nb, nb))
399 5936 : c = 0.0_dp
400 260 : IF (ALLOCATED(c_inv)) DEALLOCATE (c_inv)
401 780 : ALLOCATE (c_inv(nb, nb))
402 5936 : c_inv = 0.0_dp
403 260 : IF (ALLOCATED(alpha_c)) DEALLOCATE (alpha_c)
404 780 : ALLOCATE (alpha_c(nb))
405 1250 : alpha_c = 0.0_dp
406 260 : norm_c_inv = 0.0_dp
407 260 : IF (ALLOCATED(ev)) DEALLOCATE (ev)
408 780 : ALLOCATE (ev(nb1))
409 1510 : ev = 0.0_dp
410 :
411 : END IF
412 :
413 1250 : DO jb = 1, nb
414 990 : mixing_store%pulay_matrix(jb, ib, ispin) = 0.0_dp
415 35805150 : DO ig = 1, ng
416 35804160 : f_mix = mixing_store%special_metric(ig)
417 : mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin) + &
418 : f_mix*(REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp) &
419 : *REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
420 : AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
421 35805150 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig)))
422 : END DO
423 990 : CALL para_env%sum(mixing_store%pulay_matrix(jb, ib, ispin))
424 990 : mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)*vol*vol
425 1250 : mixing_store%pulay_matrix(ib, jb, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)
426 : END DO
427 :
428 260 : IF (nb == 1 .OR. res_norm < 1.E-14_dp) THEN
429 1306406 : DO ig = 1, ng
430 1306368 : f_mix = alpha_kerker*kerker_eff(ig)
431 : cc_mix(ig) = rho_g(ispin)%array(ig) - &
432 1306368 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
433 : rho_g(ispin)%array(ig) = f_mix*cc_mix(ig) + &
434 1306368 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
435 1306406 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
436 : END DO
437 38 : IF (mixing_store%gmix_p) THEN
438 2 : IF (gapw) THEN
439 18 : DO iatom = 1, natom
440 18 : IF (mixing_store%paw(iatom)) THEN
441 : mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
442 1850 : mixing_store%cpc_h_in(iatom, ispin)%r_coef
443 : mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
444 1850 : mixing_store%cpc_s_in(iatom, ispin)%r_coef
445 :
446 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
447 1850 : (1.0_dp - alpha_kerker)*mixing_store%cpc_h_in(iatom, ispin)%r_coef
448 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
449 1850 : (1.0_dp - alpha_kerker)*mixing_store%cpc_s_in(iatom, ispin)%r_coef
450 :
451 926 : mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
452 926 : mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
453 :
454 1850 : mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
455 1850 : mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
456 : END IF
457 : END DO
458 : END IF
459 : END IF
460 : ELSE
461 :
462 5822 : b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb, ispin)
463 5822 : c(1:nb, 1:nb) = b(1:nb, 1:nb)
464 1174 : b(nb1, 1:nb) = -1.0_dp
465 1174 : b(1:nb, nb1) = -1.0_dp
466 222 : b(nb1, nb1) = 0.0_dp
467 :
468 222 : CALL invert_matrix(c, c_inv, inv_err, improve=.TRUE.)
469 1174 : alpha_c = 0.0_dp
470 1174 : DO i = 1, nb
471 5822 : DO jb = 1, nb
472 4648 : alpha_c(i) = alpha_c(i) + c_inv(jb, i)
473 5600 : norm_c_inv = norm_c_inv + c_inv(jb, i)
474 : END DO
475 : END DO
476 1174 : alpha_c(1:nb) = alpha_c(1:nb)/norm_c_inv
477 7764126 : cc_mix = CMPLX(0._dp, 0._dp, KIND=dp)
478 1174 : DO jb = 1, nb
479 34498966 : DO ig = 1, ng
480 : cc_mix(ig) = cc_mix(ig) + &
481 : alpha_c(jb)*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) + &
482 34498744 : mixing_store%pulay_beta*mixing_store%res_buffer(jb, ispin)%cc(ig))
483 : END DO
484 : END DO
485 7764126 : mixing_store%rhoin_buffer(ibb, ispin)%cc = CMPLX(0._dp, 0._dp, KIND=dp)
486 222 : IF (alpha_pulay > 0.0_dp) THEN
487 233290 : DO ig = 1, ng
488 233280 : f_mix = alpha_pulay*kerker_eff(ig)
489 : rho_g(ispin)%array(ig) = f_mix*rho_g(ispin)%array(ig) + &
490 233280 : (1.0_dp - f_mix)*cc_mix(ig)
491 233290 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
492 : END DO
493 : ELSE
494 7530836 : DO ig = 1, ng
495 7530624 : rho_g(ispin)%array(ig) = cc_mix(ig)
496 7530836 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
497 : END DO
498 : END IF
499 :
500 222 : IF (mixing_store%gmix_p .AND. gapw) THEN
501 10 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
502 90 : DO iatom = 1, natom
503 90 : IF (mixing_store%paw(iatom)) THEN
504 10 : n1 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 1)
505 10 : n2 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 2)
506 40 : ALLOCATE (cpc_h_mix(n1, n2))
507 30 : ALLOCATE (cpc_s_mix(n1, n2))
508 :
509 : mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
510 9250 : mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef
511 : mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
512 9250 : mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef
513 4630 : cpc_h_mix = 0.0_dp
514 4630 : cpc_s_mix = 0.0_dp
515 48 : DO jb = 1, nb
516 : cpc_h_mix(:, :) = cpc_h_mix(:, :) + &
517 : alpha_c(jb)*mixing_store%cpc_h_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
518 17594 : alpha_c(jb)*beta*mixing_store%cpc_h_res_buffer(jb, iatom, ispin)%r_coef(:, :)
519 : cpc_s_mix(:, :) = cpc_s_mix(:, :) + &
520 : alpha_c(jb)*mixing_store%cpc_s_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
521 17604 : alpha_c(jb)*beta*mixing_store%cpc_s_res_buffer(jb, iatom, ispin)%r_coef(:, :)
522 : END DO
523 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
524 4630 : (1.0_dp - alpha_pulay)*cpc_h_mix
525 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
526 4630 : (1.0_dp - alpha_pulay)*cpc_s_mix
527 9250 : mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
528 9250 : mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
529 10 : DEALLOCATE (cpc_h_mix)
530 10 : DEALLOCATE (cpc_s_mix)
531 : END IF
532 : END DO
533 : END IF
534 : END IF ! nb==1
535 :
536 506 : IF (res_norm > 1.E-14_dp) THEN
537 260 : DEALLOCATE (b)
538 260 : DEALLOCATE (ev)
539 260 : DEALLOCATE (c, c_inv, alpha_c)
540 : END IF
541 :
542 : END DO ! ispin
543 :
544 246 : DEALLOCATE (cc_mix)
545 :
546 246 : IF (nspin == 2) THEN
547 14 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
548 96782 : DO ig = 1, ng
549 96782 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
550 : END DO
551 14 : IF (gapw .AND. mixing_store%gmix_p) THEN
552 0 : DO iatom = 1, natom
553 0 : IF (mixing_store%paw(iatom)) THEN
554 0 : rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
555 0 : rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
556 : mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
557 0 : mixing_store%cpc_h_in(iatom, 2)%r_coef
558 : mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
559 0 : mixing_store%cpc_s_in(iatom, 2)%r_coef
560 : END IF
561 : END DO
562 : END IF
563 :
564 : END IF
565 :
566 246 : CALL timestop(handle)
567 :
568 738 : END SUBROUTINE pulay_mixing
569 :
570 : ! **************************************************************************************************
571 : !> \brief Broyden Mixing using as metrics for the residual the Kerer damping factor
572 : !> The mixing is applied directly on the density expansion in reciprocal space
573 : !> \param qs_env ...
574 : !> \param mixing_store ...
575 : !> \param rho ...
576 : !> \param para_env ...
577 : !> \par History
578 : !> 03.2009
579 : !> \author MI
580 : ! **************************************************************************************************
581 :
582 1752 : SUBROUTINE broyden_mixing(qs_env, mixing_store, rho, para_env)
583 :
584 : TYPE(qs_environment_type), POINTER :: qs_env
585 : TYPE(mixing_storage_type), POINTER :: mixing_store
586 : TYPE(qs_rho_type), POINTER :: rho
587 : TYPE(mp_para_env_type), POINTER :: para_env
588 :
589 : CHARACTER(len=*), PARAMETER :: routineN = 'broyden_mixing'
590 :
591 : COMPLEX(dp) :: cc_mix
592 1752 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: res_rho
593 : INTEGER :: handle, i, iatom, ib, ig, ispin, j, jb, &
594 : kb, n1, n2, natom, nb, nbuffer, ng, &
595 : nspin
596 : LOGICAL :: gapw, only_kerker
597 : REAL(dp) :: alpha, alpha_eff, dcpc_h_res, &
598 : dcpc_s_res, delta_norm, f_mix, &
599 : inv_err, res_norm, rho_norm, valh, &
600 : vals, w0
601 1752 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g
602 1752 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
603 1752 : REAL(dp), DIMENSION(:), POINTER :: kerker_eff
604 : TYPE(dft_control_type), POINTER :: dft_control
605 1752 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
606 1752 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
607 :
608 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
609 1752 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
610 1752 : CPASSERT(ASSOCIATED(mixing_store%rhoin_old))
611 1752 : CPASSERT(ASSOCIATED(mixing_store%drho_buffer))
612 1752 : NULLIFY (dft_control, rho_atom, rho_g, kerker_eff)
613 1752 : CALL timeset(routineN, handle)
614 :
615 1752 : CALL get_qs_env(qs_env, dft_control=dft_control)
616 1752 : CALL qs_rho_get(rho, rho_g=rho_g)
617 :
618 1752 : nspin = SIZE(rho_g, 1)
619 1752 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
620 :
621 1752 : alpha = mixing_store%alpha
622 1752 : w0 = mixing_store%broy_w0
623 1752 : nbuffer = mixing_store%nbuffer
624 1752 : gapw = dft_control%qs_control%gapw
625 :
626 5256 : ALLOCATE (res_rho(ng))
627 :
628 1752 : mixing_store%ncall = mixing_store%ncall + 1
629 1752 : IF (mixing_store%ncall == 1) THEN
630 : only_kerker = .TRUE.
631 : ELSE
632 1536 : only_kerker = .FALSE.
633 1536 : nb = MIN(mixing_store%ncall - 1, nbuffer)
634 1536 : ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1
635 : END IF
636 :
637 1752 : IF (gapw) THEN
638 278 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
639 278 : natom = SIZE(rho_atom)
640 : ELSE
641 : natom = 0
642 : END IF
643 :
644 1752 : IF (nspin == 2) THEN
645 288 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
646 15888060 : DO ig = 1, ng
647 15888060 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
648 : END DO
649 288 : IF (gapw .AND. mixing_store%gmix_p) THEN
650 558 : DO iatom = 1, natom
651 558 : IF (mixing_store%paw(iatom)) THEN
652 375472 : rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
653 375472 : rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
654 : mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
655 375472 : mixing_store%cpc_h_in(iatom, 2)%r_coef
656 : mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
657 375472 : mixing_store%cpc_s_in(iatom, 2)%r_coef
658 : END IF
659 : END DO
660 : END IF
661 : END IF
662 :
663 3792 : DO ispin = 1, nspin
664 : ! Select spin-channel-specific mixing parameters
665 2040 : IF (ispin >= 2 .AND. nspin >= 2) THEN
666 288 : alpha_eff = mixing_store%alpha_mag
667 288 : kerker_eff => mixing_store%kerker_factor_mag
668 : ELSE
669 1752 : alpha_eff = mixing_store%alpha
670 1752 : kerker_eff => mixing_store%kerker_factor
671 : END IF
672 134962588 : res_rho = z_zero
673 134962588 : DO ig = 1, ng
674 134962588 : res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
675 : END DO
676 :
677 3792 : IF (only_kerker) THEN
678 12180866 : DO ig = 1, ng
679 12180612 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
680 12180612 : f_mix = alpha_eff*kerker_eff(ig)
681 12180612 : rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
682 12180612 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
683 12180866 : mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
684 : END DO
685 :
686 254 : IF (mixing_store%gmix_p) THEN
687 24 : IF (gapw) THEN
688 216 : DO iatom = 1, natom
689 216 : IF (mixing_store%paw(iatom)) THEN
690 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
691 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef
692 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
693 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef
694 :
695 : rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
696 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
697 : rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
698 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha_eff)
699 :
700 122972 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
701 122972 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
702 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
703 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
704 : END IF
705 : END DO
706 : END IF
707 : END IF
708 : ELSE
709 :
710 7144 : ALLOCATE (a(nb, nb))
711 65046 : a = 0.0_dp
712 5358 : ALLOCATE (b(nb, nb))
713 65046 : b = 0.0_dp
714 5358 : ALLOCATE (c(nb))
715 10320 : c = 0.0_dp
716 3572 : ALLOCATE (g(nb))
717 10320 : g = 0.0_dp
718 :
719 1786 : delta_norm = 0.0_dp
720 1786 : res_norm = 0.0_dp
721 1786 : rho_norm = 0.0_dp
722 122781722 : DO ig = 1, ng
723 122779936 : mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
724 122779936 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
725 : res_norm = res_norm + &
726 : REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + &
727 122779936 : AIMAG(res_rho(ig))*AIMAG(res_rho(ig))
728 : delta_norm = delta_norm + &
729 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)* &
730 : REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
731 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))* &
732 122779936 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
733 : rho_norm = rho_norm + &
734 : REAL(rho_g(ispin)%array(ig), dp)*REAL(rho_g(ispin)%array(ig), dp) + &
735 122781722 : AIMAG(rho_g(ispin)%array(ig))*AIMAG(rho_g(ispin)%array(ig))
736 : END DO
737 122781722 : DO ig = 1, ng
738 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
739 : mixing_store%rhoin(ispin)%cc(ig) - &
740 122781722 : mixing_store%rhoin_old(ispin)%cc(ig)
741 : END DO
742 1786 : CALL para_env%sum(delta_norm)
743 1786 : delta_norm = SQRT(delta_norm)
744 1786 : CALL para_env%sum(res_norm)
745 1786 : res_norm = SQRT(res_norm)
746 1786 : CALL para_env%sum(rho_norm)
747 1786 : rho_norm = SQRT(rho_norm)
748 :
749 1786 : IF (res_norm > 1.E-14_dp) THEN
750 122781722 : mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm
751 122781722 : mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm
752 :
753 1786 : IF (mixing_store%gmix_p .AND. gapw) THEN
754 1188 : DO iatom = 1, natom
755 1188 : IF (mixing_store%paw(iatom)) THEN
756 860 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
757 860 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
758 23912 : DO i = 1, n2
759 642788 : DO j = 1, n1
760 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
761 : (mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
762 618876 : mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
763 : dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
764 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
765 618876 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
766 : mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
767 618876 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
768 :
769 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
770 : (mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
771 618876 : mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
772 : dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
773 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
774 618876 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
775 : mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
776 618876 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
777 :
778 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
779 : alpha_eff*dcpc_h_res + &
780 618876 : mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
781 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
782 : alpha_eff*dcpc_s_res + &
783 641928 : mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
784 : END DO
785 : END DO
786 : END IF
787 : END DO
788 : END IF
789 :
790 65046 : a(:, :) = 0.0_dp
791 122781722 : DO ig = 1, ng
792 122779936 : f_mix = alpha_eff*kerker_eff(ig)
793 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
794 : f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
795 122781722 : mixing_store%drho_buffer(ib, ispin)%cc(ig)
796 : END DO
797 10320 : DO jb = 1, nb
798 41950 : DO kb = jb, nb
799 2904938942 : DO ig = 1, mixing_store%ig_max !ng
800 : a(kb, jb) = a(kb, jb) + mixing_store%p_metric(ig)*( &
801 : REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)* &
802 : REAL(mixing_store%res_buffer(kb, ispin)%cc(ig), dp) + &
803 : AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
804 2904938942 : AIMAG(mixing_store%res_buffer(kb, ispin)%cc(ig)))
805 : END DO
806 40164 : a(jb, kb) = a(kb, jb)
807 : END DO
808 : END DO
809 1786 : CALL para_env%sum(a)
810 :
811 10320 : C = 0.0_dp
812 10320 : DO jb = 1, nb
813 8534 : a(jb, jb) = w0 + a(jb, jb)
814 689264700 : DO ig = 1, mixing_store%ig_max !ng
815 : c(jb) = c(jb) + mixing_store%p_metric(ig)*( &
816 : REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)*REAL(res_rho(ig), dp) + &
817 689262914 : AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))*AIMAG(res_rho(ig)))
818 : END DO
819 : END DO
820 1786 : CALL para_env%sum(c)
821 1786 : CALL invert_matrix(a, b, inv_err)
822 :
823 1786 : CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1)
824 : ELSE
825 0 : G = 0.0_dp
826 : END IF
827 :
828 122781722 : DO ig = 1, ng
829 122779936 : cc_mix = z_zero
830 812034316 : DO jb = 1, nb
831 812034316 : cc_mix = cc_mix - G(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
832 : END DO
833 122779936 : f_mix = alpha_eff*kerker_eff(ig)
834 :
835 122779936 : IF (res_norm > 1.E-14_dp) rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
836 122779936 : f_mix*res_rho(ig) + cc_mix
837 122779936 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
838 122781722 : mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
839 : END DO
840 :
841 1786 : IF (mixing_store%gmix_p) THEN
842 132 : IF (gapw) THEN
843 1188 : DO iatom = 1, natom
844 1188 : IF (mixing_store%paw(iatom)) THEN
845 860 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
846 860 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
847 23912 : DO i = 1, n2
848 642788 : DO j = 1, n1
849 618876 : valh = 0.0_dp
850 618876 : vals = 0.0_dp
851 3304332 : DO jb = 1, nb
852 2685456 : valh = valh - G(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
853 3304332 : vals = vals - G(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
854 : END DO
855 641928 : IF (res_norm > 1.E-14_dp) THEN
856 : rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
857 : alpha_eff*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
858 618876 : mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha_eff) + valh
859 : rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
860 : alpha_eff*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
861 618876 : mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha_eff) + vals
862 : END IF
863 : END DO
864 : END DO
865 :
866 642788 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
867 642788 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
868 1284716 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
869 1284716 : mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
870 : END IF
871 : END DO
872 : END IF
873 : END IF
874 :
875 1786 : DEALLOCATE (a, b, c, g)
876 : END IF
877 :
878 : END DO ! ispin
879 1752 : IF (nspin == 2) THEN
880 288 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
881 15888060 : DO ig = 1, ng
882 15888060 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
883 : END DO
884 288 : IF (gapw .AND. mixing_store%gmix_p) THEN
885 558 : DO iatom = 1, natom
886 558 : IF (mixing_store%paw(iatom)) THEN
887 375472 : rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
888 375472 : rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
889 : mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
890 375472 : mixing_store%cpc_h_in(iatom, 2)%r_coef
891 : mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
892 375472 : mixing_store%cpc_s_in(iatom, 2)%r_coef
893 : END IF
894 : END DO
895 : END IF
896 :
897 : END IF
898 :
899 1752 : DEALLOCATE (res_rho)
900 :
901 1752 : CALL timestop(handle)
902 :
903 5256 : END SUBROUTINE broyden_mixing
904 :
905 : ! **************************************************************************************************
906 : !> \brief Multisecant scheme to perform charge density Mixing
907 : !> as preconditioner we use the Kerker damping factor
908 : !> The mixing is applied directly on the density in reciprocal space,
909 : !> therefore it affects the potentials
910 : !> on the grid but not the density matrix
911 : !> \param mixing_store ...
912 : !> \param rho ...
913 : !> \param para_env ...
914 : !> \par History
915 : !> 05.2009
916 : !> \author MI
917 : ! **************************************************************************************************
918 0 : SUBROUTINE multisecant_mixing(mixing_store, rho, para_env)
919 :
920 : TYPE(mixing_storage_type), POINTER :: mixing_store
921 : TYPE(qs_rho_type), POINTER :: rho
922 : TYPE(mp_para_env_type), POINTER :: para_env
923 :
924 : CHARACTER(len=*), PARAMETER :: routineN = 'multisecant_mixing'
925 : COMPLEX(KIND=dp), PARAMETER :: cmone = (-1.0_dp, 0.0_dp), &
926 : cone = (1.0_dp, 0.0_dp), &
927 : czero = (0.0_dp, 0.0_dp)
928 :
929 : COMPLEX(dp) :: saa, yaa
930 0 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: gn_global, pgn, res_matrix_global, &
931 0 : tmp_vec, ugn
932 0 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, res_matrix, sa, step_matrix, ya
933 0 : COMPLEX(dp), DIMENSION(:), POINTER :: gn
934 : INTEGER :: handle, ib, ib_next, ib_prev, ig, &
935 : ig_global, iig, ispin, jb, kb, nb, &
936 : nbuffer, ng, ng_global, nspin
937 : LOGICAL :: use_zgemm, use_zgemm_rev
938 : REAL(dp) :: alpha, alpha_eff, f_mix, gn_norm, gn_norm_old, inv_err, n_low, n_up, pgn_norm, &
939 : prec, r_step, reg_par, sigma_max, sigma_tilde, step_size
940 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: norm_res, norm_res_low, norm_res_up
941 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b_matrix, binv_matrix
942 0 : REAL(dp), DIMENSION(:), POINTER :: g2, kerker_eff
943 : REAL(dp), SAVE :: sigma_old = 1.0_dp
944 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
945 :
946 0 : CPASSERT(ASSOCIATED(mixing_store))
947 :
948 0 : CALL timeset(routineN, handle)
949 :
950 0 : NULLIFY (gn, rho_g)
951 :
952 0 : use_zgemm = .FALSE.
953 0 : use_zgemm_rev = .TRUE.
954 :
955 : ! prepare the parameters
956 0 : CALL qs_rho_get(rho, rho_g=rho_g)
957 :
958 0 : nspin = SIZE(rho_g, 1)
959 : ! not implemented for large grids.
960 0 : CPASSERT(rho_g(1)%pw_grid%ngpts < HUGE(ng_global))
961 0 : ng_global = INT(rho_g(1)%pw_grid%ngpts)
962 0 : ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
963 0 : alpha = mixing_store%alpha
964 :
965 0 : sigma_max = mixing_store%sigma_max
966 0 : reg_par = mixing_store%reg_par
967 0 : r_step = mixing_store%r_step
968 0 : nbuffer = mixing_store%nbuffer
969 :
970 : ! determine the step number, and multisecant iteration
971 0 : nb = MIN(mixing_store%ncall, nbuffer - 1)
972 0 : ib = MODULO(mixing_store%ncall, nbuffer) + 1
973 0 : IF (mixing_store%ncall > 0) THEN
974 0 : ib_prev = MODULO(mixing_store%ncall - 1, nbuffer) + 1
975 : ELSE
976 : ib_prev = 0
977 : END IF
978 0 : mixing_store%ncall = mixing_store%ncall + 1
979 0 : ib_next = MODULO(mixing_store%ncall, nbuffer) + 1
980 :
981 : ! compute the residual gn and its norm gn_norm
982 0 : DO ispin = 1, nspin
983 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
984 0 : gn_norm = 0.0_dp
985 0 : DO ig = 1, ng
986 0 : gn(ig) = (rho_g(ispin)%array(ig) - mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
987 : gn_norm = gn_norm + &
988 0 : REAL(gn(ig), dp)*REAL(gn(ig), dp) + AIMAG(gn(ig))*AIMAG(gn(ig))
989 : END DO
990 0 : CALL para_env%sum(gn_norm)
991 0 : gn_norm = SQRT(gn_norm)
992 0 : mixing_store%norm_res_buffer(ib, ispin) = gn_norm
993 : END DO
994 :
995 0 : IF (nb == 0) THEN
996 : !simple mixing
997 0 : DO ispin = 1, nspin
998 0 : IF (ispin >= 2 .AND. nspin >= 2) THEN
999 0 : alpha_eff = mixing_store%alpha_mag
1000 0 : kerker_eff => mixing_store%kerker_factor_mag
1001 : ELSE
1002 0 : alpha_eff = mixing_store%alpha
1003 0 : kerker_eff => mixing_store%kerker_factor
1004 : END IF
1005 0 : DO ig = 1, ng
1006 0 : f_mix = alpha_eff*kerker_eff(ig)
1007 : rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(1, ispin)%cc(ig) + &
1008 0 : f_mix*mixing_store%res_buffer(1, ispin)%cc(ig)
1009 0 : mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1010 : END DO
1011 : END DO
1012 0 : CALL timestop(handle)
1013 0 : RETURN
1014 : END IF
1015 :
1016 : ! allocate temporary arrays
1017 : ! step_matrix S ngxnb
1018 0 : ALLOCATE (step_matrix(ng, nb))
1019 : ! res_matrix Y ngxnb
1020 0 : ALLOCATE (res_matrix(ng, nb))
1021 : ! matrix A nbxnb
1022 0 : ALLOCATE (a_matrix(nb, ng_global))
1023 : ! PSI nb vector of norms
1024 0 : ALLOCATE (norm_res(nb))
1025 0 : ALLOCATE (norm_res_low(nb))
1026 0 : ALLOCATE (norm_res_up(nb))
1027 :
1028 : ! matrix B nbxnb
1029 0 : ALLOCATE (b_matrix(nb, nb))
1030 : ! matrix B_inv nbxnb
1031 0 : ALLOCATE (binv_matrix(nb, nb))
1032 :
1033 0 : ALLOCATE (gn_global(ng_global))
1034 0 : ALLOCATE (res_matrix_global(ng_global))
1035 : IF (use_zgemm) THEN
1036 : ALLOCATE (sa(ng, ng_global))
1037 : ALLOCATE (ya(ng, ng_global))
1038 : END IF
1039 : IF (use_zgemm_rev) THEN
1040 0 : ALLOCATE (tmp_vec(nb))
1041 : END IF
1042 0 : ALLOCATE (pgn(ng))
1043 0 : ALLOCATE (ugn(ng))
1044 :
1045 0 : DO ispin = 1, nspin
1046 : ! Select spin-channel-specific mixing parameters
1047 0 : IF (ispin >= 2 .AND. nspin >= 2) THEN
1048 0 : alpha_eff = mixing_store%alpha_mag
1049 0 : kerker_eff => mixing_store%kerker_factor_mag
1050 : ELSE
1051 0 : alpha_eff = mixing_store%alpha
1052 0 : kerker_eff => mixing_store%kerker_factor
1053 : END IF
1054 : ! generate the global vector with the present residual
1055 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1056 0 : gn_global = z_zero
1057 0 : DO ig = 1, ng
1058 0 : ig_global = mixing_store%ig_global_index(ig)
1059 0 : gn_global(ig_global) = gn(ig)
1060 : END DO
1061 0 : CALL para_env%sum(gn_global)
1062 :
1063 : ! compute steps (matrix S) and residual differences (matrix Y)
1064 : ! with respect to the present
1065 : ! step and the present residual (use stored rho_in and res_buffer)
1066 :
1067 : ! These quantities are pre-conditioned by means of Kerker factor multipliccation
1068 0 : g2 => rho_g(1)%pw_grid%gsq
1069 :
1070 0 : DO jb = 1, nb + 1
1071 0 : IF (jb < ib) THEN
1072 : kb = jb
1073 0 : ELSEIF (jb > ib) THEN
1074 0 : kb = jb - 1
1075 : ELSE
1076 : CYCLE
1077 : END IF
1078 0 : norm_res(kb) = 0.0_dp
1079 0 : norm_res_low(kb) = 0.0_dp
1080 0 : norm_res_up(kb) = 0.0_dp
1081 0 : DO ig = 1, ng
1082 :
1083 0 : prec = 1.0_dp
1084 :
1085 : step_matrix(ig, kb) = prec*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) - &
1086 0 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1087 : res_matrix(ig, kb) = (mixing_store%res_buffer(jb, ispin)%cc(ig) - &
1088 0 : mixing_store%res_buffer(ib, ispin)%cc(ig))
1089 : norm_res(kb) = norm_res(kb) + REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1090 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1091 0 : IF (g2(ig) < 4.0_dp) THEN
1092 : norm_res_low(kb) = norm_res_low(kb) + &
1093 : REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1094 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1095 : ELSE
1096 : norm_res_up(kb) = norm_res_up(kb) + &
1097 : REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1098 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1099 : END IF
1100 0 : res_matrix(ig, kb) = prec*res_matrix(ig, kb)
1101 : END DO
1102 : END DO !jb
1103 :
1104 : ! normalize each column of S and Y => Snorm Ynorm
1105 0 : CALL para_env%sum(norm_res)
1106 0 : CALL para_env%sum(norm_res_up)
1107 0 : CALL para_env%sum(norm_res_low)
1108 0 : norm_res(1:nb) = 1.0_dp/SQRT(norm_res(1:nb))
1109 : n_low = 0.0_dp
1110 : n_up = 0.0_dp
1111 0 : DO jb = 1, nb
1112 0 : n_low = n_low + norm_res_low(jb)/norm_res(jb)
1113 0 : n_up = n_up + norm_res_up(jb)/norm_res(jb)
1114 : END DO
1115 0 : DO ig = 1, ng
1116 0 : IF (g2(ig) > 4.0_dp) THEN
1117 0 : step_matrix(ig, 1:nb) = step_matrix(ig, 1:nb)*SQRT(n_low/n_up)
1118 0 : res_matrix(ig, 1:nb) = res_matrix(ig, 1:nb)*SQRT(n_low/n_up)
1119 : END IF
1120 : END DO
1121 0 : DO kb = 1, nb
1122 0 : step_matrix(1:ng, kb) = step_matrix(1:ng, kb)*norm_res(kb)
1123 0 : res_matrix(1:ng, kb) = res_matrix(1:ng, kb)*norm_res(kb)
1124 : END DO
1125 :
1126 : ! compute A as [(Ynorm^t Ynorm) + (alpha I)]^(-1) Ynorm^t
1127 : ! compute B
1128 0 : DO jb = 1, nb
1129 0 : DO kb = 1, nb
1130 0 : b_matrix(kb, jb) = 0.0_dp
1131 0 : DO ig = 1, ng
1132 : ! it is assumed that summing over all G vector gives a real, because
1133 : ! y(-G,kb) = (y(G,kb))*
1134 0 : b_matrix(kb, jb) = b_matrix(kb, jb) + REAL(res_matrix(ig, kb)*res_matrix(ig, jb), dp)
1135 : END DO
1136 : END DO
1137 : END DO
1138 :
1139 0 : CALL para_env%sum(b_matrix)
1140 0 : DO jb = 1, nb
1141 0 : b_matrix(jb, jb) = b_matrix(jb, jb) + reg_par
1142 : END DO
1143 : ! invert B
1144 0 : CALL invert_matrix(b_matrix, binv_matrix, inv_err)
1145 :
1146 : ! A = Binv Ynorm^t
1147 0 : a_matrix = z_zero
1148 0 : DO kb = 1, nb
1149 0 : res_matrix_global = z_zero
1150 0 : DO ig = 1, ng
1151 0 : ig_global = mixing_store%ig_global_index(ig)
1152 0 : res_matrix_global(ig_global) = res_matrix(ig, kb)
1153 : END DO
1154 0 : CALL para_env%sum(res_matrix_global)
1155 :
1156 0 : DO jb = 1, nb
1157 0 : DO ig = 1, ng
1158 0 : ig_global = mixing_store%ig_global_index(ig)
1159 : a_matrix(jb, ig_global) = a_matrix(jb, ig_global) + &
1160 0 : binv_matrix(jb, kb)*res_matrix_global(ig_global)
1161 : END DO
1162 : END DO
1163 : END DO
1164 0 : CALL para_env%sum(a_matrix)
1165 :
1166 : ! compute the two components of gn that will be used to update rho
1167 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1168 0 : pgn_norm = 0.0_dp
1169 :
1170 : IF (use_zgemm) THEN
1171 :
1172 : CALL zgemm("N", "N", ng, ng_global, nb, cmone, step_matrix(1, 1), ng, &
1173 : a_matrix(1, 1), nb, czero, sa(1, 1), ng)
1174 : CALL zgemm("N", "N", ng, ng_global, nb, cmone, res_matrix(1, 1), ng, &
1175 : a_matrix(1, 1), nb, czero, ya(1, 1), ng)
1176 : DO ig = 1, ng
1177 : ig_global = mixing_store%ig_global_index(ig)
1178 : ya(ig, ig_global) = ya(ig, ig_global) + z_one
1179 : END DO
1180 :
1181 : CALL zgemv("N", ng, ng_global, cone, sa(1, 1), &
1182 : ng, gn_global(1), 1, czero, pgn(1), 1)
1183 : CALL zgemv("N", ng, ng_global, cone, ya(1, 1), &
1184 : ng, gn_global(1), 1, czero, ugn(1), 1)
1185 :
1186 : DO ig = 1, ng
1187 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1188 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1189 : END DO
1190 : CALL para_env%sum(pgn_norm)
1191 : ELSEIF (use_zgemm_rev) THEN
1192 :
1193 : CALL zgemv("N", nb, ng_global, cone, a_matrix(1, 1), &
1194 0 : nb, gn_global(1), 1, czero, tmp_vec(1), 1)
1195 :
1196 : CALL zgemv("N", ng, nb, cmone, step_matrix(1, 1), ng, &
1197 0 : tmp_vec(1), 1, czero, pgn(1), 1)
1198 :
1199 : CALL zgemv("N", ng, nb, cmone, res_matrix(1, 1), ng, &
1200 0 : tmp_vec(1), 1, czero, ugn(1), 1)
1201 :
1202 0 : DO ig = 1, ng
1203 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1204 0 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1205 0 : ugn(ig) = ugn(ig) + gn(ig)
1206 : END DO
1207 0 : CALL para_env%sum(pgn_norm)
1208 :
1209 : ELSE
1210 : DO ig = 1, ng
1211 : pgn(ig) = z_zero
1212 : ugn(ig) = z_zero
1213 : ig_global = mixing_store%ig_global_index(ig)
1214 : DO iig = 1, ng_global
1215 : saa = z_zero
1216 : yaa = z_zero
1217 :
1218 : IF (ig_global == iig) yaa = z_one
1219 :
1220 : DO jb = 1, nb
1221 : saa = saa - step_matrix(ig, jb)*a_matrix(jb, iig)
1222 : yaa = yaa - res_matrix(ig, jb)*a_matrix(jb, iig)
1223 : END DO
1224 : pgn(ig) = pgn(ig) + saa*gn_global(iig)
1225 : ugn(ig) = ugn(ig) + yaa*gn_global(iig)
1226 : END DO
1227 : END DO
1228 : DO ig = 1, ng
1229 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1230 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1231 : END DO
1232 : CALL para_env%sum(pgn_norm)
1233 : END IF
1234 :
1235 0 : gn_norm = mixing_store%norm_res_buffer(ib, ispin)
1236 0 : gn_norm_old = mixing_store%norm_res_buffer(ib_prev, ispin)
1237 : IF (ib_prev /= 0) THEN
1238 0 : sigma_tilde = sigma_old*MAX(0.5_dp, MIN(2.0_dp, gn_norm_old/gn_norm))
1239 : ELSE
1240 : sigma_tilde = 0.5_dp
1241 : END IF
1242 0 : sigma_tilde = 0.1_dp
1243 : ! Step size for the unpredicted component
1244 0 : step_size = MIN(sigma_tilde, r_step*pgn_norm/gn_norm, sigma_max)
1245 0 : sigma_old = step_size
1246 :
1247 : ! update the density
1248 0 : DO ig = 1, ng
1249 0 : prec = kerker_eff(ig)
1250 : rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(ib, ispin)%cc(ig) &
1251 0 : - prec*step_size*ugn(ig) + prec*pgn(ig) ! - 0.1_dp * prec* gn(ig)
1252 0 : mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1253 : END DO
1254 :
1255 : END DO ! ispin
1256 :
1257 : ! Deallocate temporary arrays
1258 0 : DEALLOCATE (step_matrix, res_matrix)
1259 0 : DEALLOCATE (norm_res)
1260 0 : DEALLOCATE (norm_res_up)
1261 0 : DEALLOCATE (norm_res_low)
1262 0 : DEALLOCATE (a_matrix, b_matrix, binv_matrix)
1263 0 : DEALLOCATE (ugn, pgn)
1264 : IF (use_zgemm) THEN
1265 : DEALLOCATE (sa, ya)
1266 : END IF
1267 : IF (use_zgemm_rev) THEN
1268 0 : DEALLOCATE (tmp_vec)
1269 : END IF
1270 0 : DEALLOCATE (gn_global, res_matrix_global)
1271 :
1272 0 : CALL timestop(handle)
1273 :
1274 0 : END SUBROUTINE multisecant_mixing
1275 :
1276 : END MODULE qs_gspace_mixing
|
