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