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