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 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 2256 : 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 2256 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_q, p_in
67 :
68 2256 : CALL timeset(routineN, handle)
69 :
70 2256 : NULLIFY (matrix_q, p_in)
71 :
72 2256 : CPASSERT(ASSOCIATED(p_out))
73 2256 : NULLIFY (matrix_q, p_in)
74 2256 : p_in => rho_ao
75 2256 : matrix_q => p_delta
76 2256 : nspins = SIZE(p_in, 1)
77 2256 : nimg = SIZE(p_in, 2)
78 :
79 : ! Compute the difference (p_out - p_in)and check convergence
80 2256 : delta = 0.0_dp
81 4698 : DO ispin = 1, nspins
82 148802 : DO ic = 1, nimg
83 144104 : 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 144104 : m3=matrix_q(ispin, ic)%matrix)
87 146546 : delta = MAX(tmp, delta)
88 : END DO
89 : END DO
90 2256 : CALL timestop(handle)
91 :
92 2256 : 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 14068 : 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 14068 : 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 14068 : POINTER :: sab_orb
129 14068 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
130 :
131 14068 : CALL timeset(routineN, handle)
132 :
133 14068 : 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 14068 : dft_control=dft_control)
138 :
139 : charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
140 14068 : .OR. dft_control%qs_control%semi_empirical
141 14068 : refmatrix => matrix_s(1, 1)%matrix
142 14068 : nimg = dft_control%nimages
143 :
144 : ! *** allocate p_mix_new ***
145 14068 : IF (PRESENT(p_mix_new)) THEN
146 14064 : IF (.NOT. ASSOCIATED(p_mix_new)) THEN
147 14028 : CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
148 29604 : DO i = 1, nspins
149 141660 : DO ic = 1, nimg
150 112056 : ALLOCATE (p_mix_new(i, ic)%matrix)
151 : CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
152 112056 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
153 112056 : CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
154 127632 : 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 14068 : IF (PRESENT(p_delta)) THEN
162 14064 : IF (mixing_method >= gspace_mixing_nr) THEN
163 270 : IF (.NOT. ASSOCIATED(p_delta)) THEN
164 266 : CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
165 564 : DO i = 1, nspins
166 14950 : DO ic = 1, nimg
167 14386 : ALLOCATE (p_delta(i, ic)%matrix)
168 : CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
169 14386 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
170 14386 : CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
171 14684 : CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
172 : END DO
173 : END DO
174 : END IF
175 270 : CPASSERT(ASSOCIATED(mixing_store))
176 : END IF
177 : END IF
178 :
179 14068 : IF (charge_mixing) THEN
180 : ! *** allocate buffer for charge mixing ***
181 8198 : 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 8198 : 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 8198 : 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 22014 : 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 510 : 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 10190 : 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 10190 : mixing_store%ubroy = 0.0_dp
229 8144 : 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 5870 : IF (mixing_method >= gspace_mixing_nr) THEN
235 220 : nbuffer = mixing_store%nbuffer
236 220 : mixing_store%ncall = 0
237 220 : IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
238 790 : ALLOCATE (mixing_store%rhoin(nspins))
239 410 : DO ispin = 1, nspins
240 410 : NULLIFY (mixing_store%rhoin(ispin)%cc)
241 : END DO
242 : END IF
243 :
244 220 : 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 5870 : IF (mixing_method >= pulay_mixing_nr) THEN
264 210 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
265 2722 : ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
266 390 : DO ispin = 1, nspins
267 2182 : DO i = 1, nbuffer
268 2002 : 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 5870 : 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 394 : ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
282 72 : DO ispin = 1, nspins
283 292 : DO i = 1, nbuffer
284 258 : 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 5870 : IF (mixing_method == broyden_mixing_nr) THEN
312 172 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
313 610 : ALLOCATE (mixing_store%rhoin_old(nspins))
314 318 : DO ispin = 1, nspins
315 318 : NULLIFY (mixing_store%rhoin_old(ispin)%cc)
316 : END DO
317 : END IF
318 172 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
319 2328 : ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
320 610 : ALLOCATE (mixing_store%last_res(nspins))
321 318 : DO ispin = 1, nspins
322 1744 : DO i = 1, nbuffer
323 1744 : NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
324 : END DO
325 318 : NULLIFY (mixing_store%last_res(ispin)%cc)
326 : END DO
327 : END IF
328 172 : 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 5870 : 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 5870 : 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 14068 : CALL timestop(handle)
382 :
383 14068 : 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 220 : 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 220 : REAL(dp), DIMENSION(:), POINTER :: g2
411 220 : REAL(dp), DIMENSION(:, :), POINTER :: g_vec
412 220 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
413 220 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
414 :
415 220 : CALL timeset(routineN, handle)
416 :
417 220 : NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
418 220 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
419 :
420 220 : nspin = SIZE(rho_g)
421 220 : ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
422 220 : nimg = SIZE(rho_ao_kp, 2)
423 220 : mixing_store%ig_max = ng
424 220 : g2 => rho_g(1)%pw_grid%gsq
425 220 : g_vec => rho_g(1)%pw_grid%g
426 :
427 220 : 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 220 : IF (.NOT. ASSOCIATED(mixing_store%kerker_factor)) THEN
437 570 : ALLOCATE (mixing_store%kerker_factor(ng))
438 : END IF
439 220 : IF (.NOT. ASSOCIATED(mixing_store%special_metric)) THEN
440 570 : ALLOCATE (mixing_store%special_metric(ng))
441 : END IF
442 220 : beta = mixing_store%beta
443 220 : kmin = 0.1_dp
444 16959746 : mixing_store%kerker_factor = 1.0_dp
445 16959746 : mixing_store%special_metric = 1.0_dp
446 220 : ig1 = 1
447 220 : IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
448 16959636 : DO ig = ig1, mixing_store%ig_max
449 16959416 : mixing_store%kerker_factor(ig) = MAX(g2(ig)/(g2(ig) + beta*beta), kmin)
450 : mixing_store%special_metric(ig) = &
451 : 1.0_dp + 50.0_dp/8.0_dp*( &
452 : 1.0_dp + COS(g_vec(1, ig)) + COS(g_vec(2, ig)) + COS(g_vec(3, ig)) + &
453 : COS(g_vec(1, ig))*COS(g_vec(2, ig)) + &
454 : COS(g_vec(2, ig))*COS(g_vec(3, ig)) + &
455 : COS(g_vec(1, ig))*COS(g_vec(3, ig)) + &
456 16959636 : COS(g_vec(1, ig))*COS(g_vec(2, ig))*COS(g_vec(3, ig)))
457 : END DO
458 :
459 220 : nbuffer = mixing_store%nbuffer
460 474 : DO ispin = 1, nspin
461 254 : IF (.NOT. ASSOCIATED(mixing_store%rhoin(ispin)%cc)) THEN
462 660 : ALLOCATE (mixing_store%rhoin(ispin)%cc(ng))
463 : END IF
464 37534282 : mixing_store%rhoin(ispin)%cc = rho_g(ispin)%array
465 :
466 254 : IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
467 42 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer(1, ispin)%cc)) THEN
468 258 : DO ib = 1, nbuffer
469 698 : ALLOCATE (mixing_store%rhoin_buffer(ib, ispin)%cc(ng))
470 : END DO
471 : END IF
472 : mixing_store%rhoin_buffer(1, ispin)%cc(1:ng) = &
473 2668074 : rho_g(ispin)%array(1:ng)
474 : END IF
475 474 : IF (ASSOCIATED(mixing_store%res_buffer)) THEN
476 244 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer(1, ispin)%cc)) THEN
477 2002 : DO ib = 1, nbuffer
478 5586 : ALLOCATE (mixing_store%res_buffer(ib, ispin)%cc(ng))
479 : END DO
480 : END IF
481 : END IF
482 : END DO
483 :
484 220 : IF (nspin == 2) THEN
485 3615010 : mixing_store%rhoin(1)%cc = rho_g(1)%array + rho_g(2)%array
486 3615010 : mixing_store%rhoin(2)%cc = rho_g(1)%array - rho_g(2)%array
487 34 : IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
488 55300 : mixing_store%rhoin_buffer(1, 1)%cc = rho_g(1)%array + rho_g(2)%array
489 55300 : mixing_store%rhoin_buffer(1, 2)%cc = rho_g(1)%array - rho_g(2)%array
490 : END IF
491 : END IF
492 :
493 220 : IF (mixing_store%gmix_p) THEN
494 20 : IF (PRESENT(rho_atom)) THEN
495 20 : natom = SIZE(rho_atom)
496 50 : DO ispin = 1, nspin
497 290 : DO iat = 1, natom
498 270 : IF (ASSOCIATED(rho_atom(iat)%cpc_s(ispin)%r_coef)) THEN
499 170 : mixing_store%paw(iat) = .TRUE.
500 170 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
501 170 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
502 170 : IF (ASSOCIATED(mixing_store%cpc_s_in)) THEN
503 170 : IF (.NOT. ASSOCIATED(mixing_store%cpc_s_in(iat, ispin)%r_coef)) THEN
504 424 : ALLOCATE (mixing_store%cpc_s_in(iat, ispin)%r_coef(n1, n2))
505 318 : ALLOCATE (mixing_store%cpc_h_in(iat, ispin)%r_coef(n1, n2))
506 : END IF
507 251330 : mixing_store%cpc_h_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_h(ispin)%r_coef
508 251330 : mixing_store%cpc_s_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_s(ispin)%r_coef
509 : END IF
510 : END IF
511 : END DO
512 : END DO
513 : END IF
514 : END IF
515 :
516 220 : IF (mixing_method == gspace_mixing_nr) THEN
517 210 : ELSEIF (mixing_method == pulay_mixing_nr) THEN
518 38 : IF (mixing_store%gmix_p .AND. PRESENT(rho_atom)) THEN
519 4 : DO ispin = 1, nspin
520 20 : DO iat = 1, natom
521 18 : IF (mixing_store%paw(iat)) THEN
522 2 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
523 2 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
524 2 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer(1, iat, ispin)%r_coef)) THEN
525 12 : DO ib = 1, nbuffer
526 40 : ALLOCATE (mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
527 30 : ALLOCATE (mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
528 30 : ALLOCATE (mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
529 32 : ALLOCATE (mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
530 : END DO
531 : END IF
532 12 : DO ib = 1, nbuffer
533 4630 : mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
534 4630 : mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
535 4630 : mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
536 4632 : mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
537 : END DO
538 : END IF
539 : END DO
540 : END DO
541 : END IF
542 172 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
543 374 : DO ispin = 1, nspin
544 202 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old(ispin)%cc)) THEN
545 516 : ALLOCATE (mixing_store%rhoin_old(ispin)%cc(ng))
546 : END IF
547 202 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer(1, ispin)%cc)) THEN
548 1744 : DO ib = 1, nbuffer
549 4888 : ALLOCATE (mixing_store%drho_buffer(ib, ispin)%cc(ng))
550 : END DO
551 516 : ALLOCATE (mixing_store%last_res(ispin)%cc(ng))
552 : END IF
553 2068 : DO ib = 1, nbuffer
554 136093372 : mixing_store%drho_buffer(ib, ispin)%cc = z_zero
555 : END DO
556 17298416 : mixing_store%last_res(ispin)%cc = z_zero
557 17298588 : mixing_store%rhoin_old(ispin)%cc = z_zero
558 : END DO
559 172 : IF (mixing_store%gmix_p) THEN
560 16 : IF (PRESENT(rho_atom)) THEN
561 42 : DO ispin = 1, nspin
562 250 : DO iat = 1, natom
563 234 : IF (mixing_store%paw(iat)) THEN
564 166 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
565 166 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
566 166 : IF (.NOT. ASSOCIATED(mixing_store%cpc_s_old(iat, ispin)%r_coef)) THEN
567 408 : ALLOCATE (mixing_store%cpc_s_old(iat, ispin)%r_coef(n1, n2))
568 306 : ALLOCATE (mixing_store%cpc_h_old(iat, ispin)%r_coef(n1, n2))
569 : END IF
570 123898 : mixing_store%cpc_h_old(iat, ispin)%r_coef = 0.0_dp
571 123898 : mixing_store%cpc_s_old(iat, ispin)%r_coef = 0.0_dp
572 166 : IF (.NOT. ASSOCIATED(mixing_store%dcpc_s_in(1, iat, ispin)%r_coef)) THEN
573 912 : DO ib = 1, nbuffer
574 3240 : ALLOCATE (mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef(n1, n2))
575 2532 : ALLOCATE (mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef(n1, n2))
576 : END DO
577 408 : ALLOCATE (mixing_store%cpc_h_lastres(iat, ispin)%r_coef(n1, n2))
578 306 : ALLOCATE (mixing_store%cpc_s_lastres(iat, ispin)%r_coef(n1, n2))
579 : END IF
580 1488 : DO ib = 1, nbuffer
581 988406 : mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef = 0.0_dp
582 988572 : mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef = 0.0_dp
583 : END DO
584 123898 : mixing_store%cpc_h_lastres(iat, ispin)%r_coef = 0.0_dp
585 123898 : mixing_store%cpc_s_lastres(iat, ispin)%r_coef = 0.0_dp
586 : END IF
587 : END DO
588 : END DO
589 : END IF
590 : END IF
591 :
592 172 : IF (.NOT. ASSOCIATED(mixing_store%p_metric)) THEN
593 438 : ALLOCATE (mixing_store%p_metric(ng))
594 146 : bconst = mixing_store%bconst
595 146 : g2min = 1.E30_dp
596 15273144 : DO ig = 1, ng
597 15273144 : IF (g2(ig) > 1.E-10_dp) g2min = MIN(g2min, g2(ig))
598 : END DO
599 146 : g2max = -1.E30_dp
600 15273144 : DO ig = 1, ng
601 15273144 : g2max = MAX(g2max, g2(ig))
602 : END DO
603 146 : CALL para_env%min(g2min)
604 146 : CALL para_env%max(g2max)
605 : ! fdamp/g2 varies between (bconst-1) and 0
606 : ! i.e. p_metric varies between bconst and 1
607 : ! fdamp = (bconst-1.0_dp)*g2min
608 146 : fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
609 15273144 : DO ig = 1, ng
610 15273144 : mixing_store%p_metric(ig) = (g2(ig) + fdamp)/MAX(g2(ig), 1.E-10_dp)
611 : END DO
612 146 : IF (rho_g(1)%pw_grid%have_g0) mixing_store%p_metric(1) = bconst
613 : END IF
614 0 : ELSEIF (mixing_method == multisecant_mixing_nr) THEN
615 0 : IF (.NOT. ASSOCIATED(mixing_store%ig_global_index)) THEN
616 0 : ALLOCATE (mixing_store%ig_global_index(ng))
617 : END IF
618 0 : mixing_store%ig_global_index = 0
619 0 : ig_count = 0
620 0 : DO iproc = 0, para_env%num_pe - 1
621 0 : IF (para_env%mepos == iproc) THEN
622 0 : DO ig = 1, ng
623 0 : ig_count = ig_count + 1
624 0 : mixing_store%ig_global_index(ig) = ig_count
625 : END DO
626 : END IF
627 0 : CALL para_env%bcast(ig_count, iproc)
628 : END DO
629 : END IF
630 :
631 220 : CALL timestop(handle)
632 :
633 220 : END SUBROUTINE mixing_init
634 :
635 : ! **************************************************************************************************
636 : !> \brief initialiation needed when charge mixing is used
637 : !> \param mixing_store ...
638 : !> \par History
639 : !> 02.2019 created [JGH]
640 : !> \author JGH
641 : ! **************************************************************************************************
642 54 : ELEMENTAL SUBROUTINE charge_mixing_init(mixing_store)
643 : TYPE(mixing_storage_type), INTENT(INOUT) :: mixing_store
644 :
645 10190 : mixing_store%acharge = 0.0_dp
646 :
647 54 : END SUBROUTINE charge_mixing_init
648 :
649 : END MODULE qs_mixing_utils
|
