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