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_mixing_utils
10 : USE bibliography, ONLY: Sundararaman2017,&
11 : cite_reference
12 : USE cp_control_types, ONLY: dft_control_type
13 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
14 : dbcsr_p_type,&
15 : dbcsr_set,&
16 : dbcsr_type,&
17 : dbcsr_type_symmetric
18 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
19 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
20 : USE distribution_1d_types, ONLY: distribution_1d_type
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: z_zero
23 : USE message_passing, ONLY: mp_para_env_type
24 : USE pw_types, ONLY: pw_c1d_gs_type
25 : USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
26 : gspace_mixing_nr,&
27 : mixing_storage_type,&
28 : modified_broyden_mixing_nr,&
29 : multisecant_mixing_nr,&
30 : new_pulay_mixing_nr,&
31 : pulay_mixing_nr
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type
34 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
35 : USE qs_rho_atom_types, ONLY: rho_atom_type
36 : USE qs_rho_types, ONLY: qs_rho_get,&
37 : qs_rho_type
38 : USE qs_scf_methods, ONLY: cp_sm_mix
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mixing_utils'
46 : INTEGER, PARAMETER, PRIVATE :: max_dftb_scc_vars = 1, &
47 : max_xtb_scc_vars = 28
48 :
49 : PUBLIC :: mixing_allocate, mixing_init, charge_mixing_init, &
50 : self_consistency_check
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief ...
56 : !> \param rho_ao ...
57 : !> \param p_delta ...
58 : !> \param para_env ...
59 : !> \param p_out ...
60 : !> \param delta ...
61 : ! **************************************************************************************************
62 5724 : SUBROUTINE self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
63 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao, p_delta
64 : TYPE(mp_para_env_type), POINTER :: para_env
65 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_out
66 : REAL(KIND=dp), INTENT(INOUT) :: delta
67 :
68 : CHARACTER(len=*), PARAMETER :: routineN = 'self_consistency_check'
69 :
70 : INTEGER :: handle, ic, ispin, nimg, nspins
71 : REAL(KIND=dp) :: tmp
72 5724 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_q, p_in
73 :
74 5724 : CALL timeset(routineN, handle)
75 :
76 5724 : NULLIFY (matrix_q, p_in)
77 :
78 5724 : CPASSERT(ASSOCIATED(p_out))
79 5724 : NULLIFY (matrix_q, p_in)
80 5724 : p_in => rho_ao
81 5724 : matrix_q => p_delta
82 5724 : nspins = SIZE(p_in, 1)
83 5724 : nimg = SIZE(p_in, 2)
84 :
85 : ! Compute the difference (p_out - p_in)and check convergence
86 5724 : delta = 0.0_dp
87 12298 : DO ispin = 1, nspins
88 397456 : DO ic = 1, nimg
89 385158 : CALL dbcsr_set(matrix_q(ispin, ic)%matrix, 0.0_dp)
90 : CALL cp_sm_mix(m1=p_out(ispin, ic)%matrix, m2=p_in(ispin, ic)%matrix, &
91 : p_mix=1.0_dp, delta=tmp, para_env=para_env, &
92 385158 : m3=matrix_q(ispin, ic)%matrix)
93 391732 : delta = MAX(tmp, delta)
94 : END DO
95 : END DO
96 5724 : CALL timestop(handle)
97 :
98 5724 : END SUBROUTINE self_consistency_check
99 :
100 : ! **************************************************************************************************
101 : !> \brief allocation needed when density mixing is used
102 : !> \param qs_env ...
103 : !> \param mixing_method ...
104 : !> \param p_mix_new ...
105 : !> \param p_delta ...
106 : !> \param nspins ...
107 : !> \param mixing_store ...
108 : !> \par History
109 : !> 05.2009 created [MI]
110 : !> 08.2014 kpoints [JGH]
111 : !> 02.2015 changed input to make g-space mixing available in linear scaling SCF [Patrick Seewald]
112 : !> 02.2019 charge mixing [JGH]
113 : !> \author fawzi
114 : ! **************************************************************************************************
115 19316 : SUBROUTINE mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
116 : TYPE(qs_environment_type), POINTER :: qs_env
117 : INTEGER :: mixing_method
118 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
119 : POINTER :: p_mix_new, p_delta
120 : INTEGER, INTENT(IN) :: nspins
121 : TYPE(mixing_storage_type), POINTER :: mixing_store
122 :
123 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mixing_allocate'
124 :
125 : INTEGER :: handle, i, ia, iat, ic, ikind, ispin, &
126 : max_shell, na, natom, nbuffer, nel, &
127 : nimg, nkind
128 : LOGICAL :: charge_mixing
129 19316 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
130 : TYPE(dbcsr_type), POINTER :: refmatrix
131 : TYPE(dft_control_type), POINTER :: dft_control
132 : TYPE(distribution_1d_type), POINTER :: distribution_1d
133 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
134 19316 : POINTER :: sab_orb
135 19316 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
136 :
137 19316 : CALL timeset(routineN, handle)
138 :
139 19316 : NULLIFY (matrix_s, dft_control, sab_orb, refmatrix, rho_atom)
140 : CALL get_qs_env(qs_env, &
141 : sab_orb=sab_orb, &
142 : matrix_s_kp=matrix_s, &
143 19316 : dft_control=dft_control)
144 :
145 : charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
146 19316 : .OR. dft_control%qs_control%semi_empirical
147 19316 : refmatrix => matrix_s(1, 1)%matrix
148 19316 : nimg = dft_control%nimages
149 :
150 : ! *** allocate p_mix_new ***
151 19316 : IF (PRESENT(p_mix_new)) THEN
152 19312 : IF (.NOT. ASSOCIATED(p_mix_new)) THEN
153 19276 : CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
154 40810 : DO i = 1, nspins
155 256274 : DO ic = 1, nimg
156 215464 : ALLOCATE (p_mix_new(i, ic)%matrix)
157 : CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
158 215464 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
159 215464 : CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
160 236998 : CALL dbcsr_set(p_mix_new(i, ic)%matrix, 0.0_dp)
161 : END DO
162 : END DO
163 : END IF
164 : END IF
165 :
166 : ! *** allocate p_delta ***
167 19316 : IF (PRESENT(p_delta)) THEN
168 19312 : IF (mixing_method >= gspace_mixing_nr) THEN
169 714 : IF (.NOT. ASSOCIATED(p_delta)) THEN
170 710 : CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
171 1544 : DO i = 1, nspins
172 49748 : DO ic = 1, nimg
173 48204 : ALLOCATE (p_delta(i, ic)%matrix)
174 : CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
175 48204 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
176 48204 : CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
177 49038 : CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
178 : END DO
179 : END DO
180 : END IF
181 714 : CPASSERT(ASSOCIATED(mixing_store))
182 : END IF
183 : END IF
184 :
185 19316 : IF (charge_mixing) THEN
186 : ! *** allocate buffer for TB SCC variable mixing ***
187 12404 : IF (mixing_method >= gspace_mixing_nr) THEN
188 144 : CPASSERT(.NOT. mixing_store%gmix_p)
189 144 : IF (dft_control%qs_control%dftb) THEN
190 : max_shell = max_dftb_scc_vars
191 102 : ELSEIF (dft_control%qs_control%xtb) THEN
192 : max_shell = max_xtb_scc_vars
193 : ELSE
194 0 : CPABORT('UNKNOWN METHOD')
195 : END IF
196 144 : nbuffer = mixing_store%nbuffer
197 144 : mixing_store%ncall = 0
198 144 : CALL get_qs_env(qs_env, local_particles=distribution_1d)
199 144 : nkind = SIZE(distribution_1d%n_el)
200 292 : na = SUM(distribution_1d%n_el(:))
201 144 : IF (ASSOCIATED(mixing_store%atlist)) DEALLOCATE (mixing_store%atlist)
202 432 : ALLOCATE (mixing_store%atlist(na))
203 144 : mixing_store%nat_local = na
204 144 : mixing_store%max_shell = max_shell
205 144 : ia = 0
206 292 : DO ikind = 1, nkind
207 148 : nel = distribution_1d%n_el(ikind)
208 744 : DO iat = 1, nel
209 452 : ia = ia + 1
210 600 : mixing_store%atlist(ia) = distribution_1d%list(ikind)%array(iat)
211 : END DO
212 : END DO
213 144 : IF (ASSOCIATED(mixing_store%acharge)) DEALLOCATE (mixing_store%acharge)
214 720 : ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer))
215 144 : IF (ASSOCIATED(mixing_store%dacharge)) DEALLOCATE (mixing_store%dacharge)
216 576 : ALLOCATE (mixing_store%dacharge(na, max_shell, nbuffer))
217 : END IF
218 12404 : IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
219 0 : IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
220 0 : ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
221 0 : mixing_store%pulay_matrix = 0.0_dp
222 12404 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
223 144 : IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
224 576 : ALLOCATE (mixing_store%abroy(nbuffer, nbuffer))
225 30340 : mixing_store%abroy = 0.0_dp
226 144 : IF (ASSOCIATED(mixing_store%wbroy)) DEALLOCATE (mixing_store%wbroy)
227 432 : ALLOCATE (mixing_store%wbroy(nbuffer))
228 1420 : mixing_store%wbroy = 0.0_dp
229 144 : IF (ASSOCIATED(mixing_store%dfbroy)) DEALLOCATE (mixing_store%dfbroy)
230 720 : ALLOCATE (mixing_store%dfbroy(na, max_shell, nbuffer))
231 115852 : mixing_store%dfbroy = 0.0_dp
232 144 : IF (ASSOCIATED(mixing_store%ubroy)) DEALLOCATE (mixing_store%ubroy)
233 576 : ALLOCATE (mixing_store%ubroy(na, max_shell, nbuffer))
234 115852 : mixing_store%ubroy = 0.0_dp
235 12260 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
236 0 : CPABORT("multisecant_mixing not available")
237 : END IF
238 : ELSE
239 : ! *** allocate buffer for gspace mixing ***
240 6912 : IF (mixing_method >= gspace_mixing_nr) THEN
241 574 : nbuffer = mixing_store%nbuffer
242 574 : mixing_store%ncall = 0
243 574 : IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
244 1554 : ALLOCATE (mixing_store%rhoin(nspins))
245 806 : DO ispin = 1, nspins
246 806 : NULLIFY (mixing_store%rhoin(ispin)%cc)
247 : END DO
248 : END IF
249 :
250 574 : IF (mixing_store%gmix_p .AND. dft_control%qs_control%gapw) THEN
251 20 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
252 20 : natom = SIZE(rho_atom)
253 20 : IF (.NOT. ASSOCIATED(mixing_store%paw)) THEN
254 48 : ALLOCATE (mixing_store%paw(natom))
255 144 : mixing_store%paw = .FALSE.
256 262 : ALLOCATE (mixing_store%cpc_h_in(natom, nspins))
257 246 : ALLOCATE (mixing_store%cpc_s_in(natom, nspins))
258 38 : DO ispin = 1, nspins
259 214 : DO iat = 1, natom
260 176 : NULLIFY (mixing_store%cpc_h_in(iat, ispin)%r_coef)
261 198 : NULLIFY (mixing_store%cpc_s_in(iat, ispin)%r_coef)
262 : END DO
263 : END DO
264 : END IF
265 : END IF
266 : END IF
267 :
268 : ! *** allocare res_buffer if needed
269 6912 : IF (mixing_method >= pulay_mixing_nr) THEN
270 562 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
271 4902 : ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
272 782 : DO ispin = 1, nspins
273 3816 : DO i = 1, nbuffer
274 3454 : NULLIFY (mixing_store%res_buffer(i, ispin)%cc)
275 : END DO
276 : END DO
277 : END IF
278 : END IF
279 :
280 : ! *** allocate pulay
281 6912 : IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
282 42 : IF (.NOT. ASSOCIATED(mixing_store%pulay_matrix)) THEN
283 190 : ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
284 : END IF
285 :
286 42 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
287 476 : ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
288 84 : DO ispin = 1, nspins
289 362 : DO i = 1, nbuffer
290 324 : NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
291 : END DO
292 : END DO
293 : END IF
294 42 : IF (mixing_store%gmix_p) THEN
295 2 : IF (dft_control%qs_control%gapw) THEN
296 2 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer)) THEN
297 108 : ALLOCATE (mixing_store%cpc_h_in_buffer(nbuffer, natom, nspins))
298 106 : ALLOCATE (mixing_store%cpc_s_in_buffer(nbuffer, natom, nspins))
299 106 : ALLOCATE (mixing_store%cpc_h_res_buffer(nbuffer, natom, nspins))
300 106 : ALLOCATE (mixing_store%cpc_s_res_buffer(nbuffer, natom, nspins))
301 4 : DO ispin = 1, nspins
302 20 : DO iat = 1, natom
303 98 : DO i = 1, nbuffer
304 80 : NULLIFY (mixing_store%cpc_h_in_buffer(i, iat, ispin)%r_coef)
305 80 : NULLIFY (mixing_store%cpc_s_in_buffer(i, iat, ispin)%r_coef)
306 80 : NULLIFY (mixing_store%cpc_h_res_buffer(i, iat, ispin)%r_coef)
307 96 : NULLIFY (mixing_store%cpc_s_res_buffer(i, iat, ispin)%r_coef)
308 : END DO
309 : END DO
310 : END DO
311 : END IF
312 : END IF
313 : END IF
314 :
315 : END IF
316 : ! *** allocate broyden buffer ***
317 6912 : IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
318 520 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
319 1346 : ALLOCATE (mixing_store%rhoin_old(nspins))
320 698 : DO ispin = 1, nspins
321 698 : NULLIFY (mixing_store%rhoin_old(ispin)%cc)
322 : END DO
323 : END IF
324 520 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
325 4426 : ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
326 1346 : ALLOCATE (mixing_store%last_res(nspins))
327 698 : DO ispin = 1, nspins
328 3130 : DO i = 1, nbuffer
329 3130 : NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
330 : END DO
331 698 : NULLIFY (mixing_store%last_res(ispin)%cc)
332 : END DO
333 : END IF
334 520 : IF (mixing_method == modified_broyden_mixing_nr) THEN
335 2 : IF (.NOT. ASSOCIATED(mixing_store%weight)) THEN
336 8 : ALLOCATE (mixing_store%weight(nbuffer, nspins))
337 : END IF
338 34 : mixing_store%weight = 0.0_dp
339 : END IF
340 520 : IF (mixing_store%gmix_p) THEN
341 :
342 16 : IF (dft_control%qs_control%gapw) THEN
343 16 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_old)) THEN
344 210 : ALLOCATE (mixing_store%cpc_h_old(natom, nspins))
345 198 : ALLOCATE (mixing_store%cpc_s_old(natom, nspins))
346 30 : DO ispin = 1, nspins
347 174 : DO iat = 1, natom
348 144 : NULLIFY (mixing_store%cpc_h_old(iat, ispin)%r_coef)
349 162 : NULLIFY (mixing_store%cpc_s_old(iat, ispin)%r_coef)
350 : END DO
351 : END DO
352 : END IF
353 16 : IF (.NOT. ASSOCIATED(mixing_store%dcpc_h_in)) THEN
354 1326 : ALLOCATE (mixing_store%dcpc_h_in(nbuffer, natom, nspins))
355 1314 : ALLOCATE (mixing_store%dcpc_s_in(nbuffer, natom, nspins))
356 210 : ALLOCATE (mixing_store%cpc_h_lastres(natom, nspins))
357 198 : ALLOCATE (mixing_store%cpc_s_lastres(natom, nspins))
358 30 : DO ispin = 1, nspins
359 174 : DO iat = 1, natom
360 1248 : DO i = 1, nbuffer
361 1104 : NULLIFY (mixing_store%dcpc_h_in(i, iat, ispin)%r_coef)
362 1248 : NULLIFY (mixing_store%dcpc_s_in(i, iat, ispin)%r_coef)
363 : END DO
364 144 : NULLIFY (mixing_store%cpc_h_lastres(iat, ispin)%r_coef)
365 162 : NULLIFY (mixing_store%cpc_s_lastres(iat, ispin)%r_coef)
366 : END DO
367 : END DO
368 : END IF
369 : END IF
370 : END IF
371 : END IF
372 :
373 : ! *** allocate multisecant buffer ***
374 6912 : IF (mixing_method == multisecant_mixing_nr) THEN
375 0 : IF (.NOT. ASSOCIATED(mixing_store%norm_res_buffer)) THEN
376 0 : ALLOCATE (mixing_store%norm_res_buffer(nbuffer, nspins))
377 : END IF
378 : END IF
379 :
380 6912 : IF (mixing_method == multisecant_mixing_nr) THEN
381 0 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
382 0 : ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
383 0 : DO ispin = 1, nspins
384 0 : DO i = 1, nbuffer
385 0 : NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
386 : END DO
387 : END DO
388 : END IF
389 : END IF
390 :
391 : END IF
392 :
393 19316 : CALL timestop(handle)
394 :
395 19316 : END SUBROUTINE mixing_allocate
396 :
397 : ! **************************************************************************************************
398 : !> \brief initialiation needed when gspace mixing is used
399 : !> \param mixing_method ...
400 : !> \param rho ...
401 : !> \param mixing_store ...
402 : !> \param para_env ...
403 : !> \param rho_atom ...
404 : !> \par History
405 : !> 05.2009 created [MI]
406 : !> \author MI
407 : ! **************************************************************************************************
408 574 : SUBROUTINE mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
409 : INTEGER, INTENT(IN) :: mixing_method
410 : TYPE(qs_rho_type), POINTER :: rho
411 : TYPE(mixing_storage_type), POINTER :: mixing_store
412 : TYPE(mp_para_env_type), POINTER :: para_env
413 : TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
414 : POINTER :: rho_atom
415 :
416 : CHARACTER(len=*), PARAMETER :: routineN = 'mixing_init'
417 :
418 : INTEGER :: handle, iat, ib, icoef1, icoef2, ig, ig1, ig_count, iproc, ispin, n1, n2, natom, &
419 : nbuffer, ncoef_h1, ncoef_h2, ncoef_s1, ncoef_s2, ng, nimg, nspin
420 : REAL(dp) :: bconst, beta, cpc_h_tmp, cpc_s_tmp, &
421 : fdamp, g2max, g2min, kmin, qk, qkappa, &
422 : qm
423 574 : REAL(dp), DIMENSION(:), POINTER :: g2
424 574 : REAL(dp), DIMENSION(:, :), POINTER :: g_vec
425 574 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
426 574 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
427 :
428 574 : CALL timeset(routineN, handle)
429 :
430 574 : NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
431 574 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
432 :
433 574 : nspin = SIZE(rho_g)
434 574 : ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
435 574 : nimg = SIZE(rho_ao_kp, 2)
436 574 : mixing_store%ig_max = ng
437 574 : g2 => rho_g(1)%pw_grid%gsq
438 574 : g_vec => rho_g(1)%pw_grid%g
439 :
440 574 : IF (mixing_store%max_gvec_exp > 0._dp) THEN
441 0 : DO ig = 1, ng
442 0 : IF (g2(ig) > mixing_store%max_g2) THEN
443 0 : mixing_store%ig_max = ig
444 0 : EXIT
445 : END IF
446 : END DO
447 : END IF
448 :
449 574 : IF (.NOT. ASSOCIATED(mixing_store%kerker_factor)) THEN
450 1122 : ALLOCATE (mixing_store%kerker_factor(ng))
451 : END IF
452 574 : IF (.NOT. ASSOCIATED(mixing_store%kerker_factor_mag)) THEN
453 1122 : ALLOCATE (mixing_store%kerker_factor_mag(ng))
454 : END IF
455 574 : IF (.NOT. ASSOCIATED(mixing_store%special_metric)) THEN
456 1122 : ALLOCATE (mixing_store%special_metric(ng))
457 : END IF
458 574 : beta = mixing_store%beta
459 574 : kmin = 0.1_dp
460 20157372 : mixing_store%kerker_factor = 1.0_dp
461 20157372 : mixing_store%kerker_factor_mag = 1.0_dp
462 20157372 : mixing_store%special_metric = 1.0_dp
463 574 : ig1 = 1
464 :
465 574 : IF (.NOT. (mixing_method == new_pulay_mixing_nr)) THEN
466 570 : IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
467 19796123 : DO ig = ig1, mixing_store%ig_max
468 19795553 : mixing_store%kerker_factor(ig) = MAX(g2(ig)/(g2(ig) + beta*beta), kmin)
469 19795553 : IF (mixing_store%beta_mag > 1.0E-10_dp) THEN
470 19543322 : mixing_store%kerker_factor_mag(ig) = MAX(g2(ig)/(g2(ig) + mixing_store%beta_mag*mixing_store%beta_mag), kmin)
471 : ELSE
472 : ! beta_mag ~ 0 means no Kerker screening on the magnetization channel
473 252231 : mixing_store%kerker_factor_mag(ig) = 1.0_dp
474 : END IF
475 : mixing_store%special_metric(ig) = &
476 : 1.0_dp + 50.0_dp/8.0_dp*( &
477 : 1.0_dp + COS(g_vec(1, ig)) + COS(g_vec(2, ig)) + COS(g_vec(3, ig)) + &
478 : COS(g_vec(1, ig))*COS(g_vec(2, ig)) + &
479 : COS(g_vec(2, ig))*COS(g_vec(3, ig)) + &
480 : COS(g_vec(1, ig))*COS(g_vec(3, ig)) + &
481 19796123 : COS(g_vec(1, ig))*COS(g_vec(2, ig))*COS(g_vec(3, ig)))
482 : END DO
483 : ELSE
484 4 : CALL cite_reference(Sundararaman2017)
485 4 : qkappa = mixing_store%qkappa
486 4 : qk = mixing_store%qk
487 4 : qm = mixing_store%qm
488 : IF (rho_g(1)%pw_grid%have_g0) ig1 = 1
489 360964 : DO ig = ig1, mixing_store%ig_max
490 360960 : mixing_store%kerker_factor(ig) = (g2(ig) + qkappa*qkappa)/(g2(ig) + qkappa*qkappa + qk*qk)
491 360964 : mixing_store%special_metric(ig) = (g2(ig) + qkappa*qkappa + qm*qm)/(g2(ig) + qkappa*qkappa)
492 : END DO
493 : END IF
494 574 : nbuffer = mixing_store%nbuffer
495 1272 : DO ispin = 1, nspin
496 698 : IF (.NOT. ASSOCIATED(mixing_store%rhoin(ispin)%cc)) THEN
497 1296 : ALLOCATE (mixing_store%rhoin(ispin)%cc(ng))
498 : END IF
499 49274570 : mixing_store%rhoin(ispin)%cc = rho_g(ispin)%array
500 :
501 698 : IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
502 50 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer(1, ispin)%cc)) THEN
503 324 : DO ib = 1, nbuffer
504 880 : ALLOCATE (mixing_store%rhoin_buffer(ib, ispin)%cc(ng))
505 : END DO
506 : END IF
507 : mixing_store%rhoin_buffer(1, ispin)%cc(1:ng) = &
508 4111922 : rho_g(ispin)%array(1:ng)
509 : END IF
510 1272 : IF (ASSOCIATED(mixing_store%res_buffer)) THEN
511 686 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer(1, ispin)%cc)) THEN
512 3454 : DO ib = 1, nbuffer
513 9522 : ALLOCATE (mixing_store%res_buffer(ib, ispin)%cc(ng))
514 : END DO
515 : END IF
516 : END IF
517 : END DO
518 :
519 574 : IF (nspin == 2) THEN
520 8960400 : mixing_store%rhoin(1)%cc = rho_g(1)%array + rho_g(2)%array
521 8960400 : mixing_store%rhoin(2)%cc = rho_g(1)%array - rho_g(2)%array
522 124 : IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
523 777224 : mixing_store%rhoin_buffer(1, 1)%cc = rho_g(1)%array + rho_g(2)%array
524 777224 : mixing_store%rhoin_buffer(1, 2)%cc = rho_g(1)%array - rho_g(2)%array
525 : END IF
526 : END IF
527 :
528 574 : IF (mixing_store%gmix_p) THEN
529 20 : IF (PRESENT(rho_atom)) THEN
530 20 : natom = SIZE(rho_atom)
531 50 : DO ispin = 1, nspin
532 290 : DO iat = 1, natom
533 270 : IF (ASSOCIATED(rho_atom(iat)%cpc_s(ispin)%r_coef)) THEN
534 170 : mixing_store%paw(iat) = .TRUE.
535 170 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
536 170 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
537 170 : IF (ASSOCIATED(mixing_store%cpc_s_in)) THEN
538 170 : IF (.NOT. ASSOCIATED(mixing_store%cpc_s_in(iat, ispin)%r_coef)) THEN
539 424 : ALLOCATE (mixing_store%cpc_s_in(iat, ispin)%r_coef(n1, n2))
540 318 : ALLOCATE (mixing_store%cpc_h_in(iat, ispin)%r_coef(n1, n2))
541 : END IF
542 251330 : mixing_store%cpc_h_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_h(ispin)%r_coef
543 251330 : mixing_store%cpc_s_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_s(ispin)%r_coef
544 : END IF
545 : END IF
546 : END DO
547 : END DO
548 20 : IF (nspin == 2 .AND. ASSOCIATED(mixing_store%cpc_s_in)) THEN
549 90 : DO iat = 1, natom
550 90 : IF (mixing_store%paw(iat)) THEN
551 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_h_in(iat, 1)%r_coef))
552 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_h_in(iat, 2)%r_coef))
553 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_s_in(iat, 1)%r_coef))
554 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_s_in(iat, 2)%r_coef))
555 :
556 80 : ncoef_h1 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 1)
557 80 : ncoef_h2 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 2)
558 80 : ncoef_s1 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 1)
559 80 : ncoef_s2 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 2)
560 80 : CPASSERT(ncoef_h1 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 1))
561 80 : CPASSERT(ncoef_h2 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 2))
562 80 : CPASSERT(ncoef_s1 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 1))
563 80 : CPASSERT(ncoef_s2 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 2))
564 :
565 2240 : DO icoef2 = 1, ncoef_h2
566 60560 : DO icoef1 = 1, ncoef_h1
567 58320 : cpc_h_tmp = mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2)
568 : mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2) = &
569 58320 : cpc_h_tmp + mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
570 : mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2) = &
571 60480 : cpc_h_tmp - mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
572 : END DO
573 : END DO
574 :
575 2240 : DO icoef2 = 1, ncoef_s2
576 60560 : DO icoef1 = 1, ncoef_s1
577 58320 : cpc_s_tmp = mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2)
578 : mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2) = &
579 58320 : cpc_s_tmp + mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
580 : mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2) = &
581 60480 : cpc_s_tmp - mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
582 : END DO
583 : END DO
584 : END IF
585 : END DO
586 : END IF
587 : END IF
588 : END IF
589 :
590 : IF (mixing_method == gspace_mixing_nr) THEN
591 : ELSEIF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
592 42 : IF (mixing_store%gmix_p .AND. PRESENT(rho_atom)) THEN
593 4 : DO ispin = 1, nspin
594 20 : DO iat = 1, natom
595 18 : IF (mixing_store%paw(iat)) THEN
596 2 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
597 2 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
598 2 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer(1, iat, ispin)%r_coef)) THEN
599 12 : DO ib = 1, nbuffer
600 40 : ALLOCATE (mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
601 30 : ALLOCATE (mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
602 30 : ALLOCATE (mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
603 32 : ALLOCATE (mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
604 : END DO
605 : END IF
606 12 : DO ib = 1, nbuffer
607 4630 : mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
608 4630 : mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
609 4630 : mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
610 4632 : mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
611 : END DO
612 : END IF
613 : END DO
614 : END DO
615 : END IF
616 : ELSEIF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
617 1156 : DO ispin = 1, nspin
618 636 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old(ispin)%cc)) THEN
619 1122 : ALLOCATE (mixing_store%rhoin_old(ispin)%cc(ng))
620 : END IF
621 636 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer(1, ispin)%cc)) THEN
622 3130 : DO ib = 1, nbuffer
623 8642 : ALLOCATE (mixing_store%drho_buffer(ib, ispin)%cc(ng))
624 : END DO
625 1122 : ALLOCATE (mixing_store%last_res(ispin)%cc(ng))
626 : END IF
627 5430 : DO ib = 1, nbuffer
628 167369702 : mixing_store%drho_buffer(ib, ispin)%cc = z_zero
629 : END DO
630 22382852 : mixing_store%last_res(ispin)%cc = z_zero
631 22383372 : mixing_store%rhoin_old(ispin)%cc = z_zero
632 : END DO
633 520 : IF (mixing_store%gmix_p) THEN
634 16 : IF (PRESENT(rho_atom)) THEN
635 42 : DO ispin = 1, nspin
636 250 : DO iat = 1, natom
637 234 : IF (mixing_store%paw(iat)) THEN
638 166 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
639 166 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
640 166 : IF (.NOT. ASSOCIATED(mixing_store%cpc_s_old(iat, ispin)%r_coef)) THEN
641 408 : ALLOCATE (mixing_store%cpc_s_old(iat, ispin)%r_coef(n1, n2))
642 306 : ALLOCATE (mixing_store%cpc_h_old(iat, ispin)%r_coef(n1, n2))
643 : END IF
644 123898 : mixing_store%cpc_h_old(iat, ispin)%r_coef = 0.0_dp
645 123898 : mixing_store%cpc_s_old(iat, ispin)%r_coef = 0.0_dp
646 166 : IF (.NOT. ASSOCIATED(mixing_store%dcpc_s_in(1, iat, ispin)%r_coef)) THEN
647 912 : DO ib = 1, nbuffer
648 3240 : ALLOCATE (mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef(n1, n2))
649 2532 : ALLOCATE (mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef(n1, n2))
650 : END DO
651 408 : ALLOCATE (mixing_store%cpc_h_lastres(iat, ispin)%r_coef(n1, n2))
652 306 : ALLOCATE (mixing_store%cpc_s_lastres(iat, ispin)%r_coef(n1, n2))
653 : END IF
654 1488 : DO ib = 1, nbuffer
655 988406 : mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef = 0.0_dp
656 988572 : mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef = 0.0_dp
657 : END DO
658 123898 : mixing_store%cpc_h_lastres(iat, ispin)%r_coef = 0.0_dp
659 123898 : mixing_store%cpc_s_lastres(iat, ispin)%r_coef = 0.0_dp
660 : END IF
661 : END DO
662 : END DO
663 : END IF
664 : END IF
665 :
666 520 : IF (.NOT. ASSOCIATED(mixing_store%p_metric)) THEN
667 972 : ALLOCATE (mixing_store%p_metric(ng))
668 324 : bconst = mixing_store%bconst
669 324 : g2min = 1.E30_dp
670 15656642 : DO ig = 1, ng
671 15656642 : IF (g2(ig) > 1.E-10_dp) g2min = MIN(g2min, g2(ig))
672 : END DO
673 324 : g2max = -1.E30_dp
674 15656642 : DO ig = 1, ng
675 15656642 : g2max = MAX(g2max, g2(ig))
676 : END DO
677 324 : CALL para_env%min(g2min)
678 324 : CALL para_env%max(g2max)
679 : ! fdamp/g2 varies between (bconst-1) and 0
680 : ! i.e. p_metric varies between bconst and 1
681 : ! fdamp = (bconst-1.0_dp)*g2min
682 324 : fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
683 15656642 : DO ig = 1, ng
684 15656642 : mixing_store%p_metric(ig) = (g2(ig) + fdamp)/MAX(g2(ig), 1.E-10_dp)
685 : END DO
686 324 : IF (rho_g(1)%pw_grid%have_g0) mixing_store%p_metric(1) = bconst
687 : END IF
688 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
689 0 : IF (.NOT. ASSOCIATED(mixing_store%ig_global_index)) THEN
690 0 : ALLOCATE (mixing_store%ig_global_index(ng))
691 : END IF
692 0 : mixing_store%ig_global_index = 0
693 0 : ig_count = 0
694 0 : DO iproc = 0, para_env%num_pe - 1
695 0 : IF (para_env%mepos == iproc) THEN
696 0 : DO ig = 1, ng
697 0 : ig_count = ig_count + 1
698 0 : mixing_store%ig_global_index(ig) = ig_count
699 : END DO
700 : END IF
701 0 : CALL para_env%bcast(ig_count, iproc)
702 : END DO
703 : END IF
704 :
705 574 : CALL timestop(handle)
706 :
707 574 : END SUBROUTINE mixing_init
708 :
709 : ! **************************************************************************************************
710 : !> \brief initialiation needed when charge mixing is used
711 : !> \param mixing_store ...
712 : !> \par History
713 : !> 02.2019 created [JGH]
714 : !> \author JGH
715 : ! **************************************************************************************************
716 144 : ELEMENTAL SUBROUTINE charge_mixing_init(mixing_store)
717 : TYPE(mixing_storage_type), INTENT(INOUT) :: mixing_store
718 :
719 144 : mixing_store%ncall = 0
720 144 : mixing_store%tb_scc_mixer_error = 0.0_dp
721 144 : mixing_store%tb_scc_mixer_step = 0
722 115852 : IF (ASSOCIATED(mixing_store%acharge)) mixing_store%acharge = 0.0_dp
723 115852 : IF (ASSOCIATED(mixing_store%dacharge)) mixing_store%dacharge = 0.0_dp
724 :
725 144 : END SUBROUTINE charge_mixing_init
726 :
727 : END MODULE qs_mixing_utils
|
