(git:31876b5)
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 bibliography, ONLY: sundararaman2017,&
11 cite_reference
12 USE cp_control_types, ONLY: dft_control_type
13 USE cp_dbcsr_api, ONLY: dbcsr_create,&
14 dbcsr_p_type,&
15 dbcsr_set,&
16 dbcsr_type,&
17 dbcsr_type_symmetric
18 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
19 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
20 USE distribution_1d_types, ONLY: distribution_1d_type
21 USE kinds, ONLY: dp
22 USE mathconstants, ONLY: z_zero
23 USE message_passing, ONLY: mp_para_env_type
24 USE pw_types, ONLY: pw_c1d_gs_type
25 USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
26 gspace_mixing_nr,&
27 mixing_storage_type,&
28 modified_broyden_mixing_nr,&
29 multisecant_mixing_nr,&
30 new_pulay_mixing_nr,&
31 pulay_mixing_nr
32 USE qs_environment_types, ONLY: get_qs_env,&
33 qs_environment_type
34 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
35 USE qs_rho_atom_types, ONLY: rho_atom_type
36 USE qs_rho_types, ONLY: qs_rho_get,&
37 qs_rho_type
38 USE qs_scf_methods, ONLY: cp_sm_mix
39#include "./base/base_uses.f90"
40
41 IMPLICIT NONE
42
43 PRIVATE
44
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mixing_utils'
46 INTEGER, PARAMETER, PRIVATE :: max_dftb_scc_vars = 1, &
47 max_xtb_scc_vars = 28
48
49 PUBLIC :: mixing_allocate, mixing_init, charge_mixing_init, &
50 self_consistency_check
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief ...
56!> \param rho_ao ...
57!> \param p_delta ...
58!> \param para_env ...
59!> \param p_out ...
60!> \param delta ...
61! **************************************************************************************************
62 SUBROUTINE self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
63 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao, p_delta
64 TYPE(mp_para_env_type), POINTER :: para_env
65 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_out
66 REAL(kind=dp), INTENT(INOUT) :: delta
67
68 CHARACTER(len=*), PARAMETER :: routinen = 'self_consistency_check'
69
70 INTEGER :: handle, ic, ispin, nimg, nspins
71 REAL(kind=dp) :: tmp
72 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_q, p_in
73
74 CALL timeset(routinen, handle)
75
76 NULLIFY (matrix_q, p_in)
77
78 cpassert(ASSOCIATED(p_out))
79 NULLIFY (matrix_q, p_in)
80 p_in => rho_ao
81 matrix_q => p_delta
82 nspins = SIZE(p_in, 1)
83 nimg = SIZE(p_in, 2)
84
85 ! Compute the difference (p_out - p_in)and check convergence
86 delta = 0.0_dp
87 DO ispin = 1, nspins
88 DO ic = 1, nimg
89 CALL dbcsr_set(matrix_q(ispin, ic)%matrix, 0.0_dp)
90 CALL cp_sm_mix(m1=p_out(ispin, ic)%matrix, m2=p_in(ispin, ic)%matrix, &
91 p_mix=1.0_dp, delta=tmp, para_env=para_env, &
92 m3=matrix_q(ispin, ic)%matrix)
93 delta = max(tmp, delta)
94 END DO
95 END DO
96 CALL timestop(handle)
97
98 END SUBROUTINE self_consistency_check
99
100! **************************************************************************************************
101!> \brief allocation needed when density mixing is used
102!> \param qs_env ...
103!> \param mixing_method ...
104!> \param p_mix_new ...
105!> \param p_delta ...
106!> \param nspins ...
107!> \param mixing_store ...
108!> \par History
109!> 05.2009 created [MI]
110!> 08.2014 kpoints [JGH]
111!> 02.2015 changed input to make g-space mixing available in linear scaling SCF [Patrick Seewald]
112!> 02.2019 charge mixing [JGH]
113!> \author fawzi
114! **************************************************************************************************
115 SUBROUTINE mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
116 TYPE(qs_environment_type), POINTER :: qs_env
117 INTEGER :: mixing_method
118 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
119 POINTER :: p_mix_new, p_delta
120 INTEGER, INTENT(IN) :: nspins
121 TYPE(mixing_storage_type), POINTER :: mixing_store
122
123 CHARACTER(LEN=*), PARAMETER :: routinen = 'mixing_allocate'
124
125 INTEGER :: handle, i, ia, iat, ic, ikind, ispin, &
126 max_shell, na, natom, nbuffer, nel, &
127 nimg, nkind
128 LOGICAL :: charge_mixing
129 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
130 TYPE(dbcsr_type), POINTER :: refmatrix
131 TYPE(dft_control_type), POINTER :: dft_control
132 TYPE(distribution_1d_type), POINTER :: distribution_1d
133 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
134 POINTER :: sab_orb
135 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
136
137 CALL timeset(routinen, handle)
138
139 NULLIFY (matrix_s, dft_control, sab_orb, refmatrix, rho_atom)
140 CALL get_qs_env(qs_env, &
141 sab_orb=sab_orb, &
142 matrix_s_kp=matrix_s, &
143 dft_control=dft_control)
144
145 charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
146 .OR. dft_control%qs_control%semi_empirical
147 refmatrix => matrix_s(1, 1)%matrix
148 nimg = dft_control%nimages
149
150 ! *** allocate p_mix_new ***
151 IF (PRESENT(p_mix_new)) THEN
152 IF (.NOT. ASSOCIATED(p_mix_new)) THEN
153 CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
154 DO i = 1, nspins
155 DO ic = 1, nimg
156 ALLOCATE (p_mix_new(i, ic)%matrix)
157 CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
158 name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
159 CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
160 CALL dbcsr_set(p_mix_new(i, ic)%matrix, 0.0_dp)
161 END DO
162 END DO
163 END IF
164 END IF
165
166 ! *** allocate p_delta ***
167 IF (PRESENT(p_delta)) THEN
168 IF (mixing_method >= gspace_mixing_nr) THEN
169 IF (.NOT. ASSOCIATED(p_delta)) THEN
170 CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
171 DO i = 1, nspins
172 DO ic = 1, nimg
173 ALLOCATE (p_delta(i, ic)%matrix)
174 CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
175 name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
176 CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
177 CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
178 END DO
179 END DO
180 END IF
181 cpassert(ASSOCIATED(mixing_store))
182 END IF
183 END IF
184
185 IF (charge_mixing) THEN
186 ! *** allocate buffer for TB SCC variable mixing ***
187 IF (mixing_method >= gspace_mixing_nr) THEN
188 cpassert(.NOT. mixing_store%gmix_p)
189 IF (dft_control%qs_control%dftb) THEN
190 max_shell = max_dftb_scc_vars
191 ELSEIF (dft_control%qs_control%xtb) THEN
192 max_shell = max_xtb_scc_vars
193 ELSE
194 cpabort('UNKNOWN METHOD')
195 END IF
196 nbuffer = mixing_store%nbuffer
197 mixing_store%ncall = 0
198 CALL get_qs_env(qs_env, local_particles=distribution_1d)
199 nkind = SIZE(distribution_1d%n_el)
200 na = sum(distribution_1d%n_el(:))
201 IF (ASSOCIATED(mixing_store%atlist)) DEALLOCATE (mixing_store%atlist)
202 ALLOCATE (mixing_store%atlist(na))
203 mixing_store%nat_local = na
204 mixing_store%max_shell = max_shell
205 ia = 0
206 DO ikind = 1, nkind
207 nel = distribution_1d%n_el(ikind)
208 DO iat = 1, nel
209 ia = ia + 1
210 mixing_store%atlist(ia) = distribution_1d%list(ikind)%array(iat)
211 END DO
212 END DO
213 IF (ASSOCIATED(mixing_store%acharge)) DEALLOCATE (mixing_store%acharge)
214 ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer))
215 IF (ASSOCIATED(mixing_store%dacharge)) DEALLOCATE (mixing_store%dacharge)
216 ALLOCATE (mixing_store%dacharge(na, max_shell, nbuffer))
217 END IF
218 IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
219 IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
220 ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
221 mixing_store%pulay_matrix = 0.0_dp
222 ELSEIF (mixing_method == broyden_mixing_nr) THEN
223 IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
224 ALLOCATE (mixing_store%abroy(nbuffer, nbuffer))
225 mixing_store%abroy = 0.0_dp
226 IF (ASSOCIATED(mixing_store%wbroy)) DEALLOCATE (mixing_store%wbroy)
227 ALLOCATE (mixing_store%wbroy(nbuffer))
228 mixing_store%wbroy = 0.0_dp
229 IF (ASSOCIATED(mixing_store%dfbroy)) DEALLOCATE (mixing_store%dfbroy)
230 ALLOCATE (mixing_store%dfbroy(na, max_shell, nbuffer))
231 mixing_store%dfbroy = 0.0_dp
232 IF (ASSOCIATED(mixing_store%ubroy)) DEALLOCATE (mixing_store%ubroy)
233 ALLOCATE (mixing_store%ubroy(na, max_shell, nbuffer))
234 mixing_store%ubroy = 0.0_dp
235 ELSEIF (mixing_method == multisecant_mixing_nr) THEN
236 cpabort("multisecant_mixing not available")
237 END IF
238 ELSE
239 ! *** allocate buffer for gspace mixing ***
240 IF (mixing_method >= gspace_mixing_nr) THEN
241 nbuffer = mixing_store%nbuffer
242 mixing_store%ncall = 0
243 IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
244 ALLOCATE (mixing_store%rhoin(nspins))
245 DO ispin = 1, nspins
246 NULLIFY (mixing_store%rhoin(ispin)%cc)
247 END DO
248 END IF
249
250 IF (mixing_store%gmix_p .AND. dft_control%qs_control%gapw) THEN
251 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
252 natom = SIZE(rho_atom)
253 IF (.NOT. ASSOCIATED(mixing_store%paw)) THEN
254 ALLOCATE (mixing_store%paw(natom))
255 mixing_store%paw = .false.
256 ALLOCATE (mixing_store%cpc_h_in(natom, nspins))
257 ALLOCATE (mixing_store%cpc_s_in(natom, nspins))
258 DO ispin = 1, nspins
259 DO iat = 1, natom
260 NULLIFY (mixing_store%cpc_h_in(iat, ispin)%r_coef)
261 NULLIFY (mixing_store%cpc_s_in(iat, ispin)%r_coef)
262 END DO
263 END DO
264 END IF
265 END IF
266 END IF
267
268 ! *** allocare res_buffer if needed
269 IF (mixing_method >= pulay_mixing_nr) THEN
270 IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
271 ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
272 DO ispin = 1, nspins
273 DO i = 1, nbuffer
274 NULLIFY (mixing_store%res_buffer(i, ispin)%cc)
275 END DO
276 END DO
277 END IF
278 END IF
279
280 ! *** allocate pulay
281 IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
282 IF (.NOT. ASSOCIATED(mixing_store%pulay_matrix)) THEN
283 ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
284 END IF
285
286 IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
287 ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
288 DO ispin = 1, nspins
289 DO i = 1, nbuffer
290 NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
291 END DO
292 END DO
293 END IF
294 IF (mixing_store%gmix_p) THEN
295 IF (dft_control%qs_control%gapw) THEN
296 IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer)) THEN
297 ALLOCATE (mixing_store%cpc_h_in_buffer(nbuffer, natom, nspins))
298 ALLOCATE (mixing_store%cpc_s_in_buffer(nbuffer, natom, nspins))
299 ALLOCATE (mixing_store%cpc_h_res_buffer(nbuffer, natom, nspins))
300 ALLOCATE (mixing_store%cpc_s_res_buffer(nbuffer, natom, nspins))
301 DO ispin = 1, nspins
302 DO iat = 1, natom
303 DO i = 1, nbuffer
304 NULLIFY (mixing_store%cpc_h_in_buffer(i, iat, ispin)%r_coef)
305 NULLIFY (mixing_store%cpc_s_in_buffer(i, iat, ispin)%r_coef)
306 NULLIFY (mixing_store%cpc_h_res_buffer(i, iat, ispin)%r_coef)
307 NULLIFY (mixing_store%cpc_s_res_buffer(i, iat, ispin)%r_coef)
308 END DO
309 END DO
310 END DO
311 END IF
312 END IF
313 END IF
314
315 END IF
316 ! *** allocate broyden buffer ***
317 IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
318 IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
319 ALLOCATE (mixing_store%rhoin_old(nspins))
320 DO ispin = 1, nspins
321 NULLIFY (mixing_store%rhoin_old(ispin)%cc)
322 END DO
323 END IF
324 IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
325 ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
326 ALLOCATE (mixing_store%last_res(nspins))
327 DO ispin = 1, nspins
328 DO i = 1, nbuffer
329 NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
330 END DO
331 NULLIFY (mixing_store%last_res(ispin)%cc)
332 END DO
333 END IF
334 IF (mixing_method == modified_broyden_mixing_nr) THEN
335 IF (.NOT. ASSOCIATED(mixing_store%weight)) THEN
336 ALLOCATE (mixing_store%weight(nbuffer, nspins))
337 END IF
338 mixing_store%weight = 0.0_dp
339 END IF
340 IF (mixing_store%gmix_p) THEN
341
342 IF (dft_control%qs_control%gapw) THEN
343 IF (.NOT. ASSOCIATED(mixing_store%cpc_h_old)) THEN
344 ALLOCATE (mixing_store%cpc_h_old(natom, nspins))
345 ALLOCATE (mixing_store%cpc_s_old(natom, nspins))
346 DO ispin = 1, nspins
347 DO iat = 1, natom
348 NULLIFY (mixing_store%cpc_h_old(iat, ispin)%r_coef)
349 NULLIFY (mixing_store%cpc_s_old(iat, ispin)%r_coef)
350 END DO
351 END DO
352 END IF
353 IF (.NOT. ASSOCIATED(mixing_store%dcpc_h_in)) THEN
354 ALLOCATE (mixing_store%dcpc_h_in(nbuffer, natom, nspins))
355 ALLOCATE (mixing_store%dcpc_s_in(nbuffer, natom, nspins))
356 ALLOCATE (mixing_store%cpc_h_lastres(natom, nspins))
357 ALLOCATE (mixing_store%cpc_s_lastres(natom, nspins))
358 DO ispin = 1, nspins
359 DO iat = 1, natom
360 DO i = 1, nbuffer
361 NULLIFY (mixing_store%dcpc_h_in(i, iat, ispin)%r_coef)
362 NULLIFY (mixing_store%dcpc_s_in(i, iat, ispin)%r_coef)
363 END DO
364 NULLIFY (mixing_store%cpc_h_lastres(iat, ispin)%r_coef)
365 NULLIFY (mixing_store%cpc_s_lastres(iat, ispin)%r_coef)
366 END DO
367 END DO
368 END IF
369 END IF
370 END IF
371 END IF
372
373 ! *** allocate multisecant buffer ***
374 IF (mixing_method == multisecant_mixing_nr) THEN
375 IF (.NOT. ASSOCIATED(mixing_store%norm_res_buffer)) THEN
376 ALLOCATE (mixing_store%norm_res_buffer(nbuffer, nspins))
377 END IF
378 END IF
379
380 IF (mixing_method == multisecant_mixing_nr) THEN
381 IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
382 ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
383 DO ispin = 1, nspins
384 DO i = 1, nbuffer
385 NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
386 END DO
387 END DO
388 END IF
389 END IF
390
391 END IF
392
393 CALL timestop(handle)
394
395 END SUBROUTINE mixing_allocate
396
397! **************************************************************************************************
398!> \brief initialiation needed when gspace mixing is used
399!> \param mixing_method ...
400!> \param rho ...
401!> \param mixing_store ...
402!> \param para_env ...
403!> \param rho_atom ...
404!> \par History
405!> 05.2009 created [MI]
406!> \author MI
407! **************************************************************************************************
408 SUBROUTINE mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
409 INTEGER, INTENT(IN) :: mixing_method
410 TYPE(qs_rho_type), POINTER :: rho
411 TYPE(mixing_storage_type), POINTER :: mixing_store
412 TYPE(mp_para_env_type), POINTER :: para_env
413 TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
414 POINTER :: rho_atom
415
416 CHARACTER(len=*), PARAMETER :: routinen = 'mixing_init'
417
418 INTEGER :: handle, iat, ib, icoef1, icoef2, ig, ig1, ig_count, iproc, ispin, n1, n2, natom, &
419 nbuffer, ncoef_h1, ncoef_h2, ncoef_s1, ncoef_s2, ng, nimg, nspin
420 REAL(dp) :: bconst, beta, cpc_h_tmp, cpc_s_tmp, &
421 fdamp, g2max, g2min, kmin, qk, qkappa, &
422 qm
423 REAL(dp), DIMENSION(:), POINTER :: g2
424 REAL(dp), DIMENSION(:, :), POINTER :: g_vec
425 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
426 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
427
428 CALL timeset(routinen, handle)
429
430 NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
431 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
432
433 nspin = SIZE(rho_g)
434 ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
435 nimg = SIZE(rho_ao_kp, 2)
436 mixing_store%ig_max = ng
437 g2 => rho_g(1)%pw_grid%gsq
438 g_vec => rho_g(1)%pw_grid%g
439
440 IF (mixing_store%max_gvec_exp > 0._dp) THEN
441 DO ig = 1, ng
442 IF (g2(ig) > mixing_store%max_g2) THEN
443 mixing_store%ig_max = ig
444 EXIT
445 END IF
446 END DO
447 END IF
448
449 IF (.NOT. ASSOCIATED(mixing_store%kerker_factor)) THEN
450 ALLOCATE (mixing_store%kerker_factor(ng))
451 END IF
452 IF (.NOT. ASSOCIATED(mixing_store%kerker_factor_mag)) THEN
453 ALLOCATE (mixing_store%kerker_factor_mag(ng))
454 END IF
455 IF (.NOT. ASSOCIATED(mixing_store%special_metric)) THEN
456 ALLOCATE (mixing_store%special_metric(ng))
457 END IF
458 beta = mixing_store%beta
459 kmin = 0.1_dp
460 mixing_store%kerker_factor = 1.0_dp
461 mixing_store%kerker_factor_mag = 1.0_dp
462 mixing_store%special_metric = 1.0_dp
463 ig1 = 1
464
465 IF (.NOT. (mixing_method == new_pulay_mixing_nr)) THEN
466 IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
467 DO ig = ig1, mixing_store%ig_max
468 mixing_store%kerker_factor(ig) = max(g2(ig)/(g2(ig) + beta*beta), kmin)
469 IF (mixing_store%beta_mag > 1.0e-10_dp) THEN
470 mixing_store%kerker_factor_mag(ig) = max(g2(ig)/(g2(ig) + mixing_store%beta_mag*mixing_store%beta_mag), kmin)
471 ELSE
472 ! beta_mag ~ 0 means no Kerker screening on the magnetization channel
473 mixing_store%kerker_factor_mag(ig) = 1.0_dp
474 END IF
475 mixing_store%special_metric(ig) = &
476 1.0_dp + 50.0_dp/8.0_dp*( &
477 1.0_dp + cos(g_vec(1, ig)) + cos(g_vec(2, ig)) + cos(g_vec(3, ig)) + &
478 cos(g_vec(1, ig))*cos(g_vec(2, ig)) + &
479 cos(g_vec(2, ig))*cos(g_vec(3, ig)) + &
480 cos(g_vec(1, ig))*cos(g_vec(3, ig)) + &
481 cos(g_vec(1, ig))*cos(g_vec(2, ig))*cos(g_vec(3, ig)))
482 END DO
483 ELSE
484 CALL cite_reference(sundararaman2017)
485 qkappa = mixing_store%qkappa
486 qk = mixing_store%qk
487 qm = mixing_store%qm
488 IF (rho_g(1)%pw_grid%have_g0) ig1 = 1
489 DO ig = ig1, mixing_store%ig_max
490 mixing_store%kerker_factor(ig) = (g2(ig) + qkappa*qkappa)/(g2(ig) + qkappa*qkappa + qk*qk)
491 mixing_store%special_metric(ig) = (g2(ig) + qkappa*qkappa + qm*qm)/(g2(ig) + qkappa*qkappa)
492 END DO
493 END IF
494 nbuffer = mixing_store%nbuffer
495 DO ispin = 1, nspin
496 IF (.NOT. ASSOCIATED(mixing_store%rhoin(ispin)%cc)) THEN
497 ALLOCATE (mixing_store%rhoin(ispin)%cc(ng))
498 END IF
499 mixing_store%rhoin(ispin)%cc = rho_g(ispin)%array
500
501 IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
502 IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer(1, ispin)%cc)) THEN
503 DO ib = 1, nbuffer
504 ALLOCATE (mixing_store%rhoin_buffer(ib, ispin)%cc(ng))
505 END DO
506 END IF
507 mixing_store%rhoin_buffer(1, ispin)%cc(1:ng) = &
508 rho_g(ispin)%array(1:ng)
509 END IF
510 IF (ASSOCIATED(mixing_store%res_buffer)) THEN
511 IF (.NOT. ASSOCIATED(mixing_store%res_buffer(1, ispin)%cc)) THEN
512 DO ib = 1, nbuffer
513 ALLOCATE (mixing_store%res_buffer(ib, ispin)%cc(ng))
514 END DO
515 END IF
516 END IF
517 END DO
518
519 IF (nspin == 2) THEN
520 mixing_store%rhoin(1)%cc = rho_g(1)%array + rho_g(2)%array
521 mixing_store%rhoin(2)%cc = rho_g(1)%array - rho_g(2)%array
522 IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
523 mixing_store%rhoin_buffer(1, 1)%cc = rho_g(1)%array + rho_g(2)%array
524 mixing_store%rhoin_buffer(1, 2)%cc = rho_g(1)%array - rho_g(2)%array
525 END IF
526 END IF
527
528 IF (mixing_store%gmix_p) THEN
529 IF (PRESENT(rho_atom)) THEN
530 natom = SIZE(rho_atom)
531 DO ispin = 1, nspin
532 DO iat = 1, natom
533 IF (ASSOCIATED(rho_atom(iat)%cpc_s(ispin)%r_coef)) THEN
534 mixing_store%paw(iat) = .true.
535 n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
536 n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
537 IF (ASSOCIATED(mixing_store%cpc_s_in)) THEN
538 IF (.NOT. ASSOCIATED(mixing_store%cpc_s_in(iat, ispin)%r_coef)) THEN
539 ALLOCATE (mixing_store%cpc_s_in(iat, ispin)%r_coef(n1, n2))
540 ALLOCATE (mixing_store%cpc_h_in(iat, ispin)%r_coef(n1, n2))
541 END IF
542 mixing_store%cpc_h_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_h(ispin)%r_coef
543 mixing_store%cpc_s_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_s(ispin)%r_coef
544 END IF
545 END IF
546 END DO
547 END DO
548 IF (nspin == 2 .AND. ASSOCIATED(mixing_store%cpc_s_in)) THEN
549 DO iat = 1, natom
550 IF (mixing_store%paw(iat)) THEN
551 cpassert(ASSOCIATED(mixing_store%cpc_h_in(iat, 1)%r_coef))
552 cpassert(ASSOCIATED(mixing_store%cpc_h_in(iat, 2)%r_coef))
553 cpassert(ASSOCIATED(mixing_store%cpc_s_in(iat, 1)%r_coef))
554 cpassert(ASSOCIATED(mixing_store%cpc_s_in(iat, 2)%r_coef))
555
556 ncoef_h1 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 1)
557 ncoef_h2 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 2)
558 ncoef_s1 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 1)
559 ncoef_s2 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 2)
560 cpassert(ncoef_h1 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 1))
561 cpassert(ncoef_h2 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 2))
562 cpassert(ncoef_s1 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 1))
563 cpassert(ncoef_s2 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 2))
564
565 DO icoef2 = 1, ncoef_h2
566 DO icoef1 = 1, ncoef_h1
567 cpc_h_tmp = mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2)
568 mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2) = &
569 cpc_h_tmp + mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
570 mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2) = &
571 cpc_h_tmp - mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
572 END DO
573 END DO
574
575 DO icoef2 = 1, ncoef_s2
576 DO icoef1 = 1, ncoef_s1
577 cpc_s_tmp = mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2)
578 mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2) = &
579 cpc_s_tmp + mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
580 mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2) = &
581 cpc_s_tmp - mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
582 END DO
583 END DO
584 END IF
585 END DO
586 END IF
587 END IF
588 END IF
589
590 IF (mixing_method == gspace_mixing_nr) THEN
591 ELSEIF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
592 IF (mixing_store%gmix_p .AND. PRESENT(rho_atom)) THEN
593 DO ispin = 1, nspin
594 DO iat = 1, natom
595 IF (mixing_store%paw(iat)) THEN
596 n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
597 n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
598 IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer(1, iat, ispin)%r_coef)) THEN
599 DO ib = 1, nbuffer
600 ALLOCATE (mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
601 ALLOCATE (mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
602 ALLOCATE (mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
603 ALLOCATE (mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
604 END DO
605 END IF
606 DO ib = 1, nbuffer
607 mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
608 mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
609 mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
610 mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
611 END DO
612 END IF
613 END DO
614 END DO
615 END IF
616 ELSEIF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
617 DO ispin = 1, nspin
618 IF (.NOT. ASSOCIATED(mixing_store%rhoin_old(ispin)%cc)) THEN
619 ALLOCATE (mixing_store%rhoin_old(ispin)%cc(ng))
620 END IF
621 IF (.NOT. ASSOCIATED(mixing_store%drho_buffer(1, ispin)%cc)) THEN
622 DO ib = 1, nbuffer
623 ALLOCATE (mixing_store%drho_buffer(ib, ispin)%cc(ng))
624 END DO
625 ALLOCATE (mixing_store%last_res(ispin)%cc(ng))
626 END IF
627 DO ib = 1, nbuffer
628 mixing_store%drho_buffer(ib, ispin)%cc = z_zero
629 END DO
630 mixing_store%last_res(ispin)%cc = z_zero
631 mixing_store%rhoin_old(ispin)%cc = z_zero
632 END DO
633 IF (mixing_store%gmix_p) THEN
634 IF (PRESENT(rho_atom)) THEN
635 DO ispin = 1, nspin
636 DO iat = 1, natom
637 IF (mixing_store%paw(iat)) THEN
638 n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
639 n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
640 IF (.NOT. ASSOCIATED(mixing_store%cpc_s_old(iat, ispin)%r_coef)) THEN
641 ALLOCATE (mixing_store%cpc_s_old(iat, ispin)%r_coef(n1, n2))
642 ALLOCATE (mixing_store%cpc_h_old(iat, ispin)%r_coef(n1, n2))
643 END IF
644 mixing_store%cpc_h_old(iat, ispin)%r_coef = 0.0_dp
645 mixing_store%cpc_s_old(iat, ispin)%r_coef = 0.0_dp
646 IF (.NOT. ASSOCIATED(mixing_store%dcpc_s_in(1, iat, ispin)%r_coef)) THEN
647 DO ib = 1, nbuffer
648 ALLOCATE (mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef(n1, n2))
649 ALLOCATE (mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef(n1, n2))
650 END DO
651 ALLOCATE (mixing_store%cpc_h_lastres(iat, ispin)%r_coef(n1, n2))
652 ALLOCATE (mixing_store%cpc_s_lastres(iat, ispin)%r_coef(n1, n2))
653 END IF
654 DO ib = 1, nbuffer
655 mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef = 0.0_dp
656 mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef = 0.0_dp
657 END DO
658 mixing_store%cpc_h_lastres(iat, ispin)%r_coef = 0.0_dp
659 mixing_store%cpc_s_lastres(iat, ispin)%r_coef = 0.0_dp
660 END IF
661 END DO
662 END DO
663 END IF
664 END IF
665
666 IF (.NOT. ASSOCIATED(mixing_store%p_metric)) THEN
667 ALLOCATE (mixing_store%p_metric(ng))
668 bconst = mixing_store%bconst
669 g2min = 1.e30_dp
670 DO ig = 1, ng
671 IF (g2(ig) > 1.e-10_dp) g2min = min(g2min, g2(ig))
672 END DO
673 g2max = -1.e30_dp
674 DO ig = 1, ng
675 g2max = max(g2max, g2(ig))
676 END DO
677 CALL para_env%min(g2min)
678 CALL para_env%max(g2max)
679 ! fdamp/g2 varies between (bconst-1) and 0
680 ! i.e. p_metric varies between bconst and 1
681 ! fdamp = (bconst-1.0_dp)*g2min
682 fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
683 DO ig = 1, ng
684 mixing_store%p_metric(ig) = (g2(ig) + fdamp)/max(g2(ig), 1.e-10_dp)
685 END DO
686 IF (rho_g(1)%pw_grid%have_g0) mixing_store%p_metric(1) = bconst
687 END IF
688 ELSEIF (mixing_method == multisecant_mixing_nr) THEN
689 IF (.NOT. ASSOCIATED(mixing_store%ig_global_index)) THEN
690 ALLOCATE (mixing_store%ig_global_index(ng))
691 END IF
692 mixing_store%ig_global_index = 0
693 ig_count = 0
694 DO iproc = 0, para_env%num_pe - 1
695 IF (para_env%mepos == iproc) THEN
696 DO ig = 1, ng
697 ig_count = ig_count + 1
698 mixing_store%ig_global_index(ig) = ig_count
699 END DO
700 END IF
701 CALL para_env%bcast(ig_count, iproc)
702 END DO
703 END IF
704
705 CALL timestop(handle)
706
707 END SUBROUTINE mixing_init
708
709! **************************************************************************************************
710!> \brief initialiation needed when charge mixing is used
711!> \param mixing_store ...
712!> \par History
713!> 02.2019 created [JGH]
714!> \author JGH
715! **************************************************************************************************
716 ELEMENTAL SUBROUTINE charge_mixing_init(mixing_store)
717 TYPE(mixing_storage_type), INTENT(INOUT) :: mixing_store
718
719 mixing_store%ncall = 0
720 mixing_store%tb_scc_mixer_error = 0.0_dp
721 mixing_store%tb_scc_mixer_step = 0
722 IF (ASSOCIATED(mixing_store%acharge)) mixing_store%acharge = 0.0_dp
723 IF (ASSOCIATED(mixing_store%dacharge)) mixing_store%dacharge = 0.0_dp
724
725 END SUBROUTINE charge_mixing_init
726
727END 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
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::sundararaman2017
integer, save, public sundararaman2017
Definition bibliography.F:36
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:1203
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::new_pulay_mixing_nr
integer, parameter, public new_pulay_mixing_nr
Definition qs_density_mixing_types.F:50
qs_density_mixing_types::broyden_mixing_nr
integer, parameter, public broyden_mixing_nr
Definition qs_density_mixing_types.F:50
qs_density_mixing_types::modified_broyden_mixing_nr
integer, parameter, public modified_broyden_mixing_nr
Definition qs_density_mixing_types.F:50
qs_density_mixing_types::multisecant_mixing_nr
integer, parameter, public multisecant_mixing_nr
Definition qs_density_mixing_types.F:50
qs_density_mixing_types::pulay_mixing_nr
integer, parameter, public pulay_mixing_nr
Definition qs_density_mixing_types.F:50
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition qs_density_mixing_types.F:50
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:717
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:409
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:116
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:63
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:744
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:62
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