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