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