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