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