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