Line data Source code
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 : ! **************************************************************************************************
9 : MODULE 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 :
46 : CONTAINS
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 1880 : 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 1880 : 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 1880 : 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 1880 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
84 1880 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
85 :
86 1880 : CALL timeset(routineN, handle)
87 :
88 1880 : CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env)
89 1880 : IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
90 : dft_control%qs_control%semi_empirical)) THEN
91 :
92 1336 : CPASSERT(ASSOCIATED(mixing_store))
93 1336 : NULLIFY (auxbas_pw_pool, rho_atom, rho_g, rho_r, tot_rho_r)
94 1336 : CALL qs_rho_get(rho, rho_g=rho_g, rho_r=rho_r, tot_rho_r=tot_rho_r)
95 :
96 1336 : nspin = SIZE(rho_g, 1)
97 1336 : nimg = dft_control%nimages
98 1336 : ng = SIZE(rho_g(1)%pw_grid%gsq)
99 1336 : CPASSERT(ng == SIZE(mixing_store%rhoin(1)%cc))
100 :
101 1336 : alpha = mixing_store%alpha
102 1336 : gapw = dft_control%qs_control%gapw
103 :
104 1336 : IF (nspin == 2) THEN
105 114 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
106 114 : CALL auxbas_pw_pool%create_pw(rho_tmp)
107 114 : CALL pw_zero(rho_tmp)
108 114 : CALL pw_copy(rho_g(1), rho_tmp)
109 114 : CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
110 114 : CALL pw_axpy(rho_tmp, rho_g(2), -1.0_dp)
111 114 : CALL pw_scale(rho_g(2), -1.0_dp)
112 : END IF
113 :
114 1336 : IF (iter_count + 1 <= mixing_store%nskip_mixing) THEN
115 : ! skip mixing
116 8 : DO ispin = 1, nspin
117 27656 : DO ig = 1, ng
118 27652 : mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
119 : END DO
120 : END DO
121 4 : IF (mixing_store%gmix_p .AND. gapw) THEN
122 0 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
123 0 : natom = SIZE(rho_atom)
124 0 : DO ispin = 1, nspin
125 0 : DO iatom = 1, natom
126 0 : IF (mixing_store%paw(iatom)) THEN
127 0 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
128 0 : 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 4 : mixing_store%iter_method = "NoMix"
135 4 : IF (nspin == 2) THEN
136 0 : CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
137 0 : CALL pw_scale(rho_g(1), 0.5_dp)
138 0 : CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
139 0 : CALL pw_scale(rho_g(2), -1.0_dp)
140 0 : CALL auxbas_pw_pool%give_back_pw(rho_tmp)
141 : END IF
142 4 : CALL timestop(handle)
143 4 : RETURN
144 : END IF
145 :
146 1332 : IF ((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix) THEN
147 6 : CALL gmix_potential_only(qs_env, mixing_store, rho)
148 6 : mixing_store%iter_method = "Kerker"
149 1326 : ELSEIF (mixing_method == gspace_mixing_nr) THEN
150 40 : CALL gmix_potential_only(qs_env, mixing_store, rho)
151 40 : mixing_store%iter_method = "Kerker"
152 1286 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
153 184 : CALL pulay_mixing(qs_env, mixing_store, rho, para_env)
154 184 : mixing_store%iter_method = "Pulay"
155 1102 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
156 1102 : CALL broyden_mixing(qs_env, mixing_store, rho, para_env)
157 1102 : mixing_store%iter_method = "Broy."
158 0 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
159 0 : CPASSERT(.NOT. gapw)
160 0 : CALL multisecant_mixing(mixing_store, rho, para_env)
161 0 : mixing_store%iter_method = "MSec."
162 : END IF
163 :
164 1332 : IF (nspin == 2) THEN
165 114 : CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
166 114 : CALL pw_scale(rho_g(1), 0.5_dp)
167 114 : CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
168 114 : CALL pw_scale(rho_g(2), -1.0_dp)
169 114 : CALL auxbas_pw_pool%give_back_pw(rho_tmp)
170 : END IF
171 :
172 2778 : DO ispin = 1, nspin
173 1446 : CALL pw_transfer(rho_g(ispin), rho_r(ispin))
174 2778 : tot_rho_r(ispin) = pw_integrate_function(rho_r(ispin), isign=-1)
175 : END DO
176 : END IF
177 :
178 1876 : CALL timestop(handle)
179 :
180 1880 : 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 92 : 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 46 : 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 46 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
208 46 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
209 :
210 0 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
211 46 : CPASSERT(ASSOCIATED(mixing_store%kerker_factor))
212 :
213 46 : CALL timeset(routineN, handle)
214 46 : NULLIFY (cc_new, dft_control, rho_atom, rho_g)
215 :
216 46 : CALL get_qs_env(qs_env, dft_control=dft_control)
217 46 : CALL qs_rho_get(rho, rho_g=rho_g)
218 :
219 46 : nspin = SIZE(rho_g, 1)
220 46 : ng = SIZE(rho_g(1)%pw_grid%gsq)
221 :
222 46 : gapw = dft_control%qs_control%gapw
223 46 : alpha = mixing_store%alpha
224 :
225 92 : DO ispin = 1, nspin
226 46 : cc_new => rho_g(ispin)%array
227 580654 : DO ig = 1, mixing_store%ig_max ! ng
228 580608 : f_mix = mixing_store%alpha*mixing_store%kerker_factor(ig)
229 580608 : cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
230 580654 : mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
231 : END DO
232 92 : DO ig = mixing_store%ig_max + 1, ng
233 0 : f_mix = mixing_store%alpha
234 0 : cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
235 46 : mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
236 : END DO
237 :
238 : END DO
239 :
240 46 : IF (mixing_store%gmix_p .AND. gapw) THEN
241 : CALL get_qs_env(qs_env=qs_env, &
242 8 : rho_atom_set=rho_atom)
243 8 : natom = SIZE(rho_atom)
244 16 : DO ispin = 1, nspin
245 80 : DO iatom = 1, natom
246 72 : 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 7408 : 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 7408 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha)
251 7408 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
252 7408 : 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 46 : CALL timestop(handle)
259 :
260 46 : 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 184 : 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 184 : 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 184 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, ev
294 184 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b, c, c_inv, cpc_h_mix, cpc_s_mix
295 : TYPE(dft_control_type), POINTER :: dft_control
296 184 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
297 184 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
298 :
299 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
300 184 : CPASSERT(ASSOCIATED(mixing_store%rhoin_buffer))
301 184 : NULLIFY (dft_control, rho_atom, rho_g)
302 184 : CALL timeset(routineN, handle)
303 :
304 184 : CALL get_qs_env(qs_env, dft_control=dft_control)
305 184 : CALL qs_rho_get(rho, rho_g=rho_g)
306 184 : nspin = SIZE(rho_g, 1)
307 184 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
308 184 : vol = rho_g(1)%pw_grid%vol
309 :
310 184 : alpha_kerker = mixing_store%alpha
311 184 : beta = mixing_store%pulay_beta
312 184 : alpha_pulay = mixing_store%pulay_alpha
313 184 : nbuffer = mixing_store%nbuffer
314 184 : gapw = dft_control%qs_control%gapw
315 184 : IF (gapw) THEN
316 12 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
317 12 : natom = SIZE(rho_atom)
318 : END IF
319 :
320 552 : ALLOCATE (cc_mix(ng))
321 :
322 184 : IF (nspin == 2) THEN
323 14 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
324 96782 : DO ig = 1, ng
325 96782 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
326 : END DO
327 14 : IF (gapw .AND. mixing_store%gmix_p) THEN
328 0 : DO iatom = 1, natom
329 0 : IF (mixing_store%paw(iatom)) THEN
330 0 : 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 0 : 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 0 : 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 0 : mixing_store%cpc_s_in(iatom, 2)%r_coef
336 : END IF
337 : END DO
338 : END IF
339 : END IF
340 :
341 382 : DO ispin = 1, nspin
342 198 : ib = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
343 :
344 198 : mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) + 1
345 198 : nb = MIN(mixing_store%ncall_p(ispin), nbuffer)
346 198 : ibb = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
347 :
348 8641926 : mixing_store%res_buffer(ib, ispin)%cc(:) = CMPLX(0._dp, 0._dp, KIND=dp)
349 198 : norm = 0.0_dp
350 1306566 : IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
351 198 : res_norm = 0.0_dp
352 8641926 : DO ig = 1, ng
353 8641728 : f_mix = mixing_store%kerker_factor(ig)
354 : mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(rho_g(ispin)%array(ig) - &
355 8641728 : 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 8641926 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))*AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
359 : END DO
360 :
361 198 : CALL para_env%sum(res_norm)
362 198 : res_norm = SQRT(res_norm)
363 :
364 198 : IF (res_norm < 1.E-14_dp) THEN
365 0 : mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) - 1
366 : ELSE
367 198 : nb1 = nb + 1
368 198 : IF (ALLOCATED(b)) DEALLOCATE (b)
369 792 : ALLOCATE (b(nb1, nb1))
370 5202 : b = 0.0_dp
371 198 : IF (ALLOCATED(c)) DEALLOCATE (c)
372 792 : ALLOCATE (c(nb, nb))
373 3518 : c = 0.0_dp
374 198 : IF (ALLOCATED(c_inv)) DEALLOCATE (c_inv)
375 594 : ALLOCATE (c_inv(nb, nb))
376 3518 : c_inv = 0.0_dp
377 198 : IF (ALLOCATED(alpha_c)) DEALLOCATE (alpha_c)
378 594 : ALLOCATE (alpha_c(nb))
379 842 : alpha_c = 0.0_dp
380 198 : norm_c_inv = 0.0_dp
381 198 : IF (ALLOCATED(ev)) DEALLOCATE (ev)
382 594 : ALLOCATE (ev(nb1))
383 1040 : ev = 0.0_dp
384 :
385 : END IF
386 :
387 842 : DO jb = 1, nb
388 644 : mixing_store%pulay_matrix(jb, ib, ispin) = 0.0_dp
389 33413252 : DO ig = 1, ng
390 33412608 : 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 33413252 : AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig)))
396 : END DO
397 644 : CALL para_env%sum(mixing_store%pulay_matrix(jb, ib, ispin))
398 644 : mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)*vol*vol
399 842 : mixing_store%pulay_matrix(ib, jb, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)
400 : END DO
401 :
402 198 : IF (nb == 1 .OR. res_norm < 1.E-14_dp) THEN
403 1306406 : DO ig = 1, ng
404 1306368 : f_mix = alpha_kerker*mixing_store%kerker_factor(ig)
405 : cc_mix(ig) = rho_g(ispin)%array(ig) - &
406 1306368 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
407 : rho_g(ispin)%array(ig) = f_mix*cc_mix(ig) + &
408 1306368 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
409 1306406 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
410 : END DO
411 38 : IF (mixing_store%gmix_p) THEN
412 2 : IF (gapw) THEN
413 18 : DO iatom = 1, natom
414 18 : 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 1850 : 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 1850 : 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 1850 : (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 1850 : (1.0_dp - alpha_kerker)*mixing_store%cpc_s_in(iatom, ispin)%r_coef
424 :
425 926 : mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
426 926 : mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
427 :
428 1850 : mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
429 1850 : 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 3404 : b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb, ispin)
437 3404 : c(1:nb, 1:nb) = b(1:nb, 1:nb)
438 766 : b(nb1, 1:nb) = -1.0_dp
439 766 : b(1:nb, nb1) = -1.0_dp
440 160 : b(nb1, nb1) = 0.0_dp
441 :
442 160 : CALL invert_matrix(c, c_inv, inv_err, improve=.TRUE.)
443 766 : alpha_c = 0.0_dp
444 766 : DO i = 1, nb
445 3404 : DO jb = 1, nb
446 2638 : alpha_c(i) = alpha_c(i) + c_inv(jb, i)
447 3244 : norm_c_inv = norm_c_inv + c_inv(jb, i)
448 : END DO
449 : END DO
450 766 : alpha_c(1:nb) = alpha_c(1:nb)/norm_c_inv
451 7335520 : cc_mix = CMPLX(0._dp, 0._dp, KIND=dp)
452 766 : DO jb = 1, nb
453 32107006 : DO ig = 1, ng
454 : cc_mix(ig) = cc_mix(ig) + &
455 : alpha_c(jb)*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) + &
456 32106846 : mixing_store%pulay_beta*mixing_store%res_buffer(jb, ispin)%cc(ig))
457 : END DO
458 : END DO
459 7335520 : mixing_store%rhoin_buffer(ibb, ispin)%cc = CMPLX(0._dp, 0._dp, KIND=dp)
460 160 : IF (alpha_pulay > 0.0_dp) THEN
461 233290 : DO ig = 1, ng
462 233280 : f_mix = alpha_pulay*mixing_store%kerker_factor(ig)
463 : rho_g(ispin)%array(ig) = f_mix*rho_g(ispin)%array(ig) + &
464 233280 : (1.0_dp - f_mix)*cc_mix(ig)
465 233290 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
466 : END DO
467 : ELSE
468 7102230 : DO ig = 1, ng
469 7102080 : rho_g(ispin)%array(ig) = cc_mix(ig)
470 7102230 : mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
471 : END DO
472 : END IF
473 :
474 160 : IF (mixing_store%gmix_p .AND. gapw) THEN
475 10 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
476 90 : DO iatom = 1, natom
477 90 : IF (mixing_store%paw(iatom)) THEN
478 10 : n1 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 1)
479 10 : n2 = SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 2)
480 40 : ALLOCATE (cpc_h_mix(n1, n2))
481 30 : 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 9250 : 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 9250 : mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef
487 4630 : cpc_h_mix = 0.0_dp
488 4630 : cpc_s_mix = 0.0_dp
489 48 : 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 17594 : 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 17604 : 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 4630 : (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 4630 : (1.0_dp - alpha_pulay)*cpc_s_mix
501 9250 : mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
502 9250 : mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
503 10 : DEALLOCATE (cpc_h_mix)
504 10 : DEALLOCATE (cpc_s_mix)
505 : END IF
506 : END DO
507 : END IF
508 : END IF ! nb==1
509 :
510 382 : IF (res_norm > 1.E-14_dp) THEN
511 198 : DEALLOCATE (b)
512 198 : DEALLOCATE (ev)
513 198 : DEALLOCATE (c, c_inv, alpha_c)
514 : END IF
515 :
516 : END DO ! ispin
517 :
518 184 : DEALLOCATE (cc_mix)
519 :
520 184 : IF (nspin == 2) THEN
521 14 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
522 96782 : DO ig = 1, ng
523 96782 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
524 : END DO
525 14 : IF (gapw .AND. mixing_store%gmix_p) THEN
526 0 : DO iatom = 1, natom
527 0 : IF (mixing_store%paw(iatom)) THEN
528 0 : 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 0 : 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 0 : 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 0 : mixing_store%cpc_s_in(iatom, 2)%r_coef
534 : END IF
535 : END DO
536 : END IF
537 :
538 : END IF
539 :
540 184 : CALL timestop(handle)
541 :
542 552 : 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 1102 : 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 1102 : 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 1102 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g
575 1102 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
576 : TYPE(dft_control_type), POINTER :: dft_control
577 1102 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
578 1102 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
579 :
580 0 : CPASSERT(ASSOCIATED(mixing_store%res_buffer))
581 1102 : CPASSERT(ASSOCIATED(mixing_store%rhoin))
582 1102 : CPASSERT(ASSOCIATED(mixing_store%rhoin_old))
583 1102 : CPASSERT(ASSOCIATED(mixing_store%drho_buffer))
584 1102 : NULLIFY (dft_control, rho_atom, rho_g)
585 1102 : CALL timeset(routineN, handle)
586 :
587 1102 : CALL get_qs_env(qs_env, dft_control=dft_control)
588 1102 : CALL qs_rho_get(rho, rho_g=rho_g)
589 :
590 1102 : nspin = SIZE(rho_g, 1)
591 1102 : ng = SIZE(mixing_store%res_buffer(1, 1)%cc)
592 :
593 1102 : alpha = mixing_store%alpha
594 1102 : w0 = mixing_store%broy_w0
595 1102 : nbuffer = mixing_store%nbuffer
596 1102 : gapw = dft_control%qs_control%gapw
597 :
598 3306 : ALLOCATE (res_rho(ng))
599 :
600 1102 : mixing_store%ncall = mixing_store%ncall + 1
601 1102 : IF (mixing_store%ncall == 1) THEN
602 : only_kerker = .TRUE.
603 : ELSE
604 940 : only_kerker = .FALSE.
605 940 : nb = MIN(mixing_store%ncall - 1, nbuffer)
606 940 : ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1
607 : END IF
608 :
609 1102 : IF (gapw) THEN
610 238 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
611 238 : natom = SIZE(rho_atom)
612 : ELSE
613 : natom = 0
614 : END IF
615 :
616 1102 : IF (nspin == 2) THEN
617 100 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
618 9659620 : DO ig = 1, ng
619 9659620 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
620 : END DO
621 100 : IF (gapw .AND. mixing_store%gmix_p) THEN
622 90 : DO iatom = 1, natom
623 90 : IF (mixing_store%paw(iatom)) THEN
624 60560 : 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 60560 : 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 60560 : 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 60560 : mixing_store%cpc_s_in(iatom, 2)%r_coef
630 : END IF
631 : END DO
632 : END IF
633 : END IF
634 :
635 2304 : DO ispin = 1, nspin
636 100621384 : res_rho = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
637 100621384 : DO ig = 1, ng
638 100621384 : res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
639 : END DO
640 :
641 2304 : IF (only_kerker) THEN
642 10724966 : DO ig = 1, ng
643 10724774 : mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
644 10724774 : f_mix = alpha*mixing_store%kerker_factor(ig)
645 10724774 : rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
646 10724774 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
647 10724966 : mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
648 : END DO
649 :
650 192 : IF (mixing_store%gmix_p) THEN
651 24 : IF (gapw) THEN
652 216 : DO iatom = 1, natom
653 216 : 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 245780 : 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 245780 : 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 245780 : 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 245780 : mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha)
663 :
664 122972 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
665 122972 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
666 245780 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
667 245780 : 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 4040 : ALLOCATE (a(nb, nb))
675 32798 : a = 0.0_dp
676 3030 : ALLOCATE (b(nb, nb))
677 32798 : b = 0.0_dp
678 3030 : ALLOCATE (c(nb))
679 5440 : c = 0.0_dp
680 2020 : ALLOCATE (g(nb))
681 5440 : g = 0.0_dp
682 :
683 1010 : delta_norm = 0.0_dp
684 1010 : res_norm = 0.0_dp
685 1010 : rho_norm = 0.0_dp
686 89896418 : DO ig = 1, ng
687 89895408 : mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
688 89895408 : 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 89895408 : 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 89895408 : 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 89896418 : AIMAG(rho_g(ispin)%array(ig))*AIMAG(rho_g(ispin)%array(ig))
700 : END DO
701 89896418 : DO ig = 1, ng
702 : mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
703 : mixing_store%rhoin(ispin)%cc(ig) - &
704 89896418 : mixing_store%rhoin_old(ispin)%cc(ig)
705 : END DO
706 1010 : CALL para_env%sum(delta_norm)
707 1010 : delta_norm = SQRT(delta_norm)
708 1010 : CALL para_env%sum(res_norm)
709 1010 : res_norm = SQRT(res_norm)
710 1010 : CALL para_env%sum(rho_norm)
711 1010 : rho_norm = SQRT(rho_norm)
712 :
713 1010 : IF (res_norm > 1.E-14_dp) THEN
714 89896418 : mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm
715 89896418 : mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm
716 :
717 1010 : IF (mixing_store%gmix_p .AND. gapw) THEN
718 252 : DO iatom = 1, natom
719 252 : IF (mixing_store%paw(iatom)) THEN
720 28 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
721 28 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
722 616 : DO i = 1, n2
723 12964 : 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 12348 : 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 12348 : 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 12348 : 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 12348 : 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 12348 : 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 12348 : 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 12348 : 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 12936 : 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 32798 : a(:, :) = 0.0_dp
755 89896418 : DO ig = 1, ng
756 89895408 : 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 89896418 : mixing_store%drho_buffer(ib, ispin)%cc(ig)
760 : END DO
761 5440 : DO jb = 1, nb
762 21334 : DO kb = jb, nb
763 2014107494 : 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 2014107494 : AIMAG(mixing_store%res_buffer(kb, ispin)%cc(ig)))
769 : END DO
770 20324 : a(jb, kb) = a(kb, jb)
771 : END DO
772 : END DO
773 1010 : CALL para_env%sum(a)
774 :
775 5440 : C = 0.0_dp
776 5440 : DO jb = 1, nb
777 4430 : a(jb, jb) = w0 + a(jb, jb)
778 476123128 : 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 476122118 : AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))*AIMAG(res_rho(ig)))
782 : END DO
783 : END DO
784 1010 : CALL para_env%sum(c)
785 1010 : CALL invert_matrix(a, b, inv_err)
786 :
787 1010 : CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1)
788 : ELSE
789 0 : G = 0.0_dp
790 : END IF
791 :
792 89896418 : DO ig = 1, ng
793 89895408 : cc_mix = CMPLX(0.0_dp, 0.0_dp, kind=dp)
794 566013096 : DO jb = 1, nb
795 566013096 : cc_mix = cc_mix - G(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
796 : END DO
797 89895408 : f_mix = alpha*mixing_store%kerker_factor(ig)
798 :
799 89895408 : IF (res_norm > 1.E-14_dp) rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
800 89895408 : f_mix*res_rho(ig) + cc_mix
801 89895408 : mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
802 89896418 : mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
803 : END DO
804 :
805 1010 : IF (mixing_store%gmix_p) THEN
806 28 : IF (gapw) THEN
807 252 : DO iatom = 1, natom
808 252 : IF (mixing_store%paw(iatom)) THEN
809 28 : n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
810 28 : n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
811 616 : DO i = 1, n2
812 12964 : DO j = 1, n1
813 12348 : valh = 0.0_dp
814 12348 : vals = 0.0_dp
815 61740 : DO jb = 1, nb
816 49392 : valh = valh - G(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
817 61740 : vals = vals - G(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
818 : END DO
819 12936 : 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 12348 : 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 12348 : 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 12964 : mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
831 12964 : mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
832 25900 : mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
833 25900 : 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 1010 : DEALLOCATE (a, b, c, g)
840 : END IF
841 :
842 : END DO ! ispin
843 1102 : IF (nspin == 2) THEN
844 100 : CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
845 9659620 : DO ig = 1, ng
846 9659620 : mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
847 : END DO
848 100 : IF (gapw .AND. mixing_store%gmix_p) THEN
849 90 : DO iatom = 1, natom
850 90 : IF (mixing_store%paw(iatom)) THEN
851 60560 : 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 60560 : 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 60560 : 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 60560 : mixing_store%cpc_s_in(iatom, 2)%r_coef
857 : END IF
858 : END DO
859 : END IF
860 :
861 : END IF
862 :
863 1102 : DEALLOCATE (res_rho)
864 :
865 1102 : CALL timestop(handle)
866 :
867 3306 : 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 0 : 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 0 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: gn_global, pgn, res_matrix_global, &
895 0 : tmp_vec, ugn
896 0 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, res_matrix, sa, step_matrix, ya
897 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: norm_res, norm_res_low, norm_res_up
905 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: b_matrix, binv_matrix
906 0 : REAL(dp), DIMENSION(:), POINTER :: g2
907 : REAL(dp), SAVE :: sigma_old = 1.0_dp
908 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
909 :
910 0 : CPASSERT(ASSOCIATED(mixing_store))
911 :
912 0 : CALL timeset(routineN, handle)
913 :
914 0 : NULLIFY (gn, rho_g)
915 :
916 0 : use_zgemm = .FALSE.
917 0 : use_zgemm_rev = .TRUE.
918 :
919 : ! prepare the parameters
920 0 : CALL qs_rho_get(rho, rho_g=rho_g)
921 :
922 0 : nspin = SIZE(rho_g, 1)
923 : ! not implemented for large grids.
924 0 : CPASSERT(rho_g(1)%pw_grid%ngpts < HUGE(ng_global))
925 0 : ng_global = INT(rho_g(1)%pw_grid%ngpts)
926 0 : ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
927 0 : alpha = mixing_store%alpha
928 :
929 0 : sigma_max = mixing_store%sigma_max
930 0 : reg_par = mixing_store%reg_par
931 0 : r_step = mixing_store%r_step
932 0 : nbuffer = mixing_store%nbuffer
933 :
934 : ! determine the step number, and multisecant iteration
935 0 : nb = MIN(mixing_store%ncall, nbuffer - 1)
936 0 : ib = MODULO(mixing_store%ncall, nbuffer) + 1
937 0 : IF (mixing_store%ncall > 0) THEN
938 0 : ib_prev = MODULO(mixing_store%ncall - 1, nbuffer) + 1
939 : ELSE
940 : ib_prev = 0
941 : END IF
942 0 : mixing_store%ncall = mixing_store%ncall + 1
943 0 : ib_next = MODULO(mixing_store%ncall, nbuffer) + 1
944 :
945 : ! compute the residual gn and its norm gn_norm
946 0 : DO ispin = 1, nspin
947 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
948 0 : gn_norm = 0.0_dp
949 0 : DO ig = 1, ng
950 0 : gn(ig) = (rho_g(ispin)%array(ig) - mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
951 : gn_norm = gn_norm + &
952 0 : REAL(gn(ig), dp)*REAL(gn(ig), dp) + AIMAG(gn(ig))*AIMAG(gn(ig))
953 : END DO
954 0 : CALL para_env%sum(gn_norm)
955 0 : gn_norm = SQRT(gn_norm)
956 0 : mixing_store%norm_res_buffer(ib, ispin) = gn_norm
957 : END DO
958 :
959 0 : IF (nb == 0) THEN
960 : !simple mixing
961 0 : DO ispin = 1, nspin
962 0 : DO ig = 1, ng
963 0 : f_mix = alpha*mixing_store%kerker_factor(ig)
964 : rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(1, ispin)%cc(ig) + &
965 0 : f_mix*mixing_store%res_buffer(1, ispin)%cc(ig)
966 0 : mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
967 : END DO
968 : END DO
969 0 : CALL timestop(handle)
970 0 : RETURN
971 : END IF
972 :
973 : ! allocate temporary arrays
974 : ! step_matrix S ngxnb
975 0 : ALLOCATE (step_matrix(ng, nb))
976 : ! res_matrix Y ngxnb
977 0 : ALLOCATE (res_matrix(ng, nb))
978 : ! matrix A nbxnb
979 0 : ALLOCATE (a_matrix(nb, ng_global))
980 : ! PSI nb vector of norms
981 0 : ALLOCATE (norm_res(nb))
982 0 : ALLOCATE (norm_res_low(nb))
983 0 : ALLOCATE (norm_res_up(nb))
984 :
985 : ! matrix B nbxnb
986 0 : ALLOCATE (b_matrix(nb, nb))
987 : ! matrix B_inv nbxnb
988 0 : ALLOCATE (binv_matrix(nb, nb))
989 :
990 0 : ALLOCATE (gn_global(ng_global))
991 0 : 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 0 : ALLOCATE (tmp_vec(nb))
998 : END IF
999 0 : ALLOCATE (pgn(ng))
1000 0 : ALLOCATE (ugn(ng))
1001 :
1002 0 : DO ispin = 1, nspin
1003 : ! generate the global vector with the present residual
1004 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1005 0 : gn_global = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1006 0 : DO ig = 1, ng
1007 0 : ig_global = mixing_store%ig_global_index(ig)
1008 0 : gn_global(ig_global) = gn(ig)
1009 : END DO
1010 0 : 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 0 : g2 => rho_g(1)%pw_grid%gsq
1018 :
1019 0 : DO jb = 1, nb + 1
1020 0 : IF (jb < ib) THEN
1021 : kb = jb
1022 0 : ELSEIF (jb > ib) THEN
1023 0 : kb = jb - 1
1024 : ELSE
1025 : CYCLE
1026 : END IF
1027 0 : norm_res(kb) = 0.0_dp
1028 0 : norm_res_low(kb) = 0.0_dp
1029 0 : norm_res_up(kb) = 0.0_dp
1030 0 : DO ig = 1, ng
1031 :
1032 0 : prec = 1.0_dp
1033 :
1034 : step_matrix(ig, kb) = prec*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) - &
1035 0 : mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1036 : res_matrix(ig, kb) = (mixing_store%res_buffer(jb, ispin)%cc(ig) - &
1037 0 : 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 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1040 0 : 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 0 : 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 0 : AIMAG(res_matrix(ig, kb))*AIMAG(res_matrix(ig, kb))
1048 : END IF
1049 0 : 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 0 : CALL para_env%sum(norm_res)
1055 0 : CALL para_env%sum(norm_res_up)
1056 0 : CALL para_env%sum(norm_res_low)
1057 0 : 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 0 : DO jb = 1, nb
1061 0 : n_low = n_low + norm_res_low(jb)/norm_res(jb)
1062 0 : n_up = n_up + norm_res_up(jb)/norm_res(jb)
1063 : END DO
1064 0 : DO ig = 1, ng
1065 0 : IF (g2(ig) > 4.0_dp) THEN
1066 0 : step_matrix(ig, 1:nb) = step_matrix(ig, 1:nb)*SQRT(n_low/n_up)
1067 0 : res_matrix(ig, 1:nb) = res_matrix(ig, 1:nb)*SQRT(n_low/n_up)
1068 : END IF
1069 : END DO
1070 0 : DO kb = 1, nb
1071 0 : step_matrix(1:ng, kb) = step_matrix(1:ng, kb)*norm_res(kb)
1072 0 : 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 0 : DO jb = 1, nb
1078 0 : DO kb = 1, nb
1079 0 : b_matrix(kb, jb) = 0.0_dp
1080 0 : 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 0 : 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 0 : CALL para_env%sum(b_matrix)
1089 0 : DO jb = 1, nb
1090 0 : b_matrix(jb, jb) = b_matrix(jb, jb) + reg_par
1091 : END DO
1092 : ! invert B
1093 0 : CALL invert_matrix(b_matrix, binv_matrix, inv_err)
1094 :
1095 : ! A = Binv Ynorm^t
1096 0 : a_matrix = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1097 0 : DO kb = 1, nb
1098 0 : res_matrix_global = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1099 0 : DO ig = 1, ng
1100 0 : ig_global = mixing_store%ig_global_index(ig)
1101 0 : res_matrix_global(ig_global) = res_matrix(ig, kb)
1102 : END DO
1103 0 : CALL para_env%sum(res_matrix_global)
1104 :
1105 0 : DO jb = 1, nb
1106 0 : DO ig = 1, ng
1107 0 : ig_global = mixing_store%ig_global_index(ig)
1108 : a_matrix(jb, ig_global) = a_matrix(jb, ig_global) + &
1109 0 : binv_matrix(jb, kb)*res_matrix_global(ig_global)
1110 : END DO
1111 : END DO
1112 : END DO
1113 0 : CALL para_env%sum(a_matrix)
1114 :
1115 : ! compute the two components of gn that will be used to update rho
1116 0 : gn => mixing_store%res_buffer(ib, ispin)%cc
1117 0 : 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 0 : nb, gn_global(1), 1, czero, tmp_vec(1), 1)
1144 :
1145 : CALL zgemv("N", ng, nb, cmone, step_matrix(1, 1), ng, &
1146 0 : tmp_vec(1), 1, czero, pgn(1), 1)
1147 :
1148 : CALL zgemv("N", ng, nb, cmone, res_matrix(1, 1), ng, &
1149 0 : tmp_vec(1), 1, czero, ugn(1), 1)
1150 :
1151 0 : DO ig = 1, ng
1152 : pgn_norm = pgn_norm + REAL(pgn(ig), dp)*REAL(pgn(ig), dp) + &
1153 0 : AIMAG(pgn(ig))*AIMAG(pgn(ig))
1154 0 : ugn(ig) = ugn(ig) + gn(ig)
1155 : END DO
1156 0 : 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 0 : gn_norm = mixing_store%norm_res_buffer(ib, ispin)
1185 0 : gn_norm_old = mixing_store%norm_res_buffer(ib_prev, ispin)
1186 : IF (ib_prev /= 0) THEN
1187 0 : 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 0 : sigma_tilde = 0.1_dp
1192 : ! Step size for the unpredicted component
1193 0 : step_size = MIN(sigma_tilde, r_step*pgn_norm/gn_norm, sigma_max)
1194 0 : sigma_old = step_size
1195 :
1196 : ! update the density
1197 0 : DO ig = 1, ng
1198 0 : prec = mixing_store%kerker_factor(ig)
1199 : rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(ib, ispin)%cc(ig) &
1200 0 : - prec*step_size*ugn(ig) + prec*pgn(ig) ! - 0.1_dp * prec* gn(ig)
1201 0 : 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 0 : DEALLOCATE (step_matrix, res_matrix)
1208 0 : DEALLOCATE (norm_res)
1209 0 : DEALLOCATE (norm_res_up)
1210 0 : DEALLOCATE (norm_res_low)
1211 0 : DEALLOCATE (a_matrix, b_matrix, binv_matrix)
1212 0 : DEALLOCATE (ugn, pgn)
1213 : IF (use_zgemm) THEN
1214 : DEALLOCATE (sa, ya)
1215 : END IF
1216 : IF (use_zgemm_rev) THEN
1217 0 : DEALLOCATE (tmp_vec)
1218 : END IF
1219 0 : DEALLOCATE (gn_global, res_matrix_global)
1220 :
1221 0 : CALL timestop(handle)
1222 :
1223 0 : END SUBROUTINE multisecant_mixing
1224 :
1225 : END MODULE qs_gspace_mixing
|
