Line data Source code
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 : ! **************************************************************************************************
9 : MODULE 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 : cp_1d_z_p_type,&
27 : gspace_mixing_nr,&
28 : mixing_storage_type,&
29 : modified_broyden_mixing_nr,&
30 : multisecant_mixing_nr,&
31 : new_pulay_mixing_nr,&
32 : pulay_mixing_nr
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
36 : USE qs_rho_atom_types, ONLY: rho_atom_type
37 : USE qs_rho_types, ONLY: qs_rho_get,&
38 : qs_rho_type
39 : USE qs_scf_methods, ONLY: cp_sm_mix
40 : #include "./base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 :
44 : PRIVATE
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mixing_utils'
47 : INTEGER, PARAMETER, PRIVATE :: max_dftb_scc_vars = 1, &
48 : max_xtb_scc_vars = 28
49 :
50 : PUBLIC :: mixing_allocate, mixing_init, mixing_density, charge_mixing_init, &
51 : self_consistency_check
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief ...
57 : !> \param rho_ao ...
58 : !> \param p_delta ...
59 : !> \param para_env ...
60 : !> \param p_out ...
61 : !> \param delta ...
62 : ! **************************************************************************************************
63 11196 : SUBROUTINE self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
64 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao, p_delta
65 : TYPE(mp_para_env_type), POINTER :: para_env
66 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_out
67 : REAL(KIND=dp), INTENT(INOUT) :: delta
68 :
69 : CHARACTER(len=*), PARAMETER :: routineN = 'self_consistency_check'
70 :
71 : INTEGER :: handle, ic, ispin, nimg, nspins
72 : REAL(KIND=dp) :: tmp
73 11196 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_q, p_in
74 :
75 11196 : CALL timeset(routineN, handle)
76 :
77 11196 : NULLIFY (matrix_q, p_in)
78 :
79 11196 : CPASSERT(ASSOCIATED(p_out))
80 11196 : NULLIFY (matrix_q, p_in)
81 11196 : p_in => rho_ao
82 11196 : matrix_q => p_delta
83 11196 : nspins = SIZE(p_in, 1)
84 11196 : nimg = SIZE(p_in, 2)
85 :
86 : ! Compute the difference (p_out - p_in)and check convergence
87 11196 : delta = 0.0_dp
88 23934 : DO ispin = 1, nspins
89 671346 : DO ic = 1, nimg
90 647412 : CALL dbcsr_set(matrix_q(ispin, ic)%matrix, 0.0_dp)
91 : CALL cp_sm_mix(m1=p_out(ispin, ic)%matrix, m2=p_in(ispin, ic)%matrix, &
92 : p_mix=1.0_dp, delta=tmp, para_env=para_env, &
93 647412 : m3=matrix_q(ispin, ic)%matrix)
94 660150 : delta = MAX(tmp, delta)
95 : END DO
96 : END DO
97 11196 : CALL timestop(handle)
98 :
99 11196 : END SUBROUTINE self_consistency_check
100 :
101 : ! **************************************************************************************************
102 : !> \brief allocation needed when density mixing is used
103 : !> \param qs_env ...
104 : !> \param mixing_method ...
105 : !> \param p_mix_new ...
106 : !> \param p_delta ...
107 : !> \param nspins ...
108 : !> \param mixing_store ...
109 : !> \par History
110 : !> 05.2009 created [MI]
111 : !> 08.2014 kpoints [JGH]
112 : !> 02.2015 changed input to make g-space mixing available in linear scaling SCF [Patrick Seewald]
113 : !> 02.2019 charge mixing [JGH]
114 : !> \author fawzi
115 : ! **************************************************************************************************
116 20842 : SUBROUTINE mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
117 : TYPE(qs_environment_type), POINTER :: qs_env
118 : INTEGER :: mixing_method
119 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
120 : POINTER :: p_mix_new, p_delta
121 : INTEGER, INTENT(IN) :: nspins
122 : TYPE(mixing_storage_type), POINTER :: mixing_store
123 :
124 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mixing_allocate'
125 :
126 : INTEGER :: handle, i, ia, iat, ic, ikind, ispin, &
127 : max_shell, na, natom, nbuffer, nel, &
128 : nimg, nkind
129 : LOGICAL :: charge_mixing
130 20842 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
131 : TYPE(dbcsr_type), POINTER :: refmatrix
132 : TYPE(dft_control_type), POINTER :: dft_control
133 : TYPE(distribution_1d_type), POINTER :: distribution_1d
134 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
135 20842 : POINTER :: sab_orb
136 20842 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
137 :
138 20842 : CALL timeset(routineN, handle)
139 :
140 20842 : NULLIFY (matrix_s, dft_control, sab_orb, refmatrix, rho_atom)
141 : CALL get_qs_env(qs_env, &
142 : sab_orb=sab_orb, &
143 : matrix_s_kp=matrix_s, &
144 20842 : dft_control=dft_control)
145 :
146 : charge_mixing = dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb &
147 20842 : .OR. dft_control%qs_control%semi_empirical
148 20842 : refmatrix => matrix_s(1, 1)%matrix
149 20842 : nimg = dft_control%nimages
150 :
151 : ! *** allocate p_mix_new ***
152 20842 : IF (PRESENT(p_mix_new)) THEN
153 20838 : IF (.NOT. ASSOCIATED(p_mix_new)) THEN
154 20802 : CALL dbcsr_allocate_matrix_set(p_mix_new, nspins, nimg)
155 44350 : DO i = 1, nspins
156 328592 : DO ic = 1, nimg
157 284242 : ALLOCATE (p_mix_new(i, ic)%matrix)
158 : CALL dbcsr_create(matrix=p_mix_new(i, ic)%matrix, template=refmatrix, &
159 284242 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
160 284242 : CALL cp_dbcsr_alloc_block_from_nbl(p_mix_new(i, ic)%matrix, sab_orb)
161 307790 : CALL dbcsr_set(p_mix_new(i, ic)%matrix, 0.0_dp)
162 : END DO
163 : END DO
164 : END IF
165 : END IF
166 :
167 : ! *** allocate p_delta ***
168 20842 : IF (PRESENT(p_delta)) THEN
169 20838 : IF (mixing_method >= gspace_mixing_nr) THEN
170 1168 : IF (.NOT. ASSOCIATED(p_delta)) THEN
171 1164 : CALL dbcsr_allocate_matrix_set(p_delta, nspins, nimg)
172 2584 : DO i = 1, nspins
173 90346 : DO ic = 1, nimg
174 87762 : ALLOCATE (p_delta(i, ic)%matrix)
175 : CALL dbcsr_create(matrix=p_delta(i, ic)%matrix, template=refmatrix, &
176 87762 : name="SCF DENSITY", matrix_type=dbcsr_type_symmetric)
177 87762 : CALL cp_dbcsr_alloc_block_from_nbl(p_delta(i, ic)%matrix, sab_orb)
178 89182 : CALL dbcsr_set(p_delta(i, ic)%matrix, 0.0_dp)
179 : END DO
180 : END DO
181 : END IF
182 1168 : CPASSERT(ASSOCIATED(mixing_store))
183 : END IF
184 : END IF
185 :
186 20842 : IF (charge_mixing) THEN
187 : ! *** allocate buffer for TB SCC variable mixing ***
188 13008 : IF (mixing_method >= gspace_mixing_nr) THEN
189 180 : CPASSERT(.NOT. mixing_store%gmix_p)
190 180 : IF (dft_control%qs_control%dftb) THEN
191 : max_shell = max_dftb_scc_vars
192 136 : ELSE IF (dft_control%qs_control%xtb) THEN
193 : max_shell = max_xtb_scc_vars
194 : ELSE
195 0 : CPABORT('UNKNOWN METHOD')
196 : END IF
197 180 : nbuffer = mixing_store%nbuffer
198 180 : mixing_store%ncall = 0
199 180 : CALL get_qs_env(qs_env, local_particles=distribution_1d)
200 180 : nkind = SIZE(distribution_1d%n_el)
201 398 : na = SUM(distribution_1d%n_el(:))
202 180 : IF (ASSOCIATED(mixing_store%atlist)) DEALLOCATE (mixing_store%atlist)
203 540 : ALLOCATE (mixing_store%atlist(na))
204 180 : mixing_store%nat_local = na
205 180 : mixing_store%max_shell = max_shell
206 180 : ia = 0
207 398 : DO ikind = 1, nkind
208 218 : nel = distribution_1d%n_el(ikind)
209 914 : DO iat = 1, nel
210 516 : ia = ia + 1
211 734 : mixing_store%atlist(ia) = distribution_1d%list(ikind)%array(iat)
212 : END DO
213 : END DO
214 180 : IF (ASSOCIATED(mixing_store%acharge)) DEALLOCATE (mixing_store%acharge)
215 900 : ALLOCATE (mixing_store%acharge(na, max_shell, nbuffer))
216 180 : IF (ASSOCIATED(mixing_store%dacharge)) DEALLOCATE (mixing_store%dacharge)
217 720 : ALLOCATE (mixing_store%dacharge(na, max_shell, nbuffer))
218 : END IF
219 13008 : IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
220 0 : IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
221 0 : ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
222 0 : mixing_store%pulay_matrix = 0.0_dp
223 13008 : ELSE IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
224 180 : IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
225 720 : ALLOCATE (mixing_store%abroy(nbuffer, nbuffer))
226 32968 : mixing_store%abroy = 0.0_dp
227 180 : IF (ASSOCIATED(mixing_store%wbroy)) DEALLOCATE (mixing_store%wbroy)
228 540 : ALLOCATE (mixing_store%wbroy(nbuffer))
229 1744 : mixing_store%wbroy = 0.0_dp
230 180 : IF (ASSOCIATED(mixing_store%dfbroy)) DEALLOCATE (mixing_store%dfbroy)
231 900 : ALLOCATE (mixing_store%dfbroy(na, max_shell, nbuffer))
232 137280 : mixing_store%dfbroy = 0.0_dp
233 180 : IF (ASSOCIATED(mixing_store%ubroy)) DEALLOCATE (mixing_store%ubroy)
234 720 : ALLOCATE (mixing_store%ubroy(na, max_shell, nbuffer))
235 137280 : mixing_store%ubroy = 0.0_dp
236 12828 : ELSE IF (mixing_method == multisecant_mixing_nr) THEN
237 0 : CPABORT("multisecant_mixing not available")
238 : END IF
239 : ELSE
240 : ! *** allocate buffer for gspace mixing ***
241 7834 : IF (mixing_method >= gspace_mixing_nr) THEN
242 992 : nbuffer = mixing_store%nbuffer
243 992 : mixing_store%ncall = 0
244 2976 : mixing_store%ncall_p = 0
245 992 : IF (.NOT. ASSOCIATED(mixing_store%rhoin)) THEN
246 2548 : ALLOCATE (mixing_store%rhoin(nspins))
247 1324 : DO ispin = 1, nspins
248 1324 : NULLIFY (mixing_store%rhoin(ispin)%cc)
249 : END DO
250 : END IF
251 :
252 992 : IF (mixing_store%gmix_p .AND. dft_control%qs_control%gapw) THEN
253 20 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
254 20 : natom = SIZE(rho_atom)
255 20 : IF (.NOT. ASSOCIATED(mixing_store%paw)) THEN
256 48 : ALLOCATE (mixing_store%paw(natom))
257 144 : mixing_store%paw = .FALSE.
258 262 : ALLOCATE (mixing_store%cpc_h_in(natom, nspins))
259 246 : ALLOCATE (mixing_store%cpc_s_in(natom, nspins))
260 38 : DO ispin = 1, nspins
261 214 : DO iat = 1, natom
262 176 : NULLIFY (mixing_store%cpc_h_in(iat, ispin)%r_coef)
263 198 : NULLIFY (mixing_store%cpc_s_in(iat, ispin)%r_coef)
264 : END DO
265 : END DO
266 : END IF
267 : END IF
268 : END IF
269 :
270 : ! *** allocare res_buffer if needed
271 7834 : IF (mixing_method >= pulay_mixing_nr) THEN
272 980 : IF (.NOT. ASSOCIATED(mixing_store%res_buffer)) THEN
273 7718 : ALLOCATE (mixing_store%res_buffer(nbuffer, nspins))
274 1300 : DO ispin = 1, nspins
275 5918 : DO i = 1, nbuffer
276 5318 : NULLIFY (mixing_store%res_buffer(i, ispin)%cc)
277 : END DO
278 : END DO
279 : END IF
280 : END IF
281 :
282 : ! *** allocate pulay
283 7834 : IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
284 48 : IF (.NOT. ASSOCIATED(mixing_store%pulay_matrix)) THEN
285 220 : ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
286 : END IF
287 :
288 48 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
289 548 : ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
290 98 : DO ispin = 1, nspins
291 416 : DO i = 1, nbuffer
292 372 : NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
293 : END DO
294 : END DO
295 : END IF
296 48 : IF (mixing_store%gmix_p) THEN
297 2 : IF (dft_control%qs_control%gapw) THEN
298 2 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer)) THEN
299 108 : ALLOCATE (mixing_store%cpc_h_in_buffer(nbuffer, natom, nspins))
300 106 : ALLOCATE (mixing_store%cpc_s_in_buffer(nbuffer, natom, nspins))
301 106 : ALLOCATE (mixing_store%cpc_h_res_buffer(nbuffer, natom, nspins))
302 106 : ALLOCATE (mixing_store%cpc_s_res_buffer(nbuffer, natom, nspins))
303 4 : DO ispin = 1, nspins
304 20 : DO iat = 1, natom
305 98 : DO i = 1, nbuffer
306 80 : NULLIFY (mixing_store%cpc_h_in_buffer(i, iat, ispin)%r_coef)
307 80 : NULLIFY (mixing_store%cpc_s_in_buffer(i, iat, ispin)%r_coef)
308 80 : NULLIFY (mixing_store%cpc_h_res_buffer(i, iat, ispin)%r_coef)
309 96 : NULLIFY (mixing_store%cpc_s_res_buffer(i, iat, ispin)%r_coef)
310 : END DO
311 : END DO
312 : END DO
313 : END IF
314 : END IF
315 : END IF
316 :
317 : END IF
318 : ! *** allocate broyden buffer ***
319 7834 : IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
320 932 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_old)) THEN
321 2314 : ALLOCATE (mixing_store%rhoin_old(nspins))
322 1202 : DO ispin = 1, nspins
323 1202 : NULLIFY (mixing_store%rhoin_old(ispin)%cc)
324 : END DO
325 : END IF
326 932 : IF (.NOT. ASSOCIATED(mixing_store%drho_buffer)) THEN
327 7170 : ALLOCATE (mixing_store%drho_buffer(nbuffer, nspins))
328 2314 : ALLOCATE (mixing_store%last_res(nspins))
329 1202 : DO ispin = 1, nspins
330 4946 : DO i = 1, nbuffer
331 4946 : NULLIFY (mixing_store%drho_buffer(i, ispin)%cc)
332 : END DO
333 1202 : NULLIFY (mixing_store%last_res(ispin)%cc)
334 : END DO
335 : END IF
336 932 : IF (mixing_method == modified_broyden_mixing_nr) THEN
337 2 : IF (.NOT. ASSOCIATED(mixing_store%weight)) THEN
338 8 : ALLOCATE (mixing_store%weight(nbuffer, nspins))
339 : END IF
340 34 : mixing_store%weight = 0.0_dp
341 : END IF
342 932 : IF (mixing_store%gmix_p) THEN
343 :
344 16 : IF (dft_control%qs_control%gapw) THEN
345 16 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_old)) THEN
346 210 : ALLOCATE (mixing_store%cpc_h_old(natom, nspins))
347 198 : ALLOCATE (mixing_store%cpc_s_old(natom, nspins))
348 30 : DO ispin = 1, nspins
349 174 : DO iat = 1, natom
350 144 : NULLIFY (mixing_store%cpc_h_old(iat, ispin)%r_coef)
351 162 : NULLIFY (mixing_store%cpc_s_old(iat, ispin)%r_coef)
352 : END DO
353 : END DO
354 : END IF
355 16 : IF (.NOT. ASSOCIATED(mixing_store%dcpc_h_in)) THEN
356 1326 : ALLOCATE (mixing_store%dcpc_h_in(nbuffer, natom, nspins))
357 1314 : ALLOCATE (mixing_store%dcpc_s_in(nbuffer, natom, nspins))
358 210 : ALLOCATE (mixing_store%cpc_h_lastres(natom, nspins))
359 198 : ALLOCATE (mixing_store%cpc_s_lastres(natom, nspins))
360 30 : DO ispin = 1, nspins
361 174 : DO iat = 1, natom
362 1248 : DO i = 1, nbuffer
363 1104 : NULLIFY (mixing_store%dcpc_h_in(i, iat, ispin)%r_coef)
364 1248 : NULLIFY (mixing_store%dcpc_s_in(i, iat, ispin)%r_coef)
365 : END DO
366 144 : NULLIFY (mixing_store%cpc_h_lastres(iat, ispin)%r_coef)
367 162 : NULLIFY (mixing_store%cpc_s_lastres(iat, ispin)%r_coef)
368 : END DO
369 : END DO
370 : END IF
371 : END IF
372 : END IF
373 : END IF
374 :
375 : ! *** allocate multisecant buffer ***
376 7834 : IF (mixing_method == multisecant_mixing_nr) THEN
377 0 : IF (.NOT. ASSOCIATED(mixing_store%norm_res_buffer)) THEN
378 0 : ALLOCATE (mixing_store%norm_res_buffer(nbuffer, nspins))
379 : END IF
380 : END IF
381 :
382 7834 : IF (mixing_method == multisecant_mixing_nr) THEN
383 0 : IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN
384 0 : ALLOCATE (mixing_store%rhoin_buffer(nbuffer, nspins))
385 0 : DO ispin = 1, nspins
386 0 : DO i = 1, nbuffer
387 0 : NULLIFY (mixing_store%rhoin_buffer(i, ispin)%cc)
388 : END DO
389 : END DO
390 : END IF
391 : END IF
392 :
393 : END IF
394 :
395 20842 : CALL timestop(handle)
396 :
397 20842 : END SUBROUTINE mixing_allocate
398 :
399 : ! **************************************************************************************************
400 : !> \brief Allocate or resize a mixing vector, optionally clearing its contents.
401 : !> \param buffer Vector or array of vectors.
402 : !> \param n Local vector length, including any auxiliary coordinates.
403 : !> \param clear Whether to clear the vector; Pulay retains history across SCF runs.
404 : ! **************************************************************************************************
405 21300 : ELEMENTAL SUBROUTINE mixing_buffer_allocate(buffer, n, clear)
406 : TYPE(cp_1d_z_p_type), INTENT(INOUT) :: buffer
407 : INTEGER, INTENT(IN) :: n
408 : LOGICAL, INTENT(IN) :: clear
409 :
410 21300 : IF (ASSOCIATED(buffer%cc)) THEN
411 10060 : IF (SIZE(buffer%cc) /= n) DEALLOCATE (buffer%cc)
412 : END IF
413 43780 : IF (.NOT. ASSOCIATED(buffer%cc)) ALLOCATE (buffer%cc(n))
414 389691422 : IF (clear) buffer%cc = z_zero
415 21300 : END SUBROUTINE mixing_buffer_allocate
416 :
417 : ! **************************************************************************************************
418 : !> \brief View the physical density, appending auxiliary coordinates on their owning rank.
419 : !> \param mixing_store Owns the optional enlarged work buffer.
420 : !> \param rho_g Physical plane-wave density; its storage and grid are unchanged.
421 : !> \param density Vectors passed to the existing mixing algorithms.
422 : !> \param auxiliary Replicated density coordinates in electrons, scaled to Fourier-density units.
423 : ! **************************************************************************************************
424 9578 : SUBROUTINE mixing_density(mixing_store, rho_g, density, auxiliary)
425 : TYPE(mixing_storage_type), POINTER :: mixing_store
426 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho_g(:)
427 : TYPE(cp_1d_z_p_type), ALLOCATABLE, INTENT(OUT) :: density(:)
428 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: auxiliary(:, :)
429 :
430 : INTEGER :: ispin, ng
431 :
432 9578 : IF (mixing_store%naux > 0) THEN
433 294 : CPASSERT(PRESENT(auxiliary))
434 294 : CPASSERT(SIZE(auxiliary, 1) == mixing_store%naux)
435 294 : CPASSERT(SIZE(auxiliary, 2) == SIZE(rho_g))
436 : END IF
437 9578 : ng = SIZE(rho_g(1)%array)
438 39806 : ALLOCATE (density(SIZE(rho_g)))
439 20650 : DO ispin = 1, SIZE(rho_g)
440 11072 : density(ispin)%cc => rho_g(ispin)%array
441 20650 : IF (ALLOCATED(mixing_store%aux_buffer)) THEN
442 195 : density(ispin)%cc => mixing_store%aux_buffer(:, ispin)
443 1921920 : density(ispin)%cc(:ng) = rho_g(ispin)%array
444 1365 : density(ispin)%cc(ng + 1:) = auxiliary(:, ispin)/rho_g(ispin)%pw_grid%vol
445 : END IF
446 : END DO
447 9578 : END SUBROUTINE mixing_density
448 :
449 : ! **************************************************************************************************
450 : !> \brief initialiation needed when gspace mixing is used
451 : !> \param mixing_method ...
452 : !> \param rho ...
453 : !> \param mixing_store ...
454 : !> \param para_env ...
455 : !> \param rho_atom ...
456 : !> \param auxiliary Replicated real density coordinates, included once in the mixing metric.
457 : !> \par History
458 : !> 05.2009 created [MI]
459 : !> \author MI
460 : ! **************************************************************************************************
461 992 : SUBROUTINE mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom, auxiliary)
462 : INTEGER, INTENT(IN) :: mixing_method
463 : TYPE(qs_rho_type), POINTER :: rho
464 : TYPE(mixing_storage_type), POINTER :: mixing_store
465 : TYPE(mp_para_env_type), POINTER :: para_env
466 : TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
467 : POINTER :: rho_atom
468 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: auxiliary(:, :)
469 :
470 : CHARACTER(len=*), PARAMETER :: routineN = 'mixing_init'
471 :
472 : INTEGER :: handle, iat, ib, icoef1, icoef2, ig, ig1, ig_count, iproc, ispin, n1, n2, natom, &
473 : nbuffer, ncoef_h1, ncoef_h2, ncoef_s1, ncoef_s2, ng, nimg, nspin, nvec
474 : REAL(dp) :: bconst, beta, cpc_h_tmp, cpc_s_tmp, &
475 : fdamp, g2max, g2min, kmin, qk, qkappa, &
476 : qm
477 992 : REAL(dp), DIMENSION(:), POINTER :: g2
478 992 : REAL(dp), DIMENSION(:, :), POINTER :: g_vec
479 992 : TYPE(cp_1d_z_p_type), ALLOCATABLE :: density(:)
480 992 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
481 992 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
482 :
483 992 : CALL timeset(routineN, handle)
484 :
485 992 : NULLIFY (g2, g_vec, rho_ao_kp, rho_g)
486 992 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_g=rho_g)
487 :
488 992 : nspin = SIZE(rho_g)
489 992 : ng = SIZE(rho_g(1)%pw_grid%gsq, 1)
490 992 : nimg = SIZE(rho_ao_kp, 2)
491 992 : IF (PRESENT(auxiliary) .OR. mixing_store%naux > 0) THEN
492 76 : mixing_store%ncall = 0
493 1144 : mixing_store%ncall_p = 0
494 : END IF
495 992 : mixing_store%naux = 0
496 992 : IF (PRESENT(auxiliary)) mixing_store%naux = SIZE(auxiliary, 1)
497 992 : nvec = ng
498 992 : IF (ALLOCATED(mixing_store%aux_buffer)) DEALLOCATE (mixing_store%aux_buffer)
499 992 : IF (mixing_store%naux > 0) THEN
500 76 : IF (para_env%mepos == 0) THEN
501 38 : nvec = ng + mixing_store%naux
502 152 : ALLOCATE (mixing_store%aux_buffer(nvec, nspin))
503 : END IF
504 : END IF
505 1908 : CALL mixing_density(mixing_store, rho_g, density, auxiliary)
506 992 : mixing_store%ig_max = ng
507 992 : g2 => rho_g(1)%pw_grid%gsq
508 992 : g_vec => rho_g(1)%pw_grid%g
509 :
510 992 : IF (mixing_store%max_gvec_exp > 0._dp) THEN
511 0 : DO ig = 1, ng
512 0 : IF (g2(ig) > mixing_store%max_g2) THEN
513 0 : mixing_store%ig_max = ig
514 0 : EXIT
515 : END IF
516 : END DO
517 : END IF
518 :
519 992 : IF (ASSOCIATED(mixing_store%kerker_factor)) DEALLOCATE (mixing_store%kerker_factor)
520 2976 : ALLOCATE (mixing_store%kerker_factor(nvec))
521 992 : IF (ASSOCIATED(mixing_store%kerker_factor_mag)) DEALLOCATE (mixing_store%kerker_factor_mag)
522 1984 : ALLOCATE (mixing_store%kerker_factor_mag(nvec))
523 992 : IF (ASSOCIATED(mixing_store%special_metric)) DEALLOCATE (mixing_store%special_metric)
524 1984 : ALLOCATE (mixing_store%special_metric(nvec))
525 992 : beta = mixing_store%beta
526 992 : kmin = 0.1_dp
527 43431102 : mixing_store%kerker_factor = 1.0_dp
528 43431102 : mixing_store%kerker_factor_mag = 1.0_dp
529 43431102 : mixing_store%special_metric = 1.0_dp
530 992 : ig1 = 1
531 :
532 992 : IF (.NOT. (mixing_method == new_pulay_mixing_nr)) THEN
533 982 : IF (rho_g(1)%pw_grid%have_g0) ig1 = 2
534 42608613 : DO ig = ig1, mixing_store%ig_max
535 42607631 : mixing_store%kerker_factor(ig) = MAX(g2(ig)/(g2(ig) + beta*beta), kmin)
536 42607631 : IF (mixing_store%beta_mag > 1.0E-10_dp) THEN
537 42355400 : mixing_store%kerker_factor_mag(ig) = MAX(g2(ig)/(g2(ig) + mixing_store%beta_mag*mixing_store%beta_mag), kmin)
538 : ELSE
539 : ! beta_mag ~ 0 means no Kerker screening on the magnetization channel
540 252231 : mixing_store%kerker_factor_mag(ig) = 1.0_dp
541 : END IF
542 : mixing_store%special_metric(ig) = &
543 : 1.0_dp + 50.0_dp/8.0_dp*( &
544 : 1.0_dp + COS(g_vec(1, ig)) + COS(g_vec(2, ig)) + COS(g_vec(3, ig)) + &
545 : COS(g_vec(1, ig))*COS(g_vec(2, ig)) + &
546 : COS(g_vec(2, ig))*COS(g_vec(3, ig)) + &
547 : COS(g_vec(1, ig))*COS(g_vec(3, ig)) + &
548 42608613 : COS(g_vec(1, ig))*COS(g_vec(2, ig))*COS(g_vec(3, ig)))
549 : END DO
550 : ELSE
551 10 : CALL cite_reference(Sundararaman2017)
552 10 : qkappa = mixing_store%qkappa
553 10 : qk = mixing_store%qk
554 10 : qm = mixing_store%qm
555 : IF (rho_g(1)%pw_grid%have_g0) ig1 = 1
556 821770 : DO ig = ig1, mixing_store%ig_max
557 821760 : mixing_store%kerker_factor(ig) = (g2(ig) + qkappa*qkappa)/(g2(ig) + qkappa*qkappa + qk*qk)
558 821770 : mixing_store%special_metric(ig) = (g2(ig) + qkappa*qkappa + qm*qm)/(g2(ig) + qkappa*qkappa)
559 : END DO
560 : END IF
561 992 : nbuffer = mixing_store%nbuffer
562 2240 : CALL mixing_buffer_allocate(mixing_store%rhoin, nvec, .FALSE.)
563 1400 : IF (ASSOCIATED(mixing_store%rhoin_buffer)) CALL mixing_buffer_allocate(mixing_store%rhoin_buffer, nvec, .FALSE.)
564 11076 : IF (ASSOCIATED(mixing_store%res_buffer)) CALL mixing_buffer_allocate(mixing_store%res_buffer, nvec, .FALSE.)
565 2240 : DO ispin = 1, nspin
566 50885561 : mixing_store%rhoin(ispin)%cc = density(ispin)%cc
567 : END DO
568 992 : IF (nspin == 2) THEN
569 7453467 : mixing_store%rhoin(1)%cc = density(1)%cc + density(2)%cc
570 7453467 : mixing_store%rhoin(2)%cc = density(1)%cc - density(2)%cc
571 : END IF
572 992 : IF (ASSOCIATED(mixing_store%rhoin_buffer)) THEN
573 106 : DO ispin = 1, nspin
574 2670442 : mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
575 : END DO
576 : END IF
577 :
578 992 : IF (mixing_store%gmix_p) THEN
579 20 : IF (PRESENT(rho_atom)) THEN
580 20 : natom = SIZE(rho_atom)
581 50 : DO ispin = 1, nspin
582 290 : DO iat = 1, natom
583 270 : IF (ASSOCIATED(rho_atom(iat)%cpc_s(ispin)%r_coef)) THEN
584 170 : mixing_store%paw(iat) = .TRUE.
585 170 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
586 170 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
587 170 : IF (ASSOCIATED(mixing_store%cpc_s_in)) THEN
588 170 : IF (.NOT. ASSOCIATED(mixing_store%cpc_s_in(iat, ispin)%r_coef)) THEN
589 424 : ALLOCATE (mixing_store%cpc_s_in(iat, ispin)%r_coef(n1, n2))
590 318 : ALLOCATE (mixing_store%cpc_h_in(iat, ispin)%r_coef(n1, n2))
591 : END IF
592 251330 : mixing_store%cpc_h_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_h(ispin)%r_coef
593 251330 : mixing_store%cpc_s_in(iat, ispin)%r_coef = rho_atom(iat)%cpc_s(ispin)%r_coef
594 : END IF
595 : END IF
596 : END DO
597 : END DO
598 20 : IF (nspin == 2 .AND. ASSOCIATED(mixing_store%cpc_s_in)) THEN
599 90 : DO iat = 1, natom
600 90 : IF (mixing_store%paw(iat)) THEN
601 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_h_in(iat, 1)%r_coef))
602 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_h_in(iat, 2)%r_coef))
603 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_s_in(iat, 1)%r_coef))
604 80 : CPASSERT(ASSOCIATED(mixing_store%cpc_s_in(iat, 2)%r_coef))
605 :
606 80 : ncoef_h1 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 1)
607 80 : ncoef_h2 = SIZE(mixing_store%cpc_h_in(iat, 1)%r_coef, 2)
608 80 : ncoef_s1 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 1)
609 80 : ncoef_s2 = SIZE(mixing_store%cpc_s_in(iat, 1)%r_coef, 2)
610 80 : CPASSERT(ncoef_h1 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 1))
611 80 : CPASSERT(ncoef_h2 == SIZE(mixing_store%cpc_h_in(iat, 2)%r_coef, 2))
612 80 : CPASSERT(ncoef_s1 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 1))
613 80 : CPASSERT(ncoef_s2 == SIZE(mixing_store%cpc_s_in(iat, 2)%r_coef, 2))
614 :
615 2240 : DO icoef2 = 1, ncoef_h2
616 60560 : DO icoef1 = 1, ncoef_h1
617 58320 : cpc_h_tmp = mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2)
618 : mixing_store%cpc_h_in(iat, 1)%r_coef(icoef1, icoef2) = &
619 58320 : cpc_h_tmp + mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
620 : mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2) = &
621 60480 : cpc_h_tmp - mixing_store%cpc_h_in(iat, 2)%r_coef(icoef1, icoef2)
622 : END DO
623 : END DO
624 :
625 2240 : DO icoef2 = 1, ncoef_s2
626 60560 : DO icoef1 = 1, ncoef_s1
627 58320 : cpc_s_tmp = mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2)
628 : mixing_store%cpc_s_in(iat, 1)%r_coef(icoef1, icoef2) = &
629 58320 : cpc_s_tmp + mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
630 : mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2) = &
631 60480 : cpc_s_tmp - mixing_store%cpc_s_in(iat, 2)%r_coef(icoef1, icoef2)
632 : END DO
633 : END DO
634 : END IF
635 : END DO
636 : END IF
637 : END IF
638 : END IF
639 :
640 : IF (mixing_method == gspace_mixing_nr) THEN
641 : ELSE IF (mixing_method == pulay_mixing_nr .OR. mixing_method == new_pulay_mixing_nr) THEN
642 48 : IF (mixing_store%gmix_p .AND. PRESENT(rho_atom)) THEN
643 4 : DO ispin = 1, nspin
644 20 : DO iat = 1, natom
645 18 : IF (mixing_store%paw(iat)) THEN
646 2 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
647 2 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
648 2 : IF (.NOT. ASSOCIATED(mixing_store%cpc_h_in_buffer(1, iat, ispin)%r_coef)) THEN
649 12 : DO ib = 1, nbuffer
650 40 : ALLOCATE (mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
651 30 : ALLOCATE (mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef(n1, n2))
652 30 : ALLOCATE (mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
653 32 : ALLOCATE (mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef(n1, n2))
654 : END DO
655 : END IF
656 12 : DO ib = 1, nbuffer
657 4630 : mixing_store%cpc_h_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
658 4630 : mixing_store%cpc_s_in_buffer(ib, iat, ispin)%r_coef = 0.0_dp
659 4630 : mixing_store%cpc_h_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
660 4632 : mixing_store%cpc_s_res_buffer(ib, iat, ispin)%r_coef = 0.0_dp
661 : END DO
662 : END IF
663 : END DO
664 : END DO
665 : END IF
666 : ELSE IF (mixing_method == broyden_mixing_nr .OR. mixing_method == modified_broyden_mixing_nr) THEN
667 2110 : CALL mixing_buffer_allocate(mixing_store%rhoin_old, nvec, .TRUE.)
668 10608 : CALL mixing_buffer_allocate(mixing_store%drho_buffer, nvec, .TRUE.)
669 2110 : CALL mixing_buffer_allocate(mixing_store%last_res, nvec, .TRUE.)
670 932 : IF (mixing_store%gmix_p) THEN
671 16 : IF (PRESENT(rho_atom)) THEN
672 42 : DO ispin = 1, nspin
673 250 : DO iat = 1, natom
674 234 : IF (mixing_store%paw(iat)) THEN
675 166 : n1 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 1)
676 166 : n2 = SIZE(rho_atom(iat)%cpc_s(ispin)%r_coef, 2)
677 166 : IF (.NOT. ASSOCIATED(mixing_store%cpc_s_old(iat, ispin)%r_coef)) THEN
678 408 : ALLOCATE (mixing_store%cpc_s_old(iat, ispin)%r_coef(n1, n2))
679 306 : ALLOCATE (mixing_store%cpc_h_old(iat, ispin)%r_coef(n1, n2))
680 : END IF
681 123898 : mixing_store%cpc_h_old(iat, ispin)%r_coef = 0.0_dp
682 123898 : mixing_store%cpc_s_old(iat, ispin)%r_coef = 0.0_dp
683 166 : IF (.NOT. ASSOCIATED(mixing_store%dcpc_s_in(1, iat, ispin)%r_coef)) THEN
684 912 : DO ib = 1, nbuffer
685 3240 : ALLOCATE (mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef(n1, n2))
686 2532 : ALLOCATE (mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef(n1, n2))
687 : END DO
688 408 : ALLOCATE (mixing_store%cpc_h_lastres(iat, ispin)%r_coef(n1, n2))
689 306 : ALLOCATE (mixing_store%cpc_s_lastres(iat, ispin)%r_coef(n1, n2))
690 : END IF
691 1488 : DO ib = 1, nbuffer
692 988406 : mixing_store%dcpc_h_in(ib, iat, ispin)%r_coef = 0.0_dp
693 988572 : mixing_store%dcpc_s_in(ib, iat, ispin)%r_coef = 0.0_dp
694 : END DO
695 123898 : mixing_store%cpc_h_lastres(iat, ispin)%r_coef = 0.0_dp
696 123898 : mixing_store%cpc_s_lastres(iat, ispin)%r_coef = 0.0_dp
697 : END IF
698 : END DO
699 : END DO
700 : END IF
701 : END IF
702 :
703 932 : IF (ASSOCIATED(mixing_store%p_metric)) DEALLOCATE (mixing_store%p_metric)
704 1864 : ALLOCATE (mixing_store%p_metric(nvec))
705 41104130 : mixing_store%p_metric = 1.0_dp
706 932 : bconst = mixing_store%bconst
707 932 : g2min = 1.E30_dp
708 41103902 : DO ig = 1, ng
709 41103902 : IF (g2(ig) > 1.E-10_dp) g2min = MIN(g2min, g2(ig))
710 : END DO
711 932 : g2max = -1.E30_dp
712 41103902 : DO ig = 1, ng
713 41103902 : g2max = MAX(g2max, g2(ig))
714 : END DO
715 932 : CALL para_env%min(g2min)
716 932 : CALL para_env%max(g2max)
717 : ! fdamp/g2 varies between (bconst-1) and 0
718 : ! i.e. p_metric varies between bconst and 1
719 : ! fdamp = (bconst-1.0_dp)*g2min
720 932 : fdamp = (bconst - 1.0_dp)*g2min*g2max/(g2max - g2min*bconst)
721 41103902 : DO ig = 1, ng
722 41103902 : mixing_store%p_metric(ig) = (g2(ig) + fdamp)/MAX(g2(ig), 1.E-10_dp)
723 : END DO
724 932 : IF (rho_g(1)%pw_grid%have_g0) mixing_store%p_metric(1) = bconst
725 : ! Keep the original G cutoff; auxiliary coordinates have unit metric.
726 932 : IF (nvec > ng) mixing_store%p_metric(mixing_store%ig_max + 1:ng) = 0.0_dp
727 : ELSE IF (mixing_method == multisecant_mixing_nr) THEN
728 0 : IF (ASSOCIATED(mixing_store%ig_global_index)) DEALLOCATE (mixing_store%ig_global_index)
729 0 : ALLOCATE (mixing_store%ig_global_index(nvec))
730 0 : mixing_store%ig_global_index = 0
731 0 : ig_count = 0
732 0 : DO iproc = 0, para_env%num_pe - 1
733 0 : IF (para_env%mepos == iproc) THEN
734 0 : DO ig = 1, nvec
735 0 : ig_count = ig_count + 1
736 0 : mixing_store%ig_global_index(ig) = ig_count
737 : END DO
738 : END IF
739 0 : CALL para_env%bcast(ig_count, iproc)
740 : END DO
741 : END IF
742 :
743 992 : IF (nvec > ng) mixing_store%ig_max = nvec
744 992 : CALL timestop(handle)
745 :
746 1984 : END SUBROUTINE mixing_init
747 :
748 : ! **************************************************************************************************
749 : !> \brief initialiation needed when charge mixing is used
750 : !> \param mixing_store ...
751 : !> \par History
752 : !> 02.2019 created [JGH]
753 : !> \author JGH
754 : ! **************************************************************************************************
755 180 : ELEMENTAL SUBROUTINE charge_mixing_init(mixing_store)
756 : TYPE(mixing_storage_type), INTENT(INOUT) :: mixing_store
757 :
758 180 : mixing_store%ncall = 0
759 180 : mixing_store%tb_scc_mixer_error = 0.0_dp
760 180 : mixing_store%tb_scc_mixer_step = 0
761 137280 : IF (ASSOCIATED(mixing_store%acharge)) mixing_store%acharge = 0.0_dp
762 137280 : IF (ASSOCIATED(mixing_store%dacharge)) mixing_store%dacharge = 0.0_dp
763 :
764 180 : END SUBROUTINE charge_mixing_init
765 :
766 : END MODULE qs_mixing_utils
|
