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