(git:0de0cc2)
qs_gspace_mixing.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
10 
11  USE cp_control_types, ONLY: dft_control_type
12  USE kinds, ONLY: dp
13  USE mathlib, ONLY: invert_matrix
14  USE message_passing, ONLY: mp_para_env_type
15  USE pw_env_types, ONLY: pw_env_get,&
16  pw_env_type
17  USE pw_methods, ONLY: pw_axpy,&
18  pw_copy,&
19  pw_integrate_function,&
20  pw_scale,&
21  pw_transfer,&
22  pw_zero
23  USE pw_pool_types, ONLY: pw_pool_type
24  USE pw_types, ONLY: pw_c1d_gs_type,&
25  pw_r3d_rs_type
28  cp_1d_z_p_type,&
30  mixing_storage_type,&
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  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  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  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  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
86  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
87 
88  CALL timeset(routinen, handle)
89 
90  CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env)
91  IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
92  dft_control%qs_control%semi_empirical)) THEN
93 
94  cpassert(ASSOCIATED(mixing_store))
95  NULLIFY (auxbas_pw_pool, rho_atom, rho_g, rho_r, tot_rho_r)
96  CALL qs_rho_get(rho, rho_g=rho_g, rho_r=rho_r, tot_rho_r=tot_rho_r)
97 
98  nspin = SIZE(rho_g, 1)
99  nimg = dft_control%nimages
100  ng = SIZE(rho_g(1)%pw_grid%gsq)
101  cpassert(ng == SIZE(mixing_store%rhoin(1)%cc))
102 
103  alpha = mixing_store%alpha
104  gapw = dft_control%qs_control%gapw
105 
106  IF (nspin == 2) THEN
107  CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
108  CALL auxbas_pw_pool%create_pw(rho_tmp)
109  CALL pw_zero(rho_tmp)
110  CALL pw_copy(rho_g(1), rho_tmp)
111  CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
112  CALL pw_axpy(rho_tmp, rho_g(2), -1.0_dp)
113  CALL pw_scale(rho_g(2), -1.0_dp)
114  END IF
115 
116  IF (iter_count + 1 <= mixing_store%nskip_mixing) THEN
117  ! skip mixing
118  DO ispin = 1, nspin
119  DO ig = 1, ng
120  mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
121  END DO
122  END DO
123  IF (mixing_store%gmix_p .AND. gapw) THEN
124  CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
125  natom = SIZE(rho_atom)
126  DO ispin = 1, nspin
127  DO iatom = 1, natom
128  IF (mixing_store%paw(iatom)) THEN
129  mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
130  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  mixing_store%iter_method = "NoMix"
137  IF (nspin == 2) THEN
138  CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
139  CALL pw_scale(rho_g(1), 0.5_dp)
140  CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
141  CALL pw_scale(rho_g(2), -1.0_dp)
142  CALL auxbas_pw_pool%give_back_pw(rho_tmp)
143  END IF
144  CALL timestop(handle)
145  RETURN
146  END IF
147 
148  IF ((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) THEN
149  CALL gmix_potential_only(qs_env, mixing_store, rho)
150  mixing_store%iter_method = "Kerker"
151  ELSEIF (mixing_method == gspace_mixing_nr) THEN
152  CALL gmix_potential_only(qs_env, mixing_store, rho)
153  mixing_store%iter_method = "Kerker"
154  ELSEIF (mixing_method == pulay_mixing_nr) THEN
155  CALL pulay_mixing(qs_env, mixing_store, rho, para_env)
156  mixing_store%iter_method = "Pulay"
157  ELSEIF (mixing_method == broyden_mixing_nr) THEN
158  CALL broyden_mixing(qs_env, mixing_store, rho, para_env)
159  mixing_store%iter_method = "Broy."
160  ELSEIF (mixing_method == broyden_mixing_new_nr) THEN
161  cpassert(.NOT. gapw)
162  CALL broyden_mixing_new(mixing_store, rho, para_env)
163  mixing_store%iter_method = "Broy."
164  ELSEIF (mixing_method == multisecant_mixing_nr) THEN
165  cpassert(.NOT. gapw)
166  CALL multisecant_mixing(mixing_store, rho, para_env)
167  mixing_store%iter_method = "MSec."
168  END IF
169 
170  IF (nspin == 2) THEN
171  CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
172  CALL pw_scale(rho_g(1), 0.5_dp)
173  CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
174  CALL pw_scale(rho_g(2), -1.0_dp)
175  CALL auxbas_pw_pool%give_back_pw(rho_tmp)
176  END IF
177 
178  DO ispin = 1, nspin
179  CALL pw_transfer(rho_g(ispin), rho_r(ispin))
180  tot_rho_r(ispin) = pw_integrate_function(rho_r(ispin), isign=-1)
181  END DO
182  END IF
183 
184  CALL timestop(handle)
185 
186  END SUBROUTINE gspace_mixing
187 
188 ! **************************************************************************************************
189 !> \brief G-space mixing performed via the Kerker damping on the density on the grid
190 !> thus affecting only the caluclation of the potential, but not the denisity matrix
191 !> \param qs_env ...
192 !> \param mixing_store ...
193 !> \param rho ...
194 !> \par History
195 !> 02.2009
196 !> \author MI
197 ! **************************************************************************************************
198 
199  SUBROUTINE gmix_potential_only(qs_env, mixing_store, rho)
200 
201  TYPE(qs_environment_type), POINTER :: qs_env
202  TYPE(mixing_storage_type), POINTER :: mixing_store
203  TYPE(qs_rho_type), POINTER :: rho
204 
205  CHARACTER(len=*), PARAMETER :: routinen = 'gmix_potential_only'
206 
207  COMPLEX(dp), DIMENSION(:), POINTER :: cc_new
208  INTEGER :: handle, iatom, ig, ispin, natom, ng, &
209  nspin
210  LOGICAL :: gapw
211  REAL(dp) :: alpha, f_mix
212  TYPE(dft_control_type), POINTER :: dft_control
213  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
214  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
215 
216  cpassert(ASSOCIATED(mixing_store%rhoin))
217  cpassert(ASSOCIATED(mixing_store%kerker_factor))
218 
219  CALL timeset(routinen, handle)
220  NULLIFY (cc_new, dft_control, rho_atom, rho_g)
221 
222  CALL get_qs_env(qs_env, dft_control=dft_control)
223  CALL qs_rho_get(rho, rho_g=rho_g)
224 
225  nspin = SIZE(rho_g, 1)
226  ng = SIZE(rho_g(1)%pw_grid%gsq)
227 
228  gapw = dft_control%qs_control%gapw
229  alpha = mixing_store%alpha
230 
231  DO ispin = 1, nspin
232  cc_new => rho_g(ispin)%array
233  DO ig = 1, mixing_store%ig_max ! ng
234  f_mix = mixing_store%alpha*mixing_store%kerker_factor(ig)
235  cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
236  mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
237  END DO
238  DO ig = mixing_store%ig_max + 1, ng
239  f_mix = mixing_store%alpha
240  cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
241  mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
242  END DO
243 
244  END DO
245 
246  IF (mixing_store%gmix_p .AND. gapw) THEN
247  CALL get_qs_env(qs_env=qs_env, &
248  rho_atom_set=rho_atom)
249  natom = SIZE(rho_atom)
250  DO ispin = 1, nspin
251  DO iatom = 1, natom
252  IF (mixing_store%paw(iatom)) THEN
253  rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
254  mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha)
255  rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
256  mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha)
257  mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
258  mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
259  END IF
260  END DO
261  END DO
262  END IF
263 
264  CALL timestop(handle)
265 
266  END SUBROUTINE gmix_potential_only
267 
268 ! **************************************************************************************************
269 !> \brief Pulay Mixing using as metrics for the residual the Kerer damping factor
270 !> The mixing is applied directly on the density in reciprocal space,
271 !> therefore it affects the potentials
272 !> on the grid but not the density matrix
273 !> \param qs_env ...
274 !> \param mixing_store ...
275 !> \param rho ...
276 !> \param para_env ...
277 !> \par History
278 !> 03.2009
279 !> \author MI
280 ! **************************************************************************************************
281 
282  SUBROUTINE pulay_mixing(qs_env, mixing_store, rho, para_env)
283 
284  TYPE(qs_environment_type), POINTER :: qs_env
285  TYPE(mixing_storage_type), POINTER :: mixing_store
286  TYPE(qs_rho_type), POINTER :: rho
287  TYPE(mp_para_env_type), POINTER :: para_env
288 
289  CHARACTER(len=*), PARAMETER :: routinen = 'pulay_mixing'
290 
291  COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: cc_mix
292  INTEGER :: handle, i, iatom, ib, ibb, ig, ispin, &
293  jb, n1, n2, natom, nb, nb1, nbuffer, &
294  ng, nspin
295  LOGICAL :: gapw
296  REAL(dp) :: alpha_kerker, alpha_pulay, beta, f_mix, &
297  inv_err, norm, norm_c_inv, res_norm, &
298  vol
299  REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, ev
300  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, c, c_inv, cpc_h_mix, cpc_s_mix
301  TYPE(dft_control_type), POINTER :: dft_control
302  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
303  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
304 
305  cpassert(ASSOCIATED(mixing_store%res_buffer))
306  cpassert(ASSOCIATED(mixing_store%rhoin_buffer))
307  NULLIFY (dft_control, rho_atom, rho_g)
308  CALL timeset(routinen, handle)
309 
310  CALL get_qs_env(qs_env, dft_control=dft_control)
311  CALL qs_rho_get(rho, rho_g=rho_g)
312  nspin = SIZE(rho_g, 1)
313  ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
314  vol = rho_g(1)%pw_grid%vol
315 
316  alpha_kerker = mixing_store%alpha
317  beta = mixing_store%pulay_beta
318  alpha_pulay = mixing_store%pulay_alpha
319  nbuffer = mixing_store%nbuffer
320  gapw = dft_control%qs_control%gapw
321  IF (gapw) THEN
322  CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
323  natom = SIZE(rho_atom)
324  END IF
325 
326  ALLOCATE (cc_mix(ng))
327 
328  IF (nspin == 2) THEN
329  CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
330  DO ig = 1, ng
331  mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
332  END DO
333  IF (gapw .AND. mixing_store%gmix_p) THEN
334  DO iatom = 1, natom
335  IF (mixing_store%paw(iatom)) THEN
336  rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
337  rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
338  mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
339  mixing_store%cpc_h_in(iatom, 2)%r_coef
340  mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
341  mixing_store%cpc_s_in(iatom, 2)%r_coef
342  END IF
343  END DO
344  END IF
345  END IF
346 
347  DO ispin = 1, nspin
348  ib = modulo(mixing_store%ncall_p(ispin), nbuffer) + 1
349 
350  mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) + 1
351  nb = min(mixing_store%ncall_p(ispin), nbuffer)
352  ibb = modulo(mixing_store%ncall_p(ispin), nbuffer) + 1
353 
354  mixing_store%res_buffer(ib, ispin)%cc(:) = cmplx(0._dp, 0._dp, kind=dp)
355  norm = 0.0_dp
356  IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
357  res_norm = 0.0_dp
358  DO ig = 1, ng
359  f_mix = mixing_store%kerker_factor(ig)
360  mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(rho_g(ispin)%array(ig) - &
361  mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
362  res_norm = res_norm + &
363  REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)*REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
364  aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))*aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))
365  END DO
366 
367  CALL para_env%sum(res_norm)
368  res_norm = sqrt(res_norm)
369 
370  IF (res_norm < 1.e-14_dp) THEN
371  mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) - 1
372  ELSE
373  nb1 = nb + 1
374  ALLOCATE (a(nb1, nb1))
375  a = 0.0_dp
376  ALLOCATE (b(nb1, nb1))
377  b = 0.0_dp
378  ALLOCATE (c(nb, nb))
379  c = 0.0_dp
380  ALLOCATE (c_inv(nb, nb))
381  c_inv = 0.0_dp
382  ALLOCATE (alpha_c(nb))
383  alpha_c = 0.0_dp
384  norm_c_inv = 0.0_dp
385  ALLOCATE (ev(nb1))
386  ev = 0.0_dp
387 
388  END IF
389 
390  DO jb = 1, nb
391  mixing_store%pulay_matrix(jb, ib, ispin) = 0.0_dp
392  DO ig = 1, ng
393  f_mix = mixing_store%special_metric(ig)
394  mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin) + &
395  f_mix*(real(mixing_store%res_buffer(jb, ispin)%cc(ig), dp) &
396  *real(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
397  aimag(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
398  aimag(mixing_store%res_buffer(ib, ispin)%cc(ig)))
399  END DO
400  CALL para_env%sum(mixing_store%pulay_matrix(jb, ib, ispin))
401  mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)*vol*vol
402  mixing_store%pulay_matrix(ib, jb, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)
403  END DO
404 
405  IF (nb == 1 .OR. res_norm < 1.e-14_dp) THEN
406  DO ig = 1, ng
407  f_mix = alpha_kerker*mixing_store%kerker_factor(ig)
408  cc_mix(ig) = rho_g(ispin)%array(ig) - &
409  mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
410  rho_g(ispin)%array(ig) = f_mix*cc_mix(ig) + &
411  mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
412  mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
413  END DO
414  IF (mixing_store%gmix_p) THEN
415  IF (gapw) THEN
416  DO iatom = 1, natom
417  IF (mixing_store%paw(iatom)) THEN
418  mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
419  mixing_store%cpc_h_in(iatom, ispin)%r_coef
420  mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
421  mixing_store%cpc_s_in(iatom, ispin)%r_coef
422 
423  rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
424  (1.0_dp - alpha_kerker)*mixing_store%cpc_h_in(iatom, ispin)%r_coef
425  rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
426  (1.0_dp - alpha_kerker)*mixing_store%cpc_s_in(iatom, ispin)%r_coef
427 
428  mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
429  mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
430 
431  mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
432  mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
433  END IF
434  END DO
435  END IF
436  END IF
437  ELSE
438 
439  b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb, ispin)
440  c(1:nb, 1:nb) = b(1:nb, 1:nb)
441  b(nb1, 1:nb) = -1.0_dp
442  b(1:nb, nb1) = -1.0_dp
443  b(nb1, nb1) = 0.0_dp
444 
445  CALL invert_matrix(c, c_inv, inv_err, improve=.true.)
446  alpha_c = 0.0_dp
447  DO i = 1, nb
448  DO jb = 1, nb
449  alpha_c(i) = alpha_c(i) + c_inv(jb, i)
450  norm_c_inv = norm_c_inv + c_inv(jb, i)
451  END DO
452  END DO
453  alpha_c(1:nb) = alpha_c(1:nb)/norm_c_inv
454  cc_mix = cmplx(0._dp, 0._dp, kind=dp)
455  DO jb = 1, nb
456  DO ig = 1, ng
457  cc_mix(ig) = cc_mix(ig) + &
458  alpha_c(jb)*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) + &
459  mixing_store%pulay_beta*mixing_store%res_buffer(jb, ispin)%cc(ig))
460  END DO
461  END DO
462  mixing_store%rhoin_buffer(ibb, ispin)%cc = cmplx(0._dp, 0._dp, kind=dp)
463  IF (alpha_pulay > 0.0_dp) THEN
464  DO ig = 1, ng
465  f_mix = alpha_pulay*mixing_store%kerker_factor(ig)
466  rho_g(ispin)%array(ig) = f_mix*rho_g(ispin)%array(ig) + &
467  (1.0_dp - f_mix)*cc_mix(ig)
468  mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
469  END DO
470  ELSE
471  DO ig = 1, ng
472  rho_g(ispin)%array(ig) = cc_mix(ig)
473  mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
474  END DO
475  END IF
476 
477  IF (mixing_store%gmix_p .AND. gapw) THEN
478  CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
479  DO iatom = 1, natom
480  IF (mixing_store%paw(iatom)) THEN
481  n1 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 1)
482  n2 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 2)
483  ALLOCATE (cpc_h_mix(n1, n2))
484  ALLOCATE (cpc_s_mix(n1, n2))
485 
486  mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
487  mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef
488  mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
489  mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef
490  cpc_h_mix = 0.0_dp
491  cpc_s_mix = 0.0_dp
492  DO jb = 1, nb
493  cpc_h_mix(:, :) = cpc_h_mix(:, :) + &
494  alpha_c(jb)*mixing_store%cpc_h_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
495  alpha_c(jb)*beta*mixing_store%cpc_h_res_buffer(jb, iatom, ispin)%r_coef(:, :)
496  cpc_s_mix(:, :) = cpc_s_mix(:, :) + &
497  alpha_c(jb)*mixing_store%cpc_s_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
498  alpha_c(jb)*beta*mixing_store%cpc_s_res_buffer(jb, iatom, ispin)%r_coef(:, :)
499  END DO
500  rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
501  (1.0_dp - alpha_pulay)*cpc_h_mix
502  rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
503  (1.0_dp - alpha_pulay)*cpc_s_mix
504  mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
505  mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
506  DEALLOCATE (cpc_h_mix)
507  DEALLOCATE (cpc_s_mix)
508  END IF
509  END DO
510  END IF
511  END IF ! nb==1
512 
513  IF (res_norm > 1.e-14_dp) THEN
514  DEALLOCATE (a)
515  DEALLOCATE (b)
516  DEALLOCATE (ev)
517  DEALLOCATE (c, c_inv, alpha_c)
518  END IF
519 
520  END DO ! ispin
521 
522  DEALLOCATE (cc_mix)
523 
524  IF (nspin == 2) THEN
525  CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
526  DO ig = 1, ng
527  mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
528  END DO
529  IF (gapw .AND. mixing_store%gmix_p) THEN
530  DO iatom = 1, natom
531  IF (mixing_store%paw(iatom)) THEN
532  rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
533  rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
534  mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
535  mixing_store%cpc_h_in(iatom, 2)%r_coef
536  mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
537  mixing_store%cpc_s_in(iatom, 2)%r_coef
538  END IF
539  END DO
540  END IF
541 
542  END IF
543 
544  CALL timestop(handle)
545 
546  END SUBROUTINE pulay_mixing
547 
548 ! **************************************************************************************************
549 !> \brief Broyden Mixing using as metrics for the residual the Kerer damping factor
550 !> The mixing is applied directly on the density expansion in reciprocal space
551 !> \param qs_env ...
552 !> \param mixing_store ...
553 !> \param rho ...
554 !> \param para_env ...
555 !> \par History
556 !> 03.2009
557 !> \author MI
558 ! **************************************************************************************************
559 
560  SUBROUTINE broyden_mixing(qs_env, mixing_store, rho, para_env)
561 
562  TYPE(qs_environment_type), POINTER :: qs_env
563  TYPE(mixing_storage_type), POINTER :: mixing_store
564  TYPE(qs_rho_type), POINTER :: rho
565  TYPE(mp_para_env_type), POINTER :: para_env
566 
567  CHARACTER(len=*), PARAMETER :: routinen = 'broyden_mixing'
568 
569  COMPLEX(dp) :: cc_mix
570  COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: res_rho
571  INTEGER :: handle, i, iatom, ib, ig, ispin, j, jb, &
572  kb, n1, n2, natom, nb, nbuffer, ng, &
573  nspin
574  LOGICAL :: gapw, only_kerker
575  REAL(dp) :: alpha, dcpc_h_res, dcpc_s_res, &
576  delta_norm, f_mix, inv_err, res_norm, &
577  rho_norm, valh, vals, w0
578  REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g
579  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
580  TYPE(dft_control_type), POINTER :: dft_control
581  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
582  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
583 
584  cpassert(ASSOCIATED(mixing_store%res_buffer))
585  cpassert(ASSOCIATED(mixing_store%rhoin))
586  cpassert(ASSOCIATED(mixing_store%rhoin_old))
587  cpassert(ASSOCIATED(mixing_store%drho_buffer))
588  NULLIFY (dft_control, rho_atom, rho_g)
589  CALL timeset(routinen, handle)
590 
591  CALL get_qs_env(qs_env, dft_control=dft_control)
592  CALL qs_rho_get(rho, rho_g=rho_g)
593 
594  nspin = SIZE(rho_g, 1)
595  ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
596 
597  alpha = mixing_store%alpha
598  w0 = mixing_store%broy_w0
599  nbuffer = mixing_store%nbuffer
600  gapw = dft_control%qs_control%gapw
601 
602  ALLOCATE (res_rho(ng))
603 
604  mixing_store%ncall = mixing_store%ncall + 1
605  IF (mixing_store%ncall == 1) THEN
606  only_kerker = .true.
607  ELSE
608  only_kerker = .false.
609  nb = min(mixing_store%ncall - 1, nbuffer)
610  ib = modulo(mixing_store%ncall - 2, nbuffer) + 1
611  END IF
612 
613  IF (gapw) THEN
614  CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
615  natom = SIZE(rho_atom)
616  ELSE
617  natom = 0
618  END IF
619 
620  IF (nspin == 2) THEN
621  CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
622  DO ig = 1, ng
623  mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
624  END DO
625  IF (gapw .AND. mixing_store%gmix_p) THEN
626  DO iatom = 1, natom
627  IF (mixing_store%paw(iatom)) THEN
628  rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
629  rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
630  mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
631  mixing_store%cpc_h_in(iatom, 2)%r_coef
632  mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
633  mixing_store%cpc_s_in(iatom, 2)%r_coef
634  END IF
635  END DO
636  END IF
637  END IF
638 
639  DO ispin = 1, nspin
640  res_rho = cmplx(0.0_dp, 0.0_dp, kind=dp)
641  DO ig = 1, ng
642  res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
643  END DO
644 
645  IF (only_kerker) THEN
646  DO ig = 1, ng
647  mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
648  f_mix = alpha*mixing_store%kerker_factor(ig)
649  rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
650  mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
651  mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
652  END DO
653 
654  IF (mixing_store%gmix_p) THEN
655  IF (gapw) THEN
656  DO iatom = 1, natom
657  IF (mixing_store%paw(iatom)) THEN
658  mixing_store%cpc_h_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
659  mixing_store%cpc_h_in(iatom, ispin)%r_coef
660  mixing_store%cpc_s_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
661  mixing_store%cpc_s_in(iatom, ispin)%r_coef
662 
663  rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
664  mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha)
665  rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
666  mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha)
667 
668  mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
669  mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
670  mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
671  mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
672  END IF
673  END DO
674  END IF
675  END IF
676  ELSE
677 
678  ALLOCATE (a(nb, nb))
679  a = 0.0_dp
680  ALLOCATE (b(nb, nb))
681  b = 0.0_dp
682  ALLOCATE (c(nb))
683  c = 0.0_dp
684  ALLOCATE (g(nb))
685  g = 0.0_dp
686 
687  delta_norm = 0.0_dp
688  res_norm = 0.0_dp
689  rho_norm = 0.0_dp
690  DO ig = 1, ng
691  mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
692  mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
693  res_norm = res_norm + &
694  REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + &
695  aimag(res_rho(ig))*aimag(res_rho(ig))
696  delta_norm = delta_norm + &
697  REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)* &
698  REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
699  aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))* &
700  aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))
701  rho_norm = rho_norm + &
702  REAL(rho_g(ispin)%array(ig), dp)*REAL(rho_g(ispin)%array(ig), dp) + &
703  aimag(rho_g(ispin)%array(ig))*aimag(rho_g(ispin)%array(ig))
704  END DO
705  DO ig = 1, ng
706  mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
707  mixing_store%rhoin(ispin)%cc(ig) - &
708  mixing_store%rhoin_old(ispin)%cc(ig)
709  END DO
710  CALL para_env%sum(delta_norm)
711  delta_norm = sqrt(delta_norm)
712  CALL para_env%sum(res_norm)
713  res_norm = sqrt(res_norm)
714  CALL para_env%sum(rho_norm)
715  rho_norm = sqrt(rho_norm)
716 
717  IF (res_norm > 1.e-14_dp) THEN
718  mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm
719  mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm
720 
721  IF (mixing_store%gmix_p .AND. gapw) THEN
722  DO iatom = 1, natom
723  IF (mixing_store%paw(iatom)) THEN
724  n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
725  n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
726  DO i = 1, n2
727  DO j = 1, n1
728  mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
729  (mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
730  mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
731  dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
732  mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
733  mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
734  mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
735  mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
736 
737  mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
738  (mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
739  mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
740  dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
741  mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
742  mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
743  mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
744  mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
745 
746  mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
747  alpha*dcpc_h_res + &
748  mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
749  mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
750  alpha*dcpc_s_res + &
751  mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
752  END DO
753  END DO
754  END IF
755  END DO
756  END IF
757 
758  a(:, :) = 0.0_dp
759  DO ig = 1, ng
760  f_mix = alpha*mixing_store%kerker_factor(ig)
761  mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
762  f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
763  mixing_store%drho_buffer(ib, ispin)%cc(ig)
764  END DO
765  DO jb = 1, nb
766  DO kb = jb, nb
767  DO ig = 1, mixing_store%ig_max !ng
768  a(kb, jb) = a(kb, jb) + mixing_store%p_metric(ig)*( &
769  REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)* &
770  REAL(mixing_store%res_buffer(kb, ispin)%cc(ig), dp) + &
771  aimag(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
772  aimag(mixing_store%res_buffer(kb, ispin)%cc(ig)))
773  END DO
774  a(jb, kb) = a(kb, jb)
775  END DO
776  END DO
777  CALL para_env%sum(a)
778 
779  c = 0.0_dp
780  DO jb = 1, nb
781  a(jb, jb) = w0 + a(jb, jb)
782  DO ig = 1, mixing_store%ig_max !ng
783  c(jb) = c(jb) + mixing_store%p_metric(ig)*( &
784  REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)*REAL(res_rho(ig), dp) + &
785  aimag(mixing_store%res_buffer(jb, ispin)%cc(ig))*aimag(res_rho(ig)))
786  END DO
787  END DO
788  CALL para_env%sum(c)
789  CALL invert_matrix(a, b, inv_err)
790 
791  CALL dgemv('T', nb, nb, 1.0_dp, b, nb, c, 1, 0.0_dp, g, 1)
792  ELSE
793  g = 0.0_dp
794  END IF
795 
796  DO ig = 1, ng
797  cc_mix = cmplx(0.0_dp, 0.0_dp, kind=dp)
798  DO jb = 1, nb
799  cc_mix = cc_mix - g(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
800  END DO
801  f_mix = alpha*mixing_store%kerker_factor(ig)
802 
803  IF (res_norm > 1.e-14_dp) rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
804  f_mix*res_rho(ig) + cc_mix
805  mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
806  mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
807  END DO
808 
809  IF (mixing_store%gmix_p) THEN
810  IF (gapw) THEN
811  DO iatom = 1, natom
812  IF (mixing_store%paw(iatom)) THEN
813  n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
814  n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
815  DO i = 1, n2
816  DO j = 1, n1
817  valh = 0.0_dp
818  vals = 0.0_dp
819  DO jb = 1, nb
820  valh = valh - g(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
821  vals = vals - g(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
822  END DO
823  IF (res_norm > 1.e-14_dp) THEN
824  rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
825  alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
826  mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + valh
827  rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
828  alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
829  mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + vals
830  END IF
831  END DO
832  END DO
833 
834  mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
835  mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
836  mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
837  mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
838  END IF
839  END DO
840  END IF
841  END IF
842 
843  DEALLOCATE (a, b, c, g)
844  END IF
845 
846  END DO ! ispin
847  IF (nspin == 2) THEN
848  CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
849  DO ig = 1, ng
850  mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
851  END DO
852  IF (gapw .AND. mixing_store%gmix_p) THEN
853  DO iatom = 1, natom
854  IF (mixing_store%paw(iatom)) THEN
855  rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
856  rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
857  mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
858  mixing_store%cpc_h_in(iatom, 2)%r_coef
859  mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
860  mixing_store%cpc_s_in(iatom, 2)%r_coef
861  END IF
862  END DO
863  END IF
864 
865  END IF
866 
867  DEALLOCATE (res_rho)
868 
869  CALL timestop(handle)
870 
871  END SUBROUTINE broyden_mixing
872 
873 ! **************************************************************************************************
874 !> \brief ...
875 !> \param mixing_store ...
876 !> \param rho ...
877 !> \param para_env ...
878 ! **************************************************************************************************
879  SUBROUTINE broyden_mixing_new(mixing_store, rho, para_env)
880 
881  TYPE(mixing_storage_type), POINTER :: mixing_store
882  TYPE(qs_rho_type), POINTER :: rho
883  TYPE(mp_para_env_type), POINTER :: para_env
884 
885  CHARACTER(len=*), PARAMETER :: routinen = 'broyden_mixing_new'
886 
887  COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: delta_res_p, res_rho, res_rho_p, tmp
888  INTEGER :: handle, ib, ibb, ig, info, ispin, jb, &
889  kb, kkb, lwork, nb, nbuffer, ng, nspin
890  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
891  LOGICAL :: only_kerker, skip_bq
892  REAL(dp) :: alpha, beta, delta, delta_p, delta_rhog, delta_rhog_p, f_mix, imp, imp_j, &
893  inv_err, norm, norm_ig, rep, rep_j, sqt_uvol, sqt_vol, uvol, vol, wc, wmax
894  REAL(dp), ALLOCATABLE, DIMENSION(:) :: aval, work_dgesdd
895  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, au, av, b, bq, work
896  REAL(dp), DIMENSION(:), POINTER :: p_metric
897  REAL(dp), DIMENSION(:, :), POINTER :: fmat, weight
898  TYPE(cp_1d_z_p_type), DIMENSION(:), POINTER :: tmp_z
899  TYPE(cp_1d_z_p_type), DIMENSION(:, :), POINTER :: delta_res, u_vec, z_vec
900  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
901 
902  cpassert(ASSOCIATED(mixing_store%rhoin_buffer))
903 
904  CALL timeset(routinen, handle)
905 
906  NULLIFY (delta_res, u_vec, z_vec)
907  NULLIFY (fmat, rho_g)
908  CALL qs_rho_get(rho, rho_g=rho_g)
909 
910  nspin = SIZE(rho_g, 1)
911  ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
912  vol = rho_g(1)%pw_grid%vol
913  sqt_vol = sqrt(vol)
914  uvol = 1.0_dp/vol
915  sqt_uvol = sqrt(uvol)
916  alpha = mixing_store%alpha
917  beta = mixing_store%beta
918 
919  wc = mixing_store%wc
920  wmax = mixing_store%wmax
921  nbuffer = mixing_store%nbuffer
922 
923  mixing_store%ncall = mixing_store%ncall + 1
924  IF (mixing_store%ncall == 1) THEN
925  only_kerker = .true.
926  ELSE
927  only_kerker = .false.
928  nb = min(mixing_store%ncall - 1, nbuffer)
929  ib = modulo(mixing_store%ncall - 2, nbuffer) + 1
930 
931  ALLOCATE (a(nb, nb))
932  a = 0.0_dp
933  ALLOCATE (b(nb, nb))
934  b = 0.0_dp
935  ALLOCATE (bq(nb, nb))
936  bq = 0.0_dp
937 
938  ALLOCATE (tmp(ng))
939  ALLOCATE (delta_res_p(ng))
940  END IF
941 
942  ibb = 0
943 
944  ALLOCATE (res_rho(ng))
945  ALLOCATE (res_rho_p(ng))
946 
947  p_metric => mixing_store%p_metric
948  weight => mixing_store%weight
949  cpassert(ASSOCIATED(mixing_store%delta_res))
950  delta_res => mixing_store%delta_res
951  cpassert(ASSOCIATED(mixing_store%u_vec))
952  u_vec => mixing_store%u_vec
953  cpassert(ASSOCIATED(mixing_store%z_vec))
954  z_vec => mixing_store%z_vec
955 
956  delta_rhog = 0.0_dp
957  delta_rhog_p = 0.0_dp
958 
959  DO ispin = 1, nspin
960 
961  fmat => mixing_store%fmat(:, :, ispin)
962 
963  delta = 0.0_dp
964  delta_p = 0.0_dp
965  ! Residual at this step R_i(G) (rho_out(G)-rho_in(G))
966  ! Residual multiplied by the metrics RP_i(G) = (rho_out(G)-rho_in(G)) * P(G)
967  ! Delta is the norm of the residual, measures how far we are from convergence
968  DO ig = 1, ng
969  res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
970  res_rho_p(ig) = res_rho(ig)*p_metric(ig) !*sqt_uvol
971  norm_ig = real(res_rho(ig), dp)*real(res_rho(ig), dp) + aimag(res_rho(ig))*aimag(res_rho(ig))
972  delta = delta + norm_ig
973  delta_p = delta_p + norm_ig*p_metric(ig) !*p_metric(ig)
974  END DO
975  CALL para_env%sum(delta)
976  delta = sqrt(delta)
977  CALL para_env%sum(delta_p)
978  delta_p = sqrt(delta_p)
979  delta_rhog = delta_rhog + delta
980  delta_rhog_p = delta_rhog_p + delta_p
981 
982  weight(ib, ispin) = 1.0_dp ! wc
983  IF (wc < 0.0_dp) weight(ib, ispin) = 0.01_dp*abs(wc)/(delta_p*delta_p)
984  IF (weight(ib, ispin) == 0.0_dp) weight(ib, ispin) = 100.0_dp
985  IF (weight(ib, ispin) < 1.0_dp) weight(ib, ispin) = 1.0_dp
986 
987  IF (only_kerker) THEN
988  ! Simple Kerker damping : linear mixing rho(G) = rho_in(G) - alpha k(G)*(rho_out(G)-rho_in(G))
989  DO ig = 1, ng
990  f_mix = alpha*mixing_store%kerker_factor(ig)
991  rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
992  mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
993  mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
994  mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
995  END DO
996  ELSE
997  norm = 0.0_dp
998  ! Difference of residuals DR_{i-1)} (G) = R_i(G) - R_{i-1}(G)
999  DO ig = 1, ng
1000  delta_res(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
1001  delta_res_p(ig) = p_metric(ig)*delta_res(ib, ispin)%cc(ig)
1002  norm_ig = real(delta_res(ib, ispin)%cc(ig), dp)*real(delta_res_p(ig), dp) + &
1003  aimag(delta_res(ib, ispin)%cc(ig))*aimag(delta_res_p(ig))
1004  norm = norm + norm_ig
1005  END DO
1006  CALL para_env%sum(norm)
1007  norm = 1._dp/sqrt(norm)
1008  delta_res(ib, ispin)%cc(:) = delta_res(ib, ispin)%cc(:)*norm
1009  delta_res_p(:) = delta_res_p(:)*norm
1010 
1011  IF (para_env%is_source()) WRITE (*, *) 'norm ', norm
1012  ! Vector U_{i-1}(G) = Drho_{i-1} + k(G) * DR_{i-1}(G)
1013  DO ig = 1, ng
1014  tmp(ig) = (mixing_store%rhoin(ispin)%cc(ig) - &
1015  mixing_store%rhoin_old(ispin)%cc(ig))*norm
1016  u_vec(ib, ispin)%cc(ig) = (tmp(ig) + &
1017  mixing_store%kerker_factor(ig)*delta_res(ib, ispin)%cc(ig))
1018  END DO
1019 
1020  DO jb = 1, nb
1021  fmat(jb, ib) = 0.0_dp
1022 
1023  DO ig = 1, ng
1024  rep_j = real(delta_res(jb, ispin)%cc(ig), dp)
1025  imp_j = aimag(delta_res(jb, ispin)%cc(ig))
1026  ! < DR_{j} | DR_{i-1} >
1027  rep = real(delta_res_p(ig), dp)
1028  imp = aimag(delta_res_p(ig))
1029  fmat(jb, ib) = fmat(jb, ib) + rep_j*rep + imp_j*imp
1030  END DO
1031 
1032  END DO
1033  CALL para_env%sum(fmat(1:nb, ib))
1034 
1035  fmat(ib, ib) = 1.0_dp
1036 
1037  DO jb = 1, nb
1038  fmat(ib, jb) = fmat(jb, ib)
1039  END DO
1040 
1041  DO jb = 1, nb
1042  a(jb, jb) = weight(jb, ispin)*weight(jb, ispin)*fmat(jb, jb)
1043  DO kb = 1, jb - 1
1044  a(jb, kb) = weight(jb, ispin)*weight(kb, ispin)*fmat(jb, kb)
1045  a(kb, jb) = weight(jb, ispin)*weight(kb, ispin)*fmat(kb, jb)
1046  END DO
1047  END DO
1048  IF (.true.) THEN
1049  b = 0.0_dp
1050  CALL invert_matrix(a, b, inv_err)
1051  ELSE
1052  b = 0.0_dp
1053  ALLOCATE (work(ib - 1, ib - 1), aval(ib - 1), av(ib - 1, ib - 1), au(ib - 1, ib - 1))
1054  work(:, :) = a
1055  ALLOCATE (iwork(8*(ib - 1)), work_dgesdd(1))
1056  lwork = -1
1057  CALL dgesdd('S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info)
1058  lwork = int(work_dgesdd(1))
1059  DEALLOCATE (work_dgesdd); ALLOCATE (work_dgesdd(lwork))
1060  CALL dgesdd('S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info)
1061  ! construct the inverse
1062  DO kb = 1, ib - 1
1063  ! invert SV
1064  IF (aval(kb) < 1.e-6_dp) THEN
1065  aval(kb) = 0.0_dp
1066  ELSE
1067  aval(kb) = 1.0_dp/aval(kb)
1068  END IF
1069  av(kb, :) = av(kb, :)*aval(kb)
1070  END DO
1071  CALL dgemm('T', 'T', ib - 1, ib - 1, ib - 1, 1.0_dp, av, ib - 1, au, ib - 1, 0.0_dp, b, ib - 1)
1072  DEALLOCATE (iwork, aval, au, av, work_dgesdd, work)
1073  END IF
1074 
1075  bq = 0.0_dp
1076  ! Broyden second method requires also bq (NYI)
1077  skip_bq = .true.
1078  DO jb = 1, nb
1079  DO kb = 1, nb
1080  bq(jb, kb) = 0.0_dp
1081  DO kkb = 1, nb
1082  bq(jb, kb) = bq(jb, kb) - weight(jb, ispin)*weight(kkb, ispin)*b(jb, kkb)*fmat(kkb, kb)
1083  END DO
1084  END DO
1085  bq(jb, jb) = 1.0_dp + bq(jb, jb)
1086  END DO
1087 
1088  IF (.NOT. skip_bq) THEN
1089  ! in this case the old z_vec is needed
1090  ! a temporary array is needed to store the new one
1091  ALLOCATE (tmp_z(nb))
1092  DO jb = 1, nb
1093  ALLOCATE (tmp_z(jb)%cc(ng))
1094  END DO
1095  END IF
1096  DO jb = 1, nb
1097  tmp(:) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1098  IF (.NOT. skip_bq) THEN
1099  ! sum_{kb} bq(jb,kb) * z_vec_{kb,iter-2}
1100  ! added only for small weights
1101  DO kb = 1, nb
1102  IF (weight(kb, ispin) >= (10.0_dp*wmax)) cycle
1103  DO ig = 1, ng
1104  tmp(ig) = tmp(ig) + bq(jb, kb)*z_vec(kb, ispin)%cc(ig)
1105  END DO
1106  END DO ! kb
1107  END IF
1108 
1109  ! sum_{kb} w(jb)*w(kb)*b(jb,kb) * u_vec_{kb}
1110  DO kb = 1, ib - 1
1111  DO ig = 1, ng
1112  tmp(ig) = tmp(ig) + weight(kb, ispin)*weight(jb, ispin)*b(jb, kb)*u_vec(kb, ispin)%cc(ig)
1113  END DO
1114  END DO
1115 
1116  ! store the new z_vec(jb)
1117  IF (skip_bq .OR. (weight(jb, ispin) >= (10._dp*wmax))) THEN
1118  z_vec(jb, ispin)%cc(:) = tmp(:)
1119  ELSE
1120  ! temporary array: old z_vec may still be needed
1121  tmp_z(jb)%cc(:) = tmp(:)
1122  END IF
1123  END DO !jb
1124 
1125  IF (.NOT. skip_bq) THEN
1126  DO jb = 1, ib - 1
1127  IF (weight(jb, ispin) < (10._dp*wmax)) z_vec(jb, ispin)%cc(:) = tmp_z(jb)%cc(:)
1128  DEALLOCATE (tmp_z(jb)%cc)
1129  END DO
1130  DEALLOCATE (tmp_z)
1131  END IF
1132 
1133  ! Overwrite the density i reciprocal space
1134  rho_g(ispin)%array(:) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1135  DO jb = 1, ib - 1
1136  norm = 0.0_dp
1137  DO ig = 1, ng
1138  rep_j = real(delta_res(jb, ispin)%cc(ig), dp)
1139  imp_j = aimag(delta_res(jb, ispin)%cc(ig))
1140  rep = real(res_rho_p(ig), dp)
1141  imp = aimag(res_rho_p(ig))
1142  norm = norm + rep_j*rep + imp_j*imp
1143  END DO
1144  CALL para_env%sum(norm)
1145  ! Subtract |Z_jb)><DR_jb|P|R_{iter}>
1146  DO ig = 1, ng
1147  rho_g(ispin)%array(ig) = rho_g(ispin)%array(ig) - norm*z_vec(jb, ispin)%cc(ig)*sqt_vol
1148 
1149  END DO
1150  END DO
1151 
1152  DO ig = 1, ng
1153  f_mix = alpha*mixing_store%kerker_factor(ig)
1154  rho_g(ispin)%array(ig) = rho_g(ispin)%array(ig) + &
1155  mixing_store%rhoin_buffer(ib, ispin)%cc(ig) + f_mix*res_rho(ig)
1156  mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1157  END DO
1158 
1159  mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
1160  mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
1161  mixing_store%last_res(ispin)%cc(:) = res_rho(:)
1162  END IF ! ib
1163 
1164  END DO ! ispin
1165  IF (.NOT. only_kerker) THEN
1166  DEALLOCATE (a, b, bq)
1167  DEALLOCATE (delta_res_p, tmp)
1168  END IF
1169  DEALLOCATE (res_rho, res_rho_p)
1170 
1171  CALL timestop(handle)
1172 
1173  END SUBROUTINE broyden_mixing_new
1174 
1175 ! **************************************************************************************************
1176 !> \brief Multisecant scheme to perform charge density Mixing
1177 !> as preconditioner we use the Kerker damping factor
1178 !> The mixing is applied directly on the density in reciprocal space,
1179 !> therefore it affects the potentials
1180 !> on the grid but not the density matrix
1181 !> \param mixing_store ...
1182 !> \param rho ...
1183 !> \param para_env ...
1184 !> \par History
1185 !> 05.2009
1186 !> \author MI
1187 ! **************************************************************************************************
1188  SUBROUTINE multisecant_mixing(mixing_store, rho, para_env)
1189 
1190  TYPE(mixing_storage_type), POINTER :: mixing_store
1191  TYPE(qs_rho_type), POINTER :: rho
1192  TYPE(mp_para_env_type), POINTER :: para_env
1193 
1194  CHARACTER(len=*), PARAMETER :: routinen = 'multisecant_mixing'
1195  COMPLEX(KIND=dp), PARAMETER :: cmone = (-1.0_dp, 0.0_dp), &
1196  cone = (1.0_dp, 0.0_dp), &
1197  czero = (0.0_dp, 0.0_dp)
1198 
1199  COMPLEX(dp) :: saa, yaa
1200  COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: gn_global, pgn, res_matrix_global, &
1201  tmp_vec, ugn
1202  COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, res_matrix, sa, step_matrix, ya
1203  COMPLEX(dp), DIMENSION(:), POINTER :: gn
1204  INTEGER :: handle, ib, ib_next, ib_prev, ig, &
1205  ig_global, iig, ispin, jb, kb, nb, &
1206  nbuffer, ng, ng_global, nspin
1207  LOGICAL :: use_zgemm, use_zgemm_rev
1208  REAL(dp) :: alpha, f_mix, gn_norm, gn_norm_old, inv_err, n_low, n_up, pgn_norm, prec, &
1209  r_step, reg_par, sigma_max, sigma_tilde, step_size
1210  REAL(dp), ALLOCATABLE, DIMENSION(:) :: norm_res, norm_res_low, norm_res_up
1211  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b_matrix, binv_matrix
1212  REAL(dp), DIMENSION(:), POINTER :: g2
1213  REAL(dp), SAVE :: sigma_old = 1.0_dp
1214  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1215 
1216  cpassert(ASSOCIATED(mixing_store))
1217 
1218  CALL timeset(routinen, handle)
1219 
1220  NULLIFY (gn, rho_g)
1221 
1222  use_zgemm = .false.
1223  use_zgemm_rev = .true.
1224 
1225  ! prepare the parameters
1226  CALL qs_rho_get(rho, rho_g=rho_g)
1227 
1228  nspin = SIZE(rho_g, 1)
1229  ! not implemented for large grids.
1230  cpassert(rho_g(1)%pw_grid%ngpts < huge(ng_global))
1231  ng_global = int(rho_g(1)%pw_grid%ngpts)
1232  ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
1233  alpha = mixing_store%alpha
1234 
1235  sigma_max = mixing_store%sigma_max
1236  reg_par = mixing_store%reg_par
1237  r_step = mixing_store%r_step
1238  nbuffer = mixing_store%nbuffer
1239 
1240  ! determine the step number, and multisecant iteration
1241  nb = min(mixing_store%ncall, nbuffer - 1)
1242  ib = modulo(mixing_store%ncall, nbuffer) + 1
1243  IF (mixing_store%ncall > 0) THEN
1244  ib_prev = modulo(mixing_store%ncall - 1, nbuffer) + 1
1245  ELSE
1246  ib_prev = 0
1247  END IF
1248  mixing_store%ncall = mixing_store%ncall + 1
1249  ib_next = modulo(mixing_store%ncall, nbuffer) + 1
1250 
1251  ! compute the residual gn and its norm gn_norm
1252  DO ispin = 1, nspin
1253  gn => mixing_store%res_buffer(ib, ispin)%cc
1254  gn_norm = 0.0_dp
1255  DO ig = 1, ng
1256  gn(ig) = (rho_g(ispin)%array(ig) - mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1257  gn_norm = gn_norm + &
1258  REAL(gn(ig), dp)*REAL(gn(ig), dp) + aimag(gn(ig))*aimag(gn(ig))
1259  END DO
1260  CALL para_env%sum(gn_norm)
1261  gn_norm = sqrt(gn_norm)
1262  mixing_store%norm_res_buffer(ib, ispin) = gn_norm
1263  END DO
1264 
1265  IF (nb == 0) THEN
1266  !simple mixing
1267  DO ispin = 1, nspin
1268  DO ig = 1, ng
1269  f_mix = alpha*mixing_store%kerker_factor(ig)
1270  rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(1, ispin)%cc(ig) + &
1271  f_mix*mixing_store%res_buffer(1, ispin)%cc(ig)
1272  mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1273  END DO
1274  END DO
1275  CALL timestop(handle)
1276  RETURN
1277  END IF
1278 
1279  ! allocate temporary arrays
1280  ! step_matrix S ngxnb
1281  ALLOCATE (step_matrix(ng, nb))
1282  ! res_matrix Y ngxnb
1283  ALLOCATE (res_matrix(ng, nb))
1284  ! matrix A nbxnb
1285  ALLOCATE (a_matrix(nb, ng_global))
1286  ! PSI nb vector of norms
1287  ALLOCATE (norm_res(nb))
1288  ALLOCATE (norm_res_low(nb))
1289  ALLOCATE (norm_res_up(nb))
1290 
1291  ! matrix B nbxnb
1292  ALLOCATE (b_matrix(nb, nb))
1293  ! matrix B_inv nbxnb
1294  ALLOCATE (binv_matrix(nb, nb))
1295 
1296  ALLOCATE (gn_global(ng_global))
1297  ALLOCATE (res_matrix_global(ng_global))
1298  IF (use_zgemm) THEN
1299  ALLOCATE (sa(ng, ng_global))
1300  ALLOCATE (ya(ng, ng_global))
1301  END IF
1302  IF (use_zgemm_rev) THEN
1303  ALLOCATE (tmp_vec(nb))
1304  END IF
1305  ALLOCATE (pgn(ng))
1306  ALLOCATE (ugn(ng))
1307 
1308  DO ispin = 1, nspin
1309  ! generate the global vector with the present residual
1310  gn => mixing_store%res_buffer(ib, ispin)%cc
1311  gn_global = cmplx(0.0_dp, 0.0_dp, kind=dp)
1312  DO ig = 1, ng
1313  ig_global = mixing_store%ig_global_index(ig)
1314  gn_global(ig_global) = gn(ig)
1315  END DO
1316  CALL para_env%sum(gn_global)
1317 
1318  ! compute steps (matrix S) and residual differences (matrix Y)
1319  ! with respect to the present
1320  ! step and the present residual (use stored rho_in and res_buffer)
1321 
1322  ! These quantities are pre-conditioned by means of Kerker factor multipliccation
1323  g2 => rho_g(1)%pw_grid%gsq
1324 
1325  DO jb = 1, nb + 1
1326  IF (jb < ib) THEN
1327  kb = jb
1328  ELSEIF (jb > ib) THEN
1329  kb = jb - 1
1330  ELSE
1331  cycle
1332  END IF
1333  norm_res(kb) = 0.0_dp
1334  norm_res_low(kb) = 0.0_dp
1335  norm_res_up(kb) = 0.0_dp
1336  DO ig = 1, ng
1337 
1338  prec = 1.0_dp
1339 
1340  step_matrix(ig, kb) = prec*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) - &
1341  mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1342  res_matrix(ig, kb) = (mixing_store%res_buffer(jb, ispin)%cc(ig) - &
1343  mixing_store%res_buffer(ib, ispin)%cc(ig))
1344  norm_res(kb) = norm_res(kb) + real(res_matrix(ig, kb), dp)*real(res_matrix(ig, kb), dp) + &
1345  aimag(res_matrix(ig, kb))*aimag(res_matrix(ig, kb))
1346  IF (g2(ig) < 4.0_dp) THEN
1347  norm_res_low(kb) = norm_res_low(kb) + &
1348  REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1349  aimag(res_matrix(ig, kb))*aimag(res_matrix(ig, kb))
1350  ELSE
1351  norm_res_up(kb) = norm_res_up(kb) + &
1352  REAL(res_matrix(ig, kb), dp)*REAL(res_matrix(ig, kb), dp) + &
1353  aimag(res_matrix(ig, kb))*aimag(res_matrix(ig, kb))
1354  END IF
1355  res_matrix(ig, kb) = prec*res_matrix(ig, kb)
1356  END DO
1357  END DO !jb
1358 
1359  ! normalize each column of S and Y => Snorm Ynorm
1360  CALL para_env%sum(norm_res)
1361  CALL para_env%sum(norm_res_up)
1362  CALL para_env%sum(norm_res_low)
1363  norm_res(1:nb) = 1.0_dp/sqrt(norm_res(1:nb))
1364  n_low = 0.0_dp
1365  n_up = 0.0_dp
1366  DO jb = 1, nb
1367  n_low = n_low + norm_res_low(jb)/norm_res(jb)
1368  n_up = n_up + norm_res_up(jb)/norm_res(jb)
1369  END DO
1370  DO ig = 1, ng
1371  IF (g2(ig) > 4.0_dp) THEN
1372  step_matrix(ig, 1:nb) = step_matrix(ig, 1:nb)*sqrt(n_low/n_up)
1373  res_matrix(ig, 1:nb) = res_matrix(ig, 1:nb)*sqrt(n_low/n_up)
1374  END IF
1375  END DO
1376  DO kb = 1, nb
1377  step_matrix(1:ng, kb) = step_matrix(1:ng, kb)*norm_res(kb)
1378  res_matrix(1:ng, kb) = res_matrix(1:ng, kb)*norm_res(kb)
1379  END DO
1380 
1381  ! compute A as [(Ynorm^t Ynorm) + (alpha I)]^(-1) Ynorm^t
1382  ! compute B
1383  DO jb = 1, nb
1384  DO kb = 1, nb
1385  b_matrix(kb, jb) = 0.0_dp
1386  DO ig = 1, ng
1387  ! it is assumed that summing over all G vector gives a real, because
1388  ! y(-G,kb) = (y(G,kb))*
1389  b_matrix(kb, jb) = b_matrix(kb, jb) + real(res_matrix(ig, kb)*res_matrix(ig, jb), dp)
1390  END DO
1391  END DO
1392  END DO
1393 
1394  CALL para_env%sum(b_matrix)
1395  DO jb = 1, nb
1396  b_matrix(jb, jb) = b_matrix(jb, jb) + reg_par
1397  END DO
1398  ! invert B
1399  CALL invert_matrix(b_matrix, binv_matrix, inv_err)
1400 
1401  ! A = Binv Ynorm^t
1402  a_matrix = cmplx(0.0_dp, 0.0_dp, kind=dp)
1403  DO kb = 1, nb
1404  res_matrix_global = cmplx(0.0_dp, 0.0_dp, kind=dp)
1405  DO ig = 1, ng
1406  ig_global = mixing_store%ig_global_index(ig)
1407  res_matrix_global(ig_global) = res_matrix(ig, kb)
1408  END DO
1409  CALL para_env%sum(res_matrix_global)
1410 
1411  DO jb = 1, nb
1412  DO ig = 1, ng
1413  ig_global = mixing_store%ig_global_index(ig)
1414  a_matrix(jb, ig_global) = a_matrix(jb, ig_global) + &
1415  binv_matrix(jb, kb)*res_matrix_global(ig_global)
1416  END DO
1417  END DO
1418  END DO
1419  CALL para_env%sum(a_matrix)
1420 
1421  ! compute the two components of gn that will be used to update rho
1422  gn => mixing_store%res_buffer(ib, ispin)%cc
1423  pgn_norm = 0.0_dp
1424 
1425  IF (use_zgemm) THEN
1426 
1427  CALL zgemm("N", "N", ng, ng_global, nb, cmone, step_matrix(1, 1), ng, &
1428  a_matrix(1, 1), nb, czero, sa(1, 1), ng)
1429  CALL zgemm("N", "N", ng, ng_global, nb, cmone, res_matrix(1, 1), ng, &
1430  a_matrix(1, 1), nb, czero, ya(1, 1), ng)
1431  DO ig = 1, ng
1432  ig_global = mixing_store%ig_global_index(ig)
1433  ya(ig, ig_global) = ya(ig, ig_global) + cmplx(1.0_dp, 0.0_dp, kind=dp)
1434  END DO
1435 
1436  CALL zgemv("N", ng, ng_global, cone, sa(1, 1), &
1437  ng, gn_global(1), 1, czero, pgn(1), 1)
1438  CALL zgemv("N", ng, ng_global, cone, ya(1, 1), &
1439  ng, gn_global(1), 1, czero, ugn(1), 1)
1440 
1441  DO ig = 1, ng
1442  pgn_norm = pgn_norm + real(pgn(ig), dp)*real(pgn(ig), dp) + &
1443  aimag(pgn(ig))*aimag(pgn(ig))
1444  END DO
1445  CALL para_env%sum(pgn_norm)
1446  ELSEIF (use_zgemm_rev) THEN
1447 
1448  CALL zgemv("N", nb, ng_global, cone, a_matrix(1, 1), &
1449  nb, gn_global(1), 1, czero, tmp_vec(1), 1)
1450 
1451  CALL zgemv("N", ng, nb, cmone, step_matrix(1, 1), ng, &
1452  tmp_vec(1), 1, czero, pgn(1), 1)
1453 
1454  CALL zgemv("N", ng, nb, cmone, res_matrix(1, 1), ng, &
1455  tmp_vec(1), 1, czero, ugn(1), 1)
1456 
1457  DO ig = 1, ng
1458  pgn_norm = pgn_norm + real(pgn(ig), dp)*real(pgn(ig), dp) + &
1459  aimag(pgn(ig))*aimag(pgn(ig))
1460  ugn(ig) = ugn(ig) + gn(ig)
1461  END DO
1462  CALL para_env%sum(pgn_norm)
1463 
1464  ELSE
1465  DO ig = 1, ng
1466  pgn(ig) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1467  ugn(ig) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1468  ig_global = mixing_store%ig_global_index(ig)
1469  DO iig = 1, ng_global
1470  saa = cmplx(0.0_dp, 0.0_dp, kind=dp)
1471  yaa = cmplx(0.0_dp, 0.0_dp, kind=dp)
1472 
1473  IF (ig_global == iig) yaa = cmplx(1.0_dp, 0.0_dp, kind=dp)
1474 
1475  DO jb = 1, nb
1476  saa = saa - step_matrix(ig, jb)*a_matrix(jb, iig)
1477  yaa = yaa - res_matrix(ig, jb)*a_matrix(jb, iig)
1478  END DO
1479  pgn(ig) = pgn(ig) + saa*gn_global(iig)
1480  ugn(ig) = ugn(ig) + yaa*gn_global(iig)
1481  END DO
1482  END DO
1483  DO ig = 1, ng
1484  pgn_norm = pgn_norm + real(pgn(ig), dp)*real(pgn(ig), dp) + &
1485  aimag(pgn(ig))*aimag(pgn(ig))
1486  END DO
1487  CALL para_env%sum(pgn_norm)
1488  END IF
1489 
1490  gn_norm = mixing_store%norm_res_buffer(ib, ispin)
1491  gn_norm_old = mixing_store%norm_res_buffer(ib_prev, ispin)
1492  IF (ib_prev /= 0) THEN
1493  sigma_tilde = sigma_old*max(0.5_dp, min(2.0_dp, gn_norm_old/gn_norm))
1494  ELSE
1495  sigma_tilde = 0.5_dp
1496  END IF
1497  sigma_tilde = 0.1_dp
1498  ! Step size for the unpredicted component
1499  step_size = min(sigma_tilde, r_step*pgn_norm/gn_norm, sigma_max)
1500  sigma_old = step_size
1501 
1502  ! update the density
1503  DO ig = 1, ng
1504  prec = mixing_store%kerker_factor(ig)
1505  rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(ib, ispin)%cc(ig) &
1506  - prec*step_size*ugn(ig) + prec*pgn(ig) ! - 0.1_dp * prec* gn(ig)
1507  mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1508  END DO
1509 
1510  END DO ! ispin
1511 
1512  ! Deallocate temporary arrays
1513  DEALLOCATE (step_matrix, res_matrix)
1514  DEALLOCATE (norm_res)
1515  DEALLOCATE (norm_res_up)
1516  DEALLOCATE (norm_res_low)
1517  DEALLOCATE (a_matrix, b_matrix, binv_matrix)
1518  DEALLOCATE (ugn, pgn)
1519  IF (use_zgemm) THEN
1520  DEALLOCATE (sa, ya)
1521  END IF
1522  IF (use_zgemm_rev) THEN
1523  DEALLOCATE (tmp_vec)
1524  END IF
1525  DEALLOCATE (gn_global, res_matrix_global)
1526 
1527  CALL timestop(handle)
1528 
1529  END SUBROUTINE multisecant_mixing
1530 
1531 END MODULE qs_gspace_mixing
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
Interface to the message passing library MPI.
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
module that contains the definitions of the scf types
integer, parameter, public broyden_mixing_nr
integer, parameter, public multisecant_mixing_nr
integer, parameter, public pulay_mixing_nr
integer, parameter, public gspace_mixing_nr
integer, parameter, public broyden_mixing_new_nr
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
Driver for the g-space mixing, calls the proper routine given the requested method.
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229