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