(git:618e178)
Loading...
Searching...
No Matches
qs_mixing_utils.F
Go to the documentation of this file.
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! **************************************************************************************************
9MODULE 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 modified_broyden_mixing_nr,&
27 multisecant_mixing_nr,&
28 pulay_mixing_nr
29 USE qs_environment_types, ONLY: get_qs_env,&
30 qs_environment_type
31 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
32 USE qs_rho_atom_types, ONLY: rho_atom_type
33 USE qs_rho_types, ONLY: qs_rho_get,&
34 qs_rho_type
35 USE qs_scf_methods, ONLY: cp_sm_mix
36#include "./base/base_uses.f90"
37
38 IMPLICIT NONE
39
40 PRIVATE
41
42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mixing_utils'
43 INTEGER, PARAMETER, PRIVATE :: max_dftb_scc_vars = 1, &
44 max_xtb_scc_vars = 28
45
46 PUBLIC :: mixing_allocate, mixing_init, charge_mixing_init, &
47 self_consistency_check
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief ...
53!> \param rho_ao ...
54!> \param p_delta ...
55!> \param para_env ...
56!> \param p_out ...
57!> \param delta ...
58! **************************************************************************************************
59 SUBROUTINE self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
60 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao, p_delta
61 TYPE(mp_para_env_type), POINTER :: para_env
62 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_out
63 REAL(kind=dp), INTENT(INOUT) :: delta
64
65 CHARACTER(len=*), PARAMETER :: routinen = 'self_consistency_check'
66
67 INTEGER :: handle, ic, ispin, nimg, nspins
68 REAL(kind=dp) :: tmp
69 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_q, p_in
70
71 CALL timeset(routinen, handle)
72
73 NULLIFY (matrix_q, p_in)
74
75 cpassert(ASSOCIATED(p_out))
76 NULLIFY (matrix_q, p_in)
77 p_in => rho_ao
78 matrix_q => p_delta
79 nspins = SIZE(p_in, 1)
80 nimg = SIZE(p_in, 2)
81
82 ! Compute the difference (p_out - p_in)and check convergence
83 delta = 0.0_dp
84 DO ispin = 1, nspins
85 DO ic = 1, nimg
86 CALL dbcsr_set(matrix_q(ispin, ic)%matrix, 0.0_dp)
87 CALL cp_sm_mix(m1=p_out(ispin, ic)%matrix, m2=p_in(ispin, ic)%matrix, &
88 p_mix=1.0_dp, delta=tmp, para_env=para_env, &
89 m3=matrix_q(ispin, ic)%matrix)
90 delta = max(tmp, delta)
91 END DO
92 END DO
93 CALL timestop(handle)
94
95 END SUBROUTINE self_consistency_check
96
97! **************************************************************************************************
98!> \brief allocation needed when density mixing is used
99!> \param qs_env ...
100!> \param mixing_method ...
101!> \param p_mix_new ...
102!> \param p_delta ...
103!> \param nspins ...
104!> \param mixing_store ...
105!> \par History
106!> 05.2009 created [MI]
107!> 08.2014 kpoints [JGH]
108!> 02.2015 changed input to make g-space mixing available in linear scaling SCF [Patrick Seewald]
109!> 02.2019 charge mixing [JGH]
110!> \author fawzi
111! **************************************************************************************************
112 SUBROUTINE mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
113 TYPE(qs_environment_type), POINTER :: qs_env
114 INTEGER :: mixing_method
115 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
116 POINTER :: p_mix_new, p_delta
117 INTEGER, INTENT(IN) :: nspins
118 TYPE(mixing_storage_type), POINTER :: mixing_store
119
120 CHARACTER(LEN=*), PARAMETER :: routinen = 'mixing_allocate'
121
122 INTEGER :: handle, i, ia, iat, ic, ikind, ispin, &
123 max_shell, na, natom, nbuffer, nel, &
124 nimg, nkind
125 LOGICAL :: charge_mixing
126 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
127 TYPE(dbcsr_type), POINTER :: refmatrix
128 TYPE(dft_control_type), POINTER :: dft_control
129 TYPE(distribution_1d_type), POINTER :: distribution_1d
130 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
131 POINTER :: sab_orb
132 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
133
134 CALL timeset(routinen, handle)
135
136 NULLIFY (matrix_s, dft_control, sab_orb, refmatrix, rho_atom)
137 CALL get_qs_env(qs_env, &
138 sab_orb=sab_orb, &
139 matrix_s_kp=matrix_s, &
140 dft_control=dft_control)
141
142 charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
143 .OR. dft_control%qs_control%semi_empirical
144 refmatrix => matrix_s(1, 1)%matrix
145 nimg = dft_control%nimages
146
147 ! *** allocate p_mix_new ***
148 IF (PRESENT(p_mix_new)) THEN
149 IF (.NOT. ASSOCIATED(p_mix_new)) THEN
150 CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
151 DO i = 1, nspins
152 DO ic = 1, nimg
153 ALLOCATE (p_mix_new(i, ic)%matrix)
154 CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
155 name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
156 CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
157 CALL dbcsr_set(p_mix_new(i, ic)%matrix, 0.0_dp)
158 END DO
159 END DO
160 END IF
161 END IF
162
163 ! *** allocate p_delta ***
164 IF (PRESENT(p_delta)) THEN
165 IF (mixing_method >= gspace_mixing_nr) THEN
166 IF (.NOT. ASSOCIATED(p_delta)) THEN
167 CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
168 DO i = 1, nspins
169 DO ic = 1, nimg
170 ALLOCATE (p_delta(i, ic)%matrix)
171 CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
172 name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
173 CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
174 CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
175 END DO
176 END DO
177 END IF
178 cpassert(ASSOCIATED(mixing_store))
179 END IF
180 END IF
181
182 IF (charge_mixing) THEN
183 ! *** allocate buffer for TB SCC variable mixing ***
184 IF (mixing_method >= gspace_mixing_nr) THEN
185 cpassert(.NOT. mixing_store%gmix_p)
186 IF (dft_control%qs_control%dftb) THEN
187 max_shell = max_dftb_scc_vars
188 ELSEIF (dft_control%qs_control%xtb) THEN
189 max_shell = max_xtb_scc_vars
190 ELSE
191 cpabort('UNKNOWN METHOD')
192 END IF
193 nbuffer = mixing_store%nbuffer
194 mixing_store%ncall = 0
195 CALL get_qs_env(qs_env, local_particles=distribution_1d)
196 nkind = SIZE(distribution_1d%n_el)
197 na = sum(distribution_1d%n_el(:))
198 IF (ASSOCIATED(mixing_store%atlist)) DEALLOCATE (mixing_store%atlist)
199 ALLOCATE (mixing_store%atlist(na))
200 mixing_store%nat_local = na
201 mixing_store%max_shell = max_shell
202 ia = 0
203 DO ikind = 1, nkind
204 nel = distribution_1d%n_el(ikind)
205 DO iat = 1, nel
206 ia = ia + 1
207 mixing_store%atlist(ia) = distribution_1d%list(ikind)%array(iat)
208 END DO
209 END DO
210 IF (ASSOCIATED(mixing_store%acharge)) DEALLOCATE (mixing_store%acharge)
211 ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer))
212 IF (ASSOCIATED(mixing_store%dacharge)) DEALLOCATE (mixing_store%dacharge)
213 ALLOCATE (mixing_store%dacharge(na, max_shell, nbuffer))
214 END IF
215 IF (mixing_method == pulay_mixing_nr) THEN
216 IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
217 ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
218 mixing_store%pulay_matrix = 0.0_dp
219 ELSEIF (mixing_method == broyden_mixing_nr) THEN
220 IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
221 ALLOCATE (mixing_store%abroy(nbuffer, nbuffer))
222 mixing_store%abroy = 0.0_dp
223 IF (ASSOCIATED(mixing_store%wbroy)) DEALLOCATE (mixing_store%wbroy)
224 ALLOCATE (mixing_store%wbroy(nbuffer))
225 mixing_store%wbroy = 0.0_dp
226 IF (ASSOCIATED(mixing_store%dfbroy)) DEALLOCATE (mixing_store%dfbroy)
227 ALLOCATE (mixing_store%dfbroy(na, max_shell, nbuffer))
228 mixing_store%dfbroy = 0.0_dp
229 IF (ASSOCIATED(mixing_store%ubroy)) DEALLOCATE (mixing_store%ubroy)
230 ALLOCATE (mixing_store%ubroy(na, max_shell, nbuffer))
231 mixing_store%ubroy = 0.0_dp
232 ELSEIF (mixing_method == multisecant_mixing_nr) THEN
233 cpabort("multisecant_mixing not available")
234 END IF
235 ELSE
236 ! *** allocate buffer for gspace mixing ***
237 IF (mixing_method >= gspace_mixing_nr) THEN
238 nbuffer = mixing_store%nbuffer
239 mixing_store%ncall = 0
240 IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
241 ALLOCATE (mixing_store%rhoin(nspins))
242 DO ispin = 1, nspins
243 NULLIFY (mixing_store%rhoin(ispin)%cc)
244 END DO
245 END IF
246
247 IF (mixing_store%gmix_p .AND. dft_control%qs_control%gapw) THEN
248 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
249 natom = SIZE(rho_atom)
250 IF (.NOT. ASSOCIATED(mixing_store%paw)) THEN
251 ALLOCATE (mixing_store%paw(natom))
252 mixing_store%paw = .false.
253 ALLOCATE (mixing_store%cpc_h_in(natom, nspins))
254 ALLOCATE (mixing_store%cpc_s_in(natom, nspins))
255 DO ispin = 1, nspins
256 DO iat = 1, natom
257 NULLIFY (mixing_store%cpc_h_in(iat, ispin)%r_coef)
258 NULLIFY (mixing_store%cpc_s_in(iat, ispin)%r_coef)
259 END DO
260 END DO
261 END IF
262 END IF
263 END IF
264
265 ! *** allocare res_buffer if needed
266 IF (mixing_method >= pulay_mixing_nr) THEN
267 IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
268 ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
269 DO ispin = 1, nspins
270 DO i = 1, nbuffer
271 NULLIFY (mixing_store%res_buffer(i, ispin)%cc)
272 END DO
273 END DO
274 END IF
275 END IF
276
277 ! *** allocate pulay
278 IF (mixing_method == pulay_mixing_nr) THEN
279 IF (.NOT. ASSOCIATED(mixing_store%pulay_matrix)) THEN
280 ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
281 END IF
282
283 IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
284 ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
285 DO ispin = 1, nspins
286 DO i = 1, nbuffer
287 NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
288 END DO
289 END DO
290 END IF
291 IF (mixing_store%gmix_p) THEN
292 IF (dft_control%qs_control%gapw) THEN
293 IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer)) THEN
294 ALLOCATE (mixing_store%cpc_h_in_buffer(nbuffer, natom, nspins))
295 ALLOCATE (mixing_store%cpc_s_in_buffer(nbuffer, natom, nspins))
296 ALLOCATE (mixing_store%cpc_h_res_buffer(nbuffer, natom, nspins))
297 ALLOCATE (mixing_store%cpc_s_res_buffer(nbuffer, natom, nspins))
298 DO ispin = 1, nspins
299 DO iat = 1, natom
300 DO i = 1, nbuffer
301 NULLIFY (mixing_store%cpc_h_in_buffer(i, iat, ispin)%r_coef)
302 NULLIFY (mixing_store%cpc_s_in_buffer(i, iat, ispin)%r_coef)
303 NULLIFY (mixing_store%cpc_h_res_buffer(i, iat, ispin)%r_coef)
304 NULLIFY (mixing_store%cpc_s_res_buffer(i, iat, ispin)%r_coef)
305 END DO
306 END DO
307 END DO
308 END IF
309 END IF
310 END IF
311
312 END IF
313 ! *** allocate broyden buffer ***
314 IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
315 IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
316 ALLOCATE (mixing_store%rhoin_old(nspins))
317 DO ispin = 1, nspins
318 NULLIFY (mixing_store%rhoin_old(ispin)%cc)
319 END DO
320 END IF
321 IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
322 ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
323 ALLOCATE (mixing_store%last_res(nspins))
324 DO ispin = 1, nspins
325 DO i = 1, nbuffer
326 NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
327 END DO
328 NULLIFY (mixing_store%last_res(ispin)%cc)
329 END DO
330 END IF
331 IF (mixing_method == modified_broyden_mixing_nr) THEN
332 IF (.NOT. ASSOCIATED(mixing_store%weight)) THEN
333 ALLOCATE (mixing_store%weight(nbuffer, nspins))
334 END IF
335 mixing_store%weight = 0.0_dp
336 END IF
337 IF (mixing_store%gmix_p) THEN
338
339 IF (dft_control%qs_control%gapw) THEN
340 IF (.NOT. ASSOCIATED(mixing_store%cpc_h_old)) THEN
341 ALLOCATE (mixing_store%cpc_h_old(natom, nspins))
342 ALLOCATE (mixing_store%cpc_s_old(natom, nspins))
343 DO ispin = 1, nspins
344 DO iat = 1, natom
345 NULLIFY (mixing_store%cpc_h_old(iat, ispin)%r_coef)
346 NULLIFY (mixing_store%cpc_s_old(iat, ispin)%r_coef)
347 END DO
348 END DO
349 END IF
350 IF (.NOT. ASSOCIATED(mixing_store%dcpc_h_in)) THEN
351 ALLOCATE (mixing_store%dcpc_h_in(nbuffer, natom, nspins))
352 ALLOCATE (mixing_store%dcpc_s_in(nbuffer, natom, nspins))
353 ALLOCATE (mixing_store%cpc_h_lastres(natom, nspins))
354 ALLOCATE (mixing_store%cpc_s_lastres(natom, nspins))
355 DO ispin = 1, nspins
356 DO iat = 1, natom
357 DO i = 1, nbuffer
358 NULLIFY (mixing_store%dcpc_h_in(i, iat, ispin)%r_coef)
359 NULLIFY (mixing_store%dcpc_s_in(i, iat, ispin)%r_coef)
360 END DO
361 NULLIFY (mixing_store%cpc_h_lastres(iat, ispin)%r_coef)
362 NULLIFY (mixing_store%cpc_s_lastres(iat, ispin)%r_coef)
363 END DO
364 END DO
365 END IF
366 END IF
367 END IF
368 END IF
369
370 ! *** allocate multisecant buffer ***
371 IF (mixing_method == multisecant_mixing_nr) THEN
372 IF (.NOT. ASSOCIATED(mixing_store%norm_res_buffer)) THEN
373 ALLOCATE (mixing_store%norm_res_buffer(nbuffer, nspins))
374 END IF
375 END IF
376
377 IF (mixing_method == multisecant_mixing_nr) THEN
378 IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
379 ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
380 DO ispin = 1, nspins
381 DO i = 1, nbuffer
382 NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
383 END DO
384 END DO
385 END IF
386 END IF
387
388 END IF
389
390 CALL timestop(handle)
391
392 END SUBROUTINE mixing_allocate
393
394! **************************************************************************************************
395!> \brief initialiation needed when gspace mixing is used
396!> \param mixing_method ...
397!> \param rho ...
398!> \param mixing_store ...
399!> \param para_env ...
400!> \param rho_atom ...
401!> \par History
402!> 05.2009 created [MI]
403!> \author MI
404! **************************************************************************************************
405 SUBROUTINE mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
406 INTEGER, INTENT(IN) :: mixing_method
407 TYPE(qs_rho_type), POINTER :: rho
408 TYPE(mixing_storage_type), POINTER :: mixing_store
409 TYPE(mp_para_env_type), POINTER :: para_env
410 TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
411 POINTER :: rho_atom
412
413 CHARACTER(len=*), PARAMETER :: routinen = 'mixing_init'
414
415 INTEGER :: handle, iat, ib, icoef1, icoef2, ig, ig1, ig_count, iproc, ispin, n1, n2, natom, &
416 nbuffer, ncoef_h1, ncoef_h2, ncoef_s1, ncoef_s2, ng, nimg, nspin
417 REAL(dp) :: bconst, beta, cpc_h_tmp, cpc_s_tmp, &
418 fdamp, g2max, g2min, kmin
419 REAL(dp), DIMENSION(:), POINTER :: g2
420 REAL(dp), DIMENSION(:, :), POINTER :: g_vec
421 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
422 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
423
424 CALL timeset(routinen, handle)
425
426 NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
427 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
428
429 nspin = SIZE(rho_g)
430 ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
431 nimg = SIZE(rho_ao_kp, 2)
432 mixing_store%ig_max = ng
433 g2 => rho_g(1)%pw_grid%gsq
434 g_vec => rho_g(1)%pw_grid%g
435
436 IF (mixing_store%max_gvec_exp > 0._dp) THEN
437 DO ig = 1, ng
438 IF (g2(ig) > mixing_store%max_g2) THEN
439 mixing_store%ig_max = ig
440 EXIT
441 END IF
442 END DO
443 END IF
444
445 IF (.NOT. ASSOCIATED(mixing_store%kerker_factor)) THEN
446 ALLOCATE (mixing_store%kerker_factor(ng))
447 END IF
448 IF (.NOT. ASSOCIATED(mixing_store%kerker_factor_mag)) THEN
449 ALLOCATE (mixing_store%kerker_factor_mag(ng))
450 END IF
451 IF (.NOT. ASSOCIATED(mixing_store%special_metric)) THEN
452 ALLOCATE (mixing_store%special_metric(ng))
453 END IF
454 beta = mixing_store%beta
455 kmin = 0.1_dp
456 mixing_store%kerker_factor = 1.0_dp
457 mixing_store%kerker_factor_mag = 1.0_dp
458 mixing_store%special_metric = 1.0_dp
459 ig1 = 1
460 IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
461 DO ig = ig1, mixing_store%ig_max
462 mixing_store%kerker_factor(ig) = max(g2(ig)/(g2(ig) + beta*beta), kmin)
463 IF (mixing_store%beta_mag > 1.0e-10_dp) THEN
464 mixing_store%kerker_factor_mag(ig) = max(g2(ig)/(g2(ig) + mixing_store%beta_mag*mixing_store%beta_mag), kmin)
465 ELSE
466 ! beta_mag ~ 0 means no Kerker screening on the magnetization channel
467 mixing_store%kerker_factor_mag(ig) = 1.0_dp
468 END IF
469 mixing_store%special_metric(ig) = &
470 1.0_dp + 50.0_dp/8.0_dp*( &
471 1.0_dp + cos(g_vec(1, ig)) + cos(g_vec(2, ig)) + cos(g_vec(3, ig)) + &
472 cos(g_vec(1, ig))*cos(g_vec(2, ig)) + &
473 cos(g_vec(2, ig))*cos(g_vec(3, ig)) + &
474 cos(g_vec(1, ig))*cos(g_vec(3, ig)) + &
475 cos(g_vec(1, ig))*cos(g_vec(2, ig))*cos(g_vec(3, ig)))
476 END DO
477
478 nbuffer = mixing_store%nbuffer
479 DO ispin = 1, nspin
480 IF (.NOT. ASSOCIATED(mixing_store%rhoin(ispin)%cc)) THEN
481 ALLOCATE (mixing_store%rhoin(ispin)%cc(ng))
482 END IF
483 mixing_store%rhoin(ispin)%cc = rho_g(ispin)%array
484
485 IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
486 IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer(1, ispin)%cc)) THEN
487 DO ib = 1, nbuffer
488 ALLOCATE (mixing_store%rhoin_buffer(ib, ispin)%cc(ng))
489 END DO
490 END IF
491 mixing_store%rhoin_buffer(1, ispin)%cc(1:ng) = &
492 rho_g(ispin)%array(1:ng)
493 END IF
494 IF (ASSOCIATED(mixing_store%res_buffer)) THEN
495 IF (.NOT. ASSOCIATED(mixing_store%res_buffer(1, ispin)%cc)) THEN
496 DO ib = 1, nbuffer
497 ALLOCATE (mixing_store%res_buffer(ib, ispin)%cc(ng))
498 END DO
499 END IF
500 END IF
501 END DO
502
503 IF (nspin == 2) THEN
504 mixing_store%rhoin(1)%cc = rho_g(1)%array + rho_g(2)%array
505 mixing_store%rhoin(2)%cc = rho_g(1)%array - rho_g(2)%array
506 IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
507 mixing_store%rhoin_buffer(1, 1)%cc = rho_g(1)%array + rho_g(2)%array
508 mixing_store%rhoin_buffer(1, 2)%cc = rho_g(1)%array - rho_g(2)%array
509 END IF
510 END IF
511
512 IF (mixing_store%gmix_p) THEN
513 IF (PRESENT(rho_atom)) THEN
514 natom = SIZE(rho_atom)
515 DO ispin = 1, nspin
516 DO iat = 1, natom
517 IF (ASSOCIATED(rho_atom(iat)%cpc_s(ispin)%r_coef)) THEN
518 mixing_store%paw(iat) = .true.
519 n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
520 n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
521 IF (ASSOCIATED(mixing_store%cpc_s_in)) THEN
522 IF (.NOT. ASSOCIATED(mixing_store%cpc_s_in(iat, ispin)%r_coef)) THEN
523 ALLOCATE (mixing_store%cpc_s_in(iat, ispin)%r_coef(n1, n2))
524 ALLOCATE (mixing_store%cpc_h_in(iat, ispin)%r_coef(n1, n2))
525 END IF
526 mixing_store%cpc_h_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_h(ispin)%r_coef
527 mixing_store%cpc_s_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_s(ispin)%r_coef
528 END IF
529 END IF
530 END DO
531 END DO
532 IF (nspin == 2 .AND. ASSOCIATED(mixing_store%cpc_s_in)) THEN
533 DO iat = 1, natom
534 IF (mixing_store%paw(iat)) THEN
535 cpassert(ASSOCIATED(mixing_store%cpc_h_in(iat, 1)%r_coef))
536 cpassert(ASSOCIATED(mixing_store%cpc_h_in(iat, 2)%r_coef))
537 cpassert(ASSOCIATED(mixing_store%cpc_s_in(iat, 1)%r_coef))
538 cpassert(ASSOCIATED(mixing_store%cpc_s_in(iat, 2)%r_coef))
539
540 ncoef_h1 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 1)
541 ncoef_h2 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 2)
542 ncoef_s1 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 1)
543 ncoef_s2 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 2)
544 cpassert(ncoef_h1 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 1))
545 cpassert(ncoef_h2 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 2))
546 cpassert(ncoef_s1 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 1))
547 cpassert(ncoef_s2 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 2))
548
549 DO icoef2 = 1, ncoef_h2
550 DO icoef1 = 1, ncoef_h1
551 cpc_h_tmp = mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2)
552 mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2) = &
553 cpc_h_tmp + mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
554 mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2) = &
555 cpc_h_tmp - mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
556 END DO
557 END DO
558
559 DO icoef2 = 1, ncoef_s2
560 DO icoef1 = 1, ncoef_s1
561 cpc_s_tmp = mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2)
562 mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2) = &
563 cpc_s_tmp + mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
564 mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2) = &
565 cpc_s_tmp - mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
566 END DO
567 END DO
568 END IF
569 END DO
570 END IF
571 END IF
572 END IF
573
574 IF (mixing_method == gspace_mixing_nr) THEN
575 ELSEIF (mixing_method == pulay_mixing_nr) THEN
576 IF (mixing_store%gmix_p .AND. PRESENT(rho_atom)) THEN
577 DO ispin = 1, nspin
578 DO iat = 1, natom
579 IF (mixing_store%paw(iat)) THEN
580 n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
581 n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
582 IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer(1, iat, ispin)%r_coef)) THEN
583 DO ib = 1, nbuffer
584 ALLOCATE (mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
585 ALLOCATE (mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
586 ALLOCATE (mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
587 ALLOCATE (mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
588 END DO
589 END IF
590 DO ib = 1, nbuffer
591 mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
592 mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
593 mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
594 mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
595 END DO
596 END IF
597 END DO
598 END DO
599 END IF
600 ELSEIF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
601 DO ispin = 1, nspin
602 IF (.NOT. ASSOCIATED(mixing_store%rhoin_old(ispin)%cc)) THEN
603 ALLOCATE (mixing_store%rhoin_old(ispin)%cc(ng))
604 END IF
605 IF (.NOT. ASSOCIATED(mixing_store%drho_buffer(1, ispin)%cc)) THEN
606 DO ib = 1, nbuffer
607 ALLOCATE (mixing_store%drho_buffer(ib, ispin)%cc(ng))
608 END DO
609 ALLOCATE (mixing_store%last_res(ispin)%cc(ng))
610 END IF
611 DO ib = 1, nbuffer
612 mixing_store%drho_buffer(ib, ispin)%cc = z_zero
613 END DO
614 mixing_store%last_res(ispin)%cc = z_zero
615 mixing_store%rhoin_old(ispin)%cc = z_zero
616 END DO
617 IF (mixing_store%gmix_p) THEN
618 IF (PRESENT(rho_atom)) THEN
619 DO ispin = 1, nspin
620 DO iat = 1, natom
621 IF (mixing_store%paw(iat)) THEN
622 n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
623 n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
624 IF (.NOT. ASSOCIATED(mixing_store%cpc_s_old(iat, ispin)%r_coef)) THEN
625 ALLOCATE (mixing_store%cpc_s_old(iat, ispin)%r_coef(n1, n2))
626 ALLOCATE (mixing_store%cpc_h_old(iat, ispin)%r_coef(n1, n2))
627 END IF
628 mixing_store%cpc_h_old(iat, ispin)%r_coef = 0.0_dp
629 mixing_store%cpc_s_old(iat, ispin)%r_coef = 0.0_dp
630 IF (.NOT. ASSOCIATED(mixing_store%dcpc_s_in(1, iat, ispin)%r_coef)) THEN
631 DO ib = 1, nbuffer
632 ALLOCATE (mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef(n1, n2))
633 ALLOCATE (mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef(n1, n2))
634 END DO
635 ALLOCATE (mixing_store%cpc_h_lastres(iat, ispin)%r_coef(n1, n2))
636 ALLOCATE (mixing_store%cpc_s_lastres(iat, ispin)%r_coef(n1, n2))
637 END IF
638 DO ib = 1, nbuffer
639 mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef = 0.0_dp
640 mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef = 0.0_dp
641 END DO
642 mixing_store%cpc_h_lastres(iat, ispin)%r_coef = 0.0_dp
643 mixing_store%cpc_s_lastres(iat, ispin)%r_coef = 0.0_dp
644 END IF
645 END DO
646 END DO
647 END IF
648 END IF
649
650 IF (.NOT. ASSOCIATED(mixing_store%p_metric)) THEN
651 ALLOCATE (mixing_store%p_metric(ng))
652 bconst = mixing_store%bconst
653 g2min = 1.e30_dp
654 DO ig = 1, ng
655 IF (g2(ig) > 1.e-10_dp) g2min = min(g2min, g2(ig))
656 END DO
657 g2max = -1.e30_dp
658 DO ig = 1, ng
659 g2max = max(g2max, g2(ig))
660 END DO
661 CALL para_env%min(g2min)
662 CALL para_env%max(g2max)
663 ! fdamp/g2 varies between (bconst-1) and 0
664 ! i.e. p_metric varies between bconst and 1
665 ! fdamp = (bconst-1.0_dp)*g2min
666 fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
667 DO ig = 1, ng
668 mixing_store%p_metric(ig) = (g2(ig) + fdamp)/max(g2(ig), 1.e-10_dp)
669 END DO
670 IF (rho_g(1)%pw_grid%have_g0) mixing_store%p_metric(1) = bconst
671 END IF
672 ELSEIF (mixing_method == multisecant_mixing_nr) THEN
673 IF (.NOT. ASSOCIATED(mixing_store%ig_global_index)) THEN
674 ALLOCATE (mixing_store%ig_global_index(ng))
675 END IF
676 mixing_store%ig_global_index = 0
677 ig_count = 0
678 DO iproc = 0, para_env%num_pe - 1
679 IF (para_env%mepos == iproc) THEN
680 DO ig = 1, ng
681 ig_count = ig_count + 1
682 mixing_store%ig_global_index(ig) = ig_count
683 END DO
684 END IF
685 CALL para_env%bcast(ig_count, iproc)
686 END DO
687 END IF
688
689 CALL timestop(handle)
690
691 END SUBROUTINE mixing_init
692
693! **************************************************************************************************
694!> \brief initialiation needed when charge mixing is used
695!> \param mixing_store ...
696!> \par History
697!> 02.2019 created [JGH]
698!> \author JGH
699! **************************************************************************************************
700 ELEMENTAL SUBROUTINE charge_mixing_init(mixing_store)
701 TYPE(mixing_storage_type), INTENT(INOUT) :: mixing_store
702
703 mixing_store%ncall = 0
704 mixing_store%tb_scc_mixer_error = 0.0_dp
705 mixing_store%tb_scc_mixer_step = 0
706 IF (ASSOCIATED(mixing_store%acharge)) mixing_store%acharge = 0.0_dp
707 IF (ASSOCIATED(mixing_store%dacharge)) mixing_store%dacharge = 0.0_dp
708
709 END SUBROUTINE charge_mixing_init
710
711END MODULE qs_mixing_utils
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
pw_types
Definition pw_types.F:24
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::broyden_mixing_nr
integer, parameter, public broyden_mixing_nr
Definition qs_density_mixing_types.F:49
qs_density_mixing_types::modified_broyden_mixing_nr
integer, parameter, public modified_broyden_mixing_nr
Definition qs_density_mixing_types.F:49
qs_density_mixing_types::multisecant_mixing_nr
integer, parameter, public multisecant_mixing_nr
Definition qs_density_mixing_types.F:49
qs_density_mixing_types::pulay_mixing_nr
integer, parameter, public pulay_mixing_nr
Definition qs_density_mixing_types.F:49
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition qs_density_mixing_types.F:49
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_mixing_utils
Definition qs_mixing_utils.F:9
qs_mixing_utils::charge_mixing_init
elemental subroutine, public charge_mixing_init(mixing_store)
initialiation needed when charge mixing is used
Definition qs_mixing_utils.F:701
qs_mixing_utils::mixing_init
subroutine, public mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
initialiation needed when gspace mixing is used
Definition qs_mixing_utils.F:406
qs_mixing_utils::mixing_allocate
subroutine, public mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
allocation needed when density mixing is used
Definition qs_mixing_utils.F:113
qs_mixing_utils::self_consistency_check
subroutine, public self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
...
Definition qs_mixing_utils.F:60
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition qs_scf_methods.F:19
qs_scf_methods::cp_sm_mix
subroutine, public cp_sm_mix(m1, m2, p_mix, delta, para_env, m3)
Perform a mixing of the given matrixes into the first matrix m1 = m2 + p_mix (m1-m2)
Definition qs_scf_methods.F:547
cp_control_types::dft_control_type
Definition cp_control_types.F:726
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
qs_density_mixing_types::mixing_storage_type
Definition qs_density_mixing_types.F:60
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62