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 2046 : 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 2046 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_q, p_in
66 :
67 2046 : CALL timeset(routineN, handle)
68 :
69 2046 : NULLIFY (matrix_q, p_in)
70 :
71 2046 : CPASSERT(ASSOCIATED(p_out))
72 2046 : NULLIFY (matrix_q, p_in)
73 2046 : p_in => rho_ao
74 2046 : matrix_q => p_delta
75 2046 : nspins = SIZE(p_in, 1)
76 2046 : nimg = SIZE(p_in, 2)
77 :
78 : ! Compute the difference (p_out - p_in)and check convergence
79 2046 : delta = 0.0_dp
80 4278 : DO ispin = 1, nspins
81 146476 : DO ic = 1, nimg
82 142198 : 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 142198 : m3=matrix_q(ispin, ic)%matrix)
86 144430 : delta = MAX(tmp, delta)
87 : END DO
88 : END DO
89 2046 : CALL timestop(handle)
90 :
91 2046 : 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 13966 : 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 13966 : 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 13966 : POINTER :: sab_orb
128 13966 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
129 :
130 13966 : CALL timeset(routineN, handle)
131 :
132 13966 : 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 13966 : dft_control=dft_control)
137 :
138 : charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
139 13966 : .OR. dft_control%qs_control%semi_empirical
140 13966 : refmatrix => matrix_s(1, 1)%matrix
141 13966 : nimg = dft_control%nimages
142 :
143 : ! *** allocate p_mix_new ***
144 13966 : IF (PRESENT(p_mix_new)) THEN
145 13962 : IF (.NOT. ASSOCIATED(p_mix_new)) THEN
146 13926 : CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
147 29384 : DO i = 1, nspins
148 140850 : DO ic = 1, nimg
149 111466 : ALLOCATE (p_mix_new(i, ic)%matrix)
150 : CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
151 111466 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric, nze=0)
152 111466 : CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
153 126924 : 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 13966 : IF (PRESENT(p_delta)) THEN
161 13962 : IF (mixing_method >= gspace_mixing_nr) THEN
162 246 : IF (.NOT. ASSOCIATED(p_delta)) THEN
163 242 : CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
164 516 : DO i = 1, nspins
165 15062 : DO ic = 1, nimg
166 14546 : ALLOCATE (p_delta(i, ic)%matrix)
167 : CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
168 14546 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric, nze=0)
169 14546 : CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
170 14820 : CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
171 : END DO
172 : END DO
173 : END IF
174 246 : CPASSERT(ASSOCIATED(mixing_store))
175 : END IF
176 : END IF
177 :
178 13966 : 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 5786 : IF (mixing_method >= gspace_mixing_nr) THEN
234 214 : nbuffer = mixing_store%nbuffer
235 214 : mixing_store%ncall = 0
236 214 : IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
237 766 : ALLOCATE (mixing_store%rhoin(nspins))
238 398 : DO ispin = 1, nspins
239 398 : NULLIFY (mixing_store%rhoin(ispin)%cc)
240 : END DO
241 : END IF
242 :
243 214 : 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 5786 : IF (mixing_method >= pulay_mixing_nr) THEN
263 204 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
264 2644 : ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
265 378 : DO ispin = 1, nspins
266 2122 : DO i = 1, nbuffer
267 1948 : 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 5786 : 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 5786 : IF (mixing_method == broyden_mixing_nr) THEN
311 166 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
312 586 : ALLOCATE (mixing_store%rhoin_old(nspins))
313 306 : DO ispin = 1, nspins
314 306 : NULLIFY (mixing_store%rhoin_old(ispin)%cc)
315 : END DO
316 : END IF
317 166 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
318 2250 : ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
319 586 : ALLOCATE (mixing_store%last_res(nspins))
320 306 : DO ispin = 1, nspins
321 1690 : DO i = 1, nbuffer
322 1690 : NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
323 : END DO
324 306 : NULLIFY (mixing_store%last_res(ispin)%cc)
325 : END DO
326 : END IF
327 166 : 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 5786 : 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 5786 : 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 13966 : CALL timestop(handle)
381 :
382 13966 : 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 214 : 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 214 : REAL(dp), DIMENSION(:), POINTER :: g2
410 214 : REAL(dp), DIMENSION(:, :), POINTER :: g_vec
411 214 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
412 214 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
413 :
414 214 : CALL timeset(routineN, handle)
415 :
416 214 : NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
417 214 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
418 :
419 214 : nspin = SIZE(rho_g)
420 214 : ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
421 214 : nimg = SIZE(rho_ao_kp, 2)
422 214 : mixing_store%ig_max = ng
423 214 : g2 => rho_g(1)%pw_grid%gsq
424 214 : g_vec => rho_g(1)%pw_grid%g
425 :
426 214 : 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 214 : IF (.NOT. ASSOCIATED(mixing_store%kerker_factor)) THEN
436 552 : ALLOCATE (mixing_store%kerker_factor(ng))
437 : END IF
438 214 : IF (.NOT. ASSOCIATED(mixing_store%special_metric)) THEN
439 552 : ALLOCATE (mixing_store%special_metric(ng))
440 : END IF
441 214 : beta = mixing_store%beta
442 214 : kmin = 0.1_dp
443 10465882 : mixing_store%kerker_factor = 1.0_dp
444 10465882 : mixing_store%special_metric = 1.0_dp
445 214 : ig1 = 1
446 214 : IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
447 10465775 : DO ig = ig1, mixing_store%ig_max
448 10465561 : 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 10465775 : COS(g_vec(1, ig))*COS(g_vec(2, ig))*COS(g_vec(3, ig)))
456 : END DO
457 :
458 214 : nbuffer = mixing_store%nbuffer
459 462 : DO ispin = 1, nspin
460 248 : IF (.NOT. ASSOCIATED(mixing_store%rhoin(ispin)%cc)) THEN
461 642 : ALLOCATE (mixing_store%rhoin(ispin)%cc(ng))
462 : END IF
463 24546560 : mixing_store%rhoin(ispin)%cc = rho_g(ispin)%array
464 :
465 248 : 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 462 : IF (ASSOCIATED(mixing_store%res_buffer)) THEN
475 238 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer(1, ispin)%cc)) THEN
476 1948 : DO ib = 1, nbuffer
477 5436 : ALLOCATE (mixing_store%res_buffer(ib, ispin)%cc(ng))
478 : END DO
479 : END IF
480 : END IF
481 : END DO
482 :
483 214 : 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 214 : 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 214 : IF (mixing_method == gspace_mixing_nr) THEN
516 204 : 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 166 : ELSEIF (mixing_method == broyden_mixing_nr) THEN
542 362 : DO ispin = 1, nspin
543 196 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old(ispin)%cc)) THEN
544 498 : ALLOCATE (mixing_store%rhoin_old(ispin)%cc(ng))
545 : END IF
546 196 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer(1, ispin)%cc)) THEN
547 1690 : DO ib = 1, nbuffer
548 4738 : ALLOCATE (mixing_store%drho_buffer(ib, ispin)%cc(ng))
549 : END DO
550 498 : ALLOCATE (mixing_store%last_res(ispin)%cc(ng))
551 : END IF
552 2014 : DO ib = 1, nbuffer
553 83224950 : mixing_store%drho_buffer(ib, ispin)%cc = CMPLX(0.0_dp, 0.0_dp, kind=dp)
554 : END DO
555 10804552 : mixing_store%last_res(ispin)%cc = CMPLX(0.0_dp, 0.0_dp, kind=dp)
556 10804718 : mixing_store%rhoin_old(ispin)%cc = CMPLX(0.0_dp, 0.0_dp, kind=dp)
557 : END DO
558 166 : 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 166 : IF (.NOT. ASSOCIATED(mixing_store%p_metric)) THEN
592 420 : ALLOCATE (mixing_store%p_metric(ng))
593 140 : bconst = mixing_store%bconst
594 140 : g2min = 1.E30_dp
595 8779280 : DO ig = 1, ng
596 8779280 : IF (g2(ig) > 1.E-10_dp) g2min = MIN(g2min, g2(ig))
597 : END DO
598 140 : g2max = -1.E30_dp
599 8779280 : DO ig = 1, ng
600 8779280 : g2max = MAX(g2max, g2(ig))
601 : END DO
602 140 : CALL para_env%min(g2min)
603 140 : 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 140 : fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
608 8779280 : DO ig = 1, ng
609 8779280 : mixing_store%p_metric(ig) = (g2(ig) + fdamp)/MAX(g2(ig), 1.E-10_dp)
610 : END DO
611 140 : 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 214 : CALL timestop(handle)
631 :
632 214 : 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
|
