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