13 USE mathlib,
ONLY: invert_matrix
19 pw_integrate_function,&
38 #include "./base/base_uses.f90"
44 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_gspace_mixing'
65 SUBROUTINE gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
66 TYPE(qs_environment_type),
POINTER :: qs_env
67 INTEGER :: mixing_method
68 TYPE(mixing_storage_type),
POINTER :: mixing_store
69 TYPE(qs_rho_type),
POINTER :: rho
70 TYPE(mp_para_env_type),
POINTER :: para_env
73 CHARACTER(len=*),
PARAMETER :: routinen =
'gspace_mixing'
75 INTEGER :: handle, iatom, ig, ispin, natom, ng, &
79 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tot_rho_r
80 TYPE(dft_control_type),
POINTER :: dft_control
81 TYPE(pw_c1d_gs_type) :: rho_tmp
82 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g
83 TYPE(pw_env_type),
POINTER :: pw_env
84 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
85 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rho_r
86 TYPE(rho_atom_type),
DIMENSION(:),
POINTER :: rho_atom
88 CALL timeset(routinen, handle)
90 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env)
91 IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb .OR. &
92 dft_control%qs_control%semi_empirical))
THEN
94 cpassert(
ASSOCIATED(mixing_store))
95 NULLIFY (auxbas_pw_pool, rho_atom, rho_g, rho_r, tot_rho_r)
96 CALL qs_rho_get(rho, rho_g=rho_g, rho_r=rho_r, tot_rho_r=tot_rho_r)
98 nspin =
SIZE(rho_g, 1)
99 nimg = dft_control%nimages
100 ng =
SIZE(rho_g(1)%pw_grid%gsq)
101 cpassert(ng ==
SIZE(mixing_store%rhoin(1)%cc))
103 alpha = mixing_store%alpha
104 gapw = dft_control%qs_control%gapw
107 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
108 CALL auxbas_pw_pool%create_pw(rho_tmp)
109 CALL pw_zero(rho_tmp)
110 CALL pw_copy(rho_g(1), rho_tmp)
111 CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
112 CALL pw_axpy(rho_tmp, rho_g(2), -1.0_dp)
113 CALL pw_scale(rho_g(2), -1.0_dp)
116 IF (iter_count + 1 <= mixing_store%nskip_mixing)
THEN
120 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
123 IF (mixing_store%gmix_p .AND. gapw)
THEN
124 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
125 natom =
SIZE(rho_atom)
128 IF (mixing_store%paw(iatom))
THEN
129 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
130 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
136 mixing_store%iter_method =
"NoMix"
138 CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
139 CALL pw_scale(rho_g(1), 0.5_dp)
140 CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
141 CALL pw_scale(rho_g(2), -1.0_dp)
142 CALL auxbas_pw_pool%give_back_pw(rho_tmp)
144 CALL timestop(handle)
148 IF ((iter_count + 1 - mixing_store%nskip_mixing) <= mixing_store%n_simple_mix)
THEN
149 CALL gmix_potential_only(qs_env, mixing_store, rho)
150 mixing_store%iter_method =
"Kerker"
152 CALL gmix_potential_only(qs_env, mixing_store, rho)
153 mixing_store%iter_method =
"Kerker"
155 CALL pulay_mixing(qs_env, mixing_store, rho, para_env)
156 mixing_store%iter_method =
"Pulay"
158 CALL broyden_mixing(qs_env, mixing_store, rho, para_env)
159 mixing_store%iter_method =
"Broy."
162 CALL broyden_mixing_new(mixing_store, rho, para_env)
163 mixing_store%iter_method =
"Broy."
166 CALL multisecant_mixing(mixing_store, rho, para_env)
167 mixing_store%iter_method =
"MSec."
171 CALL pw_axpy(rho_g(2), rho_g(1), 1.0_dp)
172 CALL pw_scale(rho_g(1), 0.5_dp)
173 CALL pw_axpy(rho_g(1), rho_g(2), -1.0_dp)
174 CALL pw_scale(rho_g(2), -1.0_dp)
175 CALL auxbas_pw_pool%give_back_pw(rho_tmp)
179 CALL pw_transfer(rho_g(ispin), rho_r(ispin))
180 tot_rho_r(ispin) = pw_integrate_function(rho_r(ispin), isign=-1)
184 CALL timestop(handle)
199 SUBROUTINE gmix_potential_only(qs_env, mixing_store, rho)
201 TYPE(qs_environment_type),
POINTER :: qs_env
202 TYPE(mixing_storage_type),
POINTER :: mixing_store
203 TYPE(qs_rho_type),
POINTER :: rho
205 CHARACTER(len=*),
PARAMETER :: routinen =
'gmix_potential_only'
207 COMPLEX(dp),
DIMENSION(:),
POINTER :: cc_new
208 INTEGER :: handle, iatom, ig, ispin, natom, ng, &
211 REAL(
dp) :: alpha, f_mix
212 TYPE(dft_control_type),
POINTER :: dft_control
213 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g
214 TYPE(rho_atom_type),
DIMENSION(:),
POINTER :: rho_atom
216 cpassert(
ASSOCIATED(mixing_store%rhoin))
217 cpassert(
ASSOCIATED(mixing_store%kerker_factor))
219 CALL timeset(routinen, handle)
220 NULLIFY (cc_new, dft_control, rho_atom, rho_g)
222 CALL get_qs_env(qs_env, dft_control=dft_control)
225 nspin =
SIZE(rho_g, 1)
226 ng =
SIZE(rho_g(1)%pw_grid%gsq)
228 gapw = dft_control%qs_control%gapw
229 alpha = mixing_store%alpha
232 cc_new => rho_g(ispin)%array
233 DO ig = 1, mixing_store%ig_max
234 f_mix = mixing_store%alpha*mixing_store%kerker_factor(ig)
235 cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
236 mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
238 DO ig = mixing_store%ig_max + 1, ng
239 f_mix = mixing_store%alpha
240 cc_new(ig) = (1.0_dp - f_mix)*mixing_store%rhoin(ispin)%cc(ig) + f_mix*cc_new(ig)
241 mixing_store%rhoin(ispin)%cc(ig) = cc_new(ig)
246 IF (mixing_store%gmix_p .AND. gapw)
THEN
248 rho_atom_set=rho_atom)
249 natom =
SIZE(rho_atom)
252 IF (mixing_store%paw(iatom))
THEN
253 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
254 mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha)
255 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
256 mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha)
257 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
258 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
264 CALL timestop(handle)
266 END SUBROUTINE gmix_potential_only
282 SUBROUTINE pulay_mixing(qs_env, mixing_store, rho, para_env)
284 TYPE(qs_environment_type),
POINTER :: qs_env
285 TYPE(mixing_storage_type),
POINTER :: mixing_store
286 TYPE(qs_rho_type),
POINTER :: rho
287 TYPE(mp_para_env_type),
POINTER :: para_env
289 CHARACTER(len=*),
PARAMETER :: routinen =
'pulay_mixing'
291 COMPLEX(dp),
ALLOCATABLE,
DIMENSION(:) :: cc_mix
292 INTEGER :: handle, i, iatom, ib, ibb, ig, ispin, &
293 jb, n1, n2, natom, nb, nb1, nbuffer, &
296 REAL(
dp) :: alpha_kerker, alpha_pulay, beta, f_mix, &
297 inv_err, norm, norm_c_inv, res_norm, &
299 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: alpha_c, ev
300 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: a, b, c, c_inv, cpc_h_mix, cpc_s_mix
301 TYPE(dft_control_type),
POINTER :: dft_control
302 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g
303 TYPE(rho_atom_type),
DIMENSION(:),
POINTER :: rho_atom
305 cpassert(
ASSOCIATED(mixing_store%res_buffer))
306 cpassert(
ASSOCIATED(mixing_store%rhoin_buffer))
307 NULLIFY (dft_control, rho_atom, rho_g)
308 CALL timeset(routinen, handle)
310 CALL get_qs_env(qs_env, dft_control=dft_control)
312 nspin =
SIZE(rho_g, 1)
313 ng =
SIZE(mixing_store%res_buffer(1, 1)%cc)
314 vol = rho_g(1)%pw_grid%vol
316 alpha_kerker = mixing_store%alpha
317 beta = mixing_store%pulay_beta
318 alpha_pulay = mixing_store%pulay_alpha
319 nbuffer = mixing_store%nbuffer
320 gapw = dft_control%qs_control%gapw
322 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
323 natom =
SIZE(rho_atom)
326 ALLOCATE (cc_mix(ng))
329 CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
331 mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
333 IF (gapw .AND. mixing_store%gmix_p)
THEN
335 IF (mixing_store%paw(iatom))
THEN
336 rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
337 rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
338 mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
339 mixing_store%cpc_h_in(iatom, 2)%r_coef
340 mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
341 mixing_store%cpc_s_in(iatom, 2)%r_coef
348 ib =
modulo(mixing_store%ncall_p(ispin), nbuffer) + 1
350 mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) + 1
351 nb = min(mixing_store%ncall_p(ispin), nbuffer)
352 ibb =
modulo(mixing_store%ncall_p(ispin), nbuffer) + 1
354 mixing_store%res_buffer(ib, ispin)%cc(:) = cmplx(0._dp, 0._dp, kind=
dp)
356 IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
359 f_mix = mixing_store%kerker_factor(ig)
360 mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(rho_g(ispin)%array(ig) - &
361 mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
362 res_norm = res_norm + &
363 REAL(mixing_store%res_buffer(ib, ispin)%cc(ig),
dp)*
REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
364 aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))*aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))
367 CALL para_env%sum(res_norm)
368 res_norm = sqrt(res_norm)
370 IF (res_norm < 1.e-14_dp)
THEN
371 mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) - 1
374 ALLOCATE (a(nb1, nb1))
376 ALLOCATE (b(nb1, nb1))
380 ALLOCATE (c_inv(nb, nb))
382 ALLOCATE (alpha_c(nb))
391 mixing_store%pulay_matrix(jb, ib, ispin) = 0.0_dp
393 f_mix = mixing_store%special_metric(ig)
394 mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin) + &
395 f_mix*(real(mixing_store%res_buffer(jb, ispin)%cc(ig),
dp) &
396 *real(mixing_store%res_buffer(ib, ispin)%cc(ig),
dp) + &
397 aimag(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
398 aimag(mixing_store%res_buffer(ib, ispin)%cc(ig)))
400 CALL para_env%sum(mixing_store%pulay_matrix(jb, ib, ispin))
401 mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)*vol*vol
402 mixing_store%pulay_matrix(ib, jb, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)
405 IF (nb == 1 .OR. res_norm < 1.e-14_dp)
THEN
407 f_mix = alpha_kerker*mixing_store%kerker_factor(ig)
408 cc_mix(ig) = rho_g(ispin)%array(ig) - &
409 mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
410 rho_g(ispin)%array(ig) = f_mix*cc_mix(ig) + &
411 mixing_store%rhoin_buffer(ib, ispin)%cc(ig)
412 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
414 IF (mixing_store%gmix_p)
THEN
417 IF (mixing_store%paw(iatom))
THEN
418 mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
419 mixing_store%cpc_h_in(iatom, ispin)%r_coef
420 mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
421 mixing_store%cpc_s_in(iatom, ispin)%r_coef
423 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
424 (1.0_dp - alpha_kerker)*mixing_store%cpc_h_in(iatom, ispin)%r_coef
425 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_kerker*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
426 (1.0_dp - alpha_kerker)*mixing_store%cpc_s_in(iatom, ispin)%r_coef
428 mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
429 mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
431 mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
432 mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
439 b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb, ispin)
440 c(1:nb, 1:nb) = b(1:nb, 1:nb)
441 b(nb1, 1:nb) = -1.0_dp
442 b(1:nb, nb1) = -1.0_dp
445 CALL invert_matrix(c, c_inv, inv_err, improve=.true.)
449 alpha_c(i) = alpha_c(i) + c_inv(jb, i)
450 norm_c_inv = norm_c_inv + c_inv(jb, i)
453 alpha_c(1:nb) = alpha_c(1:nb)/norm_c_inv
454 cc_mix = cmplx(0._dp, 0._dp, kind=
dp)
457 cc_mix(ig) = cc_mix(ig) + &
458 alpha_c(jb)*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) + &
459 mixing_store%pulay_beta*mixing_store%res_buffer(jb, ispin)%cc(ig))
462 mixing_store%rhoin_buffer(ibb, ispin)%cc = cmplx(0._dp, 0._dp, kind=
dp)
463 IF (alpha_pulay > 0.0_dp)
THEN
465 f_mix = alpha_pulay*mixing_store%kerker_factor(ig)
466 rho_g(ispin)%array(ig) = f_mix*rho_g(ispin)%array(ig) + &
467 (1.0_dp - f_mix)*cc_mix(ig)
468 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
472 rho_g(ispin)%array(ig) = cc_mix(ig)
473 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
477 IF (mixing_store%gmix_p .AND. gapw)
THEN
478 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
480 IF (mixing_store%paw(iatom))
THEN
481 n1 =
SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 1)
482 n2 =
SIZE(rho_atom(iatom)%cpc_h(ispin)%r_coef, 2)
483 ALLOCATE (cpc_h_mix(n1, n2))
484 ALLOCATE (cpc_s_mix(n1, n2))
486 mixing_store%cpc_h_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
487 mixing_store%cpc_h_in_buffer(ib, iatom, ispin)%r_coef
488 mixing_store%cpc_s_res_buffer(ib, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
489 mixing_store%cpc_s_in_buffer(ib, iatom, ispin)%r_coef
493 cpc_h_mix(:, :) = cpc_h_mix(:, :) + &
494 alpha_c(jb)*mixing_store%cpc_h_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
495 alpha_c(jb)*beta*mixing_store%cpc_h_res_buffer(jb, iatom, ispin)%r_coef(:, :)
496 cpc_s_mix(:, :) = cpc_s_mix(:, :) + &
497 alpha_c(jb)*mixing_store%cpc_s_in_buffer(jb, iatom, ispin)%r_coef(:, :) + &
498 alpha_c(jb)*beta*mixing_store%cpc_s_res_buffer(jb, iatom, ispin)%r_coef(:, :)
500 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
501 (1.0_dp - alpha_pulay)*cpc_h_mix
502 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha_pulay*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
503 (1.0_dp - alpha_pulay)*cpc_s_mix
504 mixing_store%cpc_h_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
505 mixing_store%cpc_s_in_buffer(ibb, iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
506 DEALLOCATE (cpc_h_mix)
507 DEALLOCATE (cpc_s_mix)
513 IF (res_norm > 1.e-14_dp)
THEN
517 DEALLOCATE (c, c_inv, alpha_c)
525 CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
527 mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
529 IF (gapw .AND. mixing_store%gmix_p)
THEN
531 IF (mixing_store%paw(iatom))
THEN
532 rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
533 rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
534 mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
535 mixing_store%cpc_h_in(iatom, 2)%r_coef
536 mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
537 mixing_store%cpc_s_in(iatom, 2)%r_coef
544 CALL timestop(handle)
546 END SUBROUTINE pulay_mixing
560 SUBROUTINE broyden_mixing(qs_env, mixing_store, rho, para_env)
562 TYPE(qs_environment_type),
POINTER :: qs_env
563 TYPE(mixing_storage_type),
POINTER :: mixing_store
564 TYPE(qs_rho_type),
POINTER :: rho
565 TYPE(mp_para_env_type),
POINTER :: para_env
567 CHARACTER(len=*),
PARAMETER :: routinen =
'broyden_mixing'
569 COMPLEX(dp) :: cc_mix
570 COMPLEX(dp),
ALLOCATABLE,
DIMENSION(:) :: res_rho
571 INTEGER :: handle, i, iatom, ib, ig, ispin, j, jb, &
572 kb, n1, n2, natom, nb, nbuffer, ng, &
574 LOGICAL :: gapw, only_kerker
575 REAL(
dp) :: alpha, dcpc_h_res, dcpc_s_res, &
576 delta_norm, f_mix, inv_err, res_norm, &
577 rho_norm, valh, vals, w0
578 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: c, g
579 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: a, b
580 TYPE(dft_control_type),
POINTER :: dft_control
581 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g
582 TYPE(rho_atom_type),
DIMENSION(:),
POINTER :: rho_atom
584 cpassert(
ASSOCIATED(mixing_store%res_buffer))
585 cpassert(
ASSOCIATED(mixing_store%rhoin))
586 cpassert(
ASSOCIATED(mixing_store%rhoin_old))
587 cpassert(
ASSOCIATED(mixing_store%drho_buffer))
588 NULLIFY (dft_control, rho_atom, rho_g)
589 CALL timeset(routinen, handle)
591 CALL get_qs_env(qs_env, dft_control=dft_control)
594 nspin =
SIZE(rho_g, 1)
595 ng =
SIZE(mixing_store%res_buffer(1, 1)%cc)
597 alpha = mixing_store%alpha
598 w0 = mixing_store%broy_w0
599 nbuffer = mixing_store%nbuffer
600 gapw = dft_control%qs_control%gapw
602 ALLOCATE (res_rho(ng))
604 mixing_store%ncall = mixing_store%ncall + 1
605 IF (mixing_store%ncall == 1)
THEN
608 only_kerker = .false.
609 nb = min(mixing_store%ncall - 1, nbuffer)
610 ib =
modulo(mixing_store%ncall - 2, nbuffer) + 1
614 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
615 natom =
SIZE(rho_atom)
621 CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
623 mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
625 IF (gapw .AND. mixing_store%gmix_p)
THEN
627 IF (mixing_store%paw(iatom))
THEN
628 rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
629 rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
630 mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
631 mixing_store%cpc_h_in(iatom, 2)%r_coef
632 mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
633 mixing_store%cpc_s_in(iatom, 2)%r_coef
640 res_rho = cmplx(0.0_dp, 0.0_dp, kind=
dp)
642 res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
645 IF (only_kerker)
THEN
647 mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
648 f_mix = alpha*mixing_store%kerker_factor(ig)
649 rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
650 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
651 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
654 IF (mixing_store%gmix_p)
THEN
657 IF (mixing_store%paw(iatom))
THEN
658 mixing_store%cpc_h_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef - &
659 mixing_store%cpc_h_in(iatom, ispin)%r_coef
660 mixing_store%cpc_s_lastres(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef - &
661 mixing_store%cpc_s_in(iatom, ispin)%r_coef
663 rho_atom(iatom)%cpc_h(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef + &
664 mixing_store%cpc_h_in(iatom, ispin)%r_coef*(1._dp - alpha)
665 rho_atom(iatom)%cpc_s(ispin)%r_coef = alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef + &
666 mixing_store%cpc_s_in(iatom, ispin)%r_coef*(1._dp - alpha)
668 mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
669 mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
670 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
671 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
691 mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
692 mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
693 res_norm = res_norm + &
694 REAL(res_rho(ig),
dp)*
REAL(res_rho(ig), dp) + &
695 aimag(res_rho(ig))*aimag(res_rho(ig))
696 delta_norm = delta_norm + &
697 REAL(mixing_store%res_buffer(ib, ispin)%cc(ig),
dp)* &
698 REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
699 aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))* &
700 aimag(mixing_store%res_buffer(ib, ispin)%cc(ig))
701 rho_norm = rho_norm + &
702 REAL(rho_g(ispin)%array(ig),
dp)*
REAL(rho_g(ispin)%array(ig), dp) + &
703 aimag(rho_g(ispin)%array(ig))*aimag(rho_g(ispin)%array(ig))
706 mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
707 mixing_store%rhoin(ispin)%cc(ig) - &
708 mixing_store%rhoin_old(ispin)%cc(ig)
710 CALL para_env%sum(delta_norm)
711 delta_norm = sqrt(delta_norm)
712 CALL para_env%sum(res_norm)
713 res_norm = sqrt(res_norm)
714 CALL para_env%sum(rho_norm)
715 rho_norm = sqrt(rho_norm)
717 IF (res_norm > 1.e-14_dp)
THEN
718 mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm
719 mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm
721 IF (mixing_store%gmix_p .AND. gapw)
THEN
723 IF (mixing_store%paw(iatom))
THEN
724 n1 =
SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
725 n2 =
SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
728 mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
729 (mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
730 mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
731 dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
732 mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
733 mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
734 mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
735 mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
737 mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
738 (mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
739 mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
740 dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
741 mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
742 mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
743 mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
744 mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
746 mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
748 mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
749 mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
751 mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
760 f_mix = alpha*mixing_store%kerker_factor(ig)
761 mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
762 f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
763 mixing_store%drho_buffer(ib, ispin)%cc(ig)
767 DO ig = 1, mixing_store%ig_max
768 a(kb, jb) = a(kb, jb) + mixing_store%p_metric(ig)*( &
769 REAL(mixing_store%res_buffer(jb, ispin)%cc(ig),
dp)* &
770 REAL(mixing_store%res_buffer(kb, ispin)%cc(ig), dp) + &
771 aimag(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
772 aimag(mixing_store%res_buffer(kb, ispin)%cc(ig)))
774 a(jb, kb) = a(kb, jb)
781 a(jb, jb) = w0 + a(jb, jb)
782 DO ig = 1, mixing_store%ig_max
783 c(jb) = c(jb) + mixing_store%p_metric(ig)*( &
784 REAL(mixing_store%res_buffer(jb, ispin)%cc(ig),
dp)*
REAL(res_rho(ig), dp) + &
785 aimag(mixing_store%res_buffer(jb, ispin)%cc(ig))*aimag(res_rho(ig)))
789 CALL invert_matrix(a, b, inv_err)
791 CALL dgemv(
'T', nb, nb, 1.0_dp, b, nb, c, 1, 0.0_dp, g, 1)
797 cc_mix = cmplx(0.0_dp, 0.0_dp, kind=
dp)
799 cc_mix = cc_mix - g(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
801 f_mix = alpha*mixing_store%kerker_factor(ig)
803 IF (res_norm > 1.e-14_dp) rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
804 f_mix*res_rho(ig) + cc_mix
805 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
806 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
809 IF (mixing_store%gmix_p)
THEN
812 IF (mixing_store%paw(iatom))
THEN
813 n1 =
SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
814 n2 =
SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
820 valh = valh - g(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
821 vals = vals - g(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
823 IF (res_norm > 1.e-14_dp)
THEN
824 rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
825 alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
826 mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + valh
827 rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
828 alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
829 mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + vals
834 mixing_store%cpc_h_old(iatom, ispin)%r_coef = mixing_store%cpc_h_in(iatom, ispin)%r_coef
835 mixing_store%cpc_s_old(iatom, ispin)%r_coef = mixing_store%cpc_s_in(iatom, ispin)%r_coef
836 mixing_store%cpc_h_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_h(ispin)%r_coef
837 mixing_store%cpc_s_in(iatom, ispin)%r_coef = rho_atom(iatom)%cpc_s(ispin)%r_coef
843 DEALLOCATE (a, b, c, g)
848 CALL pw_axpy(rho_g(1), rho_g(2), 1.0_dp, -1.0_dp)
850 mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
852 IF (gapw .AND. mixing_store%gmix_p)
THEN
854 IF (mixing_store%paw(iatom))
THEN
855 rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
856 rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
857 mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
858 mixing_store%cpc_h_in(iatom, 2)%r_coef
859 mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
860 mixing_store%cpc_s_in(iatom, 2)%r_coef
869 CALL timestop(handle)
871 END SUBROUTINE broyden_mixing
879 SUBROUTINE broyden_mixing_new(mixing_store, rho, para_env)
881 TYPE(mixing_storage_type),
POINTER :: mixing_store
882 TYPE(qs_rho_type),
POINTER :: rho
883 TYPE(mp_para_env_type),
POINTER :: para_env
885 CHARACTER(len=*),
PARAMETER :: routinen =
'broyden_mixing_new'
887 COMPLEX(dp),
ALLOCATABLE,
DIMENSION(:) :: delta_res_p, res_rho, res_rho_p, tmp
888 INTEGER :: handle, ib, ibb, ig, info, ispin, jb, &
889 kb, kkb, lwork, nb, nbuffer, ng, nspin
890 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: iwork
891 LOGICAL :: only_kerker, skip_bq
892 REAL(
dp) :: alpha, beta, delta, delta_p, delta_rhog, delta_rhog_p, f_mix, imp, imp_j, &
893 inv_err, norm, norm_ig, rep, rep_j, sqt_uvol, sqt_vol, uvol, vol, wc, wmax
894 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: aval, work_dgesdd
895 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: a, au, av, b, bq, work
896 REAL(
dp),
DIMENSION(:),
POINTER :: p_metric
897 REAL(
dp),
DIMENSION(:, :),
POINTER :: fmat, weight
898 TYPE(cp_1d_z_p_type),
DIMENSION(:),
POINTER :: tmp_z
899 TYPE(cp_1d_z_p_type),
DIMENSION(:, :),
POINTER :: delta_res, u_vec, z_vec
900 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g
902 cpassert(
ASSOCIATED(mixing_store%rhoin_buffer))
904 CALL timeset(routinen, handle)
906 NULLIFY (delta_res, u_vec, z_vec)
907 NULLIFY (fmat, rho_g)
910 nspin =
SIZE(rho_g, 1)
911 ng =
SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
912 vol = rho_g(1)%pw_grid%vol
915 sqt_uvol = sqrt(uvol)
916 alpha = mixing_store%alpha
917 beta = mixing_store%beta
920 wmax = mixing_store%wmax
921 nbuffer = mixing_store%nbuffer
923 mixing_store%ncall = mixing_store%ncall + 1
924 IF (mixing_store%ncall == 1)
THEN
927 only_kerker = .false.
928 nb = min(mixing_store%ncall - 1, nbuffer)
929 ib =
modulo(mixing_store%ncall - 2, nbuffer) + 1
935 ALLOCATE (bq(nb, nb))
939 ALLOCATE (delta_res_p(ng))
944 ALLOCATE (res_rho(ng))
945 ALLOCATE (res_rho_p(ng))
947 p_metric => mixing_store%p_metric
948 weight => mixing_store%weight
949 cpassert(
ASSOCIATED(mixing_store%delta_res))
950 delta_res => mixing_store%delta_res
951 cpassert(
ASSOCIATED(mixing_store%u_vec))
952 u_vec => mixing_store%u_vec
953 cpassert(
ASSOCIATED(mixing_store%z_vec))
954 z_vec => mixing_store%z_vec
957 delta_rhog_p = 0.0_dp
961 fmat => mixing_store%fmat(:, :, ispin)
969 res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
970 res_rho_p(ig) = res_rho(ig)*p_metric(ig)
971 norm_ig = real(res_rho(ig),
dp)*real(res_rho(ig),
dp) + aimag(res_rho(ig))*aimag(res_rho(ig))
972 delta = delta + norm_ig
973 delta_p = delta_p + norm_ig*p_metric(ig)
975 CALL para_env%sum(delta)
977 CALL para_env%sum(delta_p)
978 delta_p = sqrt(delta_p)
979 delta_rhog = delta_rhog + delta
980 delta_rhog_p = delta_rhog_p + delta_p
982 weight(ib, ispin) = 1.0_dp
983 IF (wc < 0.0_dp) weight(ib, ispin) = 0.01_dp*abs(wc)/(delta_p*delta_p)
984 IF (weight(ib, ispin) == 0.0_dp) weight(ib, ispin) = 100.0_dp
985 IF (weight(ib, ispin) < 1.0_dp) weight(ib, ispin) = 1.0_dp
987 IF (only_kerker)
THEN
990 f_mix = alpha*mixing_store%kerker_factor(ig)
991 rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
992 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
993 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
994 mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
1000 delta_res(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
1001 delta_res_p(ig) = p_metric(ig)*delta_res(ib, ispin)%cc(ig)
1002 norm_ig = real(delta_res(ib, ispin)%cc(ig),
dp)*real(delta_res_p(ig),
dp) + &
1003 aimag(delta_res(ib, ispin)%cc(ig))*aimag(delta_res_p(ig))
1004 norm = norm + norm_ig
1006 CALL para_env%sum(norm)
1007 norm = 1._dp/sqrt(norm)
1008 delta_res(ib, ispin)%cc(:) = delta_res(ib, ispin)%cc(:)*norm
1009 delta_res_p(:) = delta_res_p(:)*norm
1011 IF (para_env%is_source())
WRITE (*, *)
'norm ', norm
1014 tmp(ig) = (mixing_store%rhoin(ispin)%cc(ig) - &
1015 mixing_store%rhoin_old(ispin)%cc(ig))*norm
1016 u_vec(ib, ispin)%cc(ig) = (tmp(ig) + &
1017 mixing_store%kerker_factor(ig)*delta_res(ib, ispin)%cc(ig))
1021 fmat(jb, ib) = 0.0_dp
1024 rep_j = real(delta_res(jb, ispin)%cc(ig),
dp)
1025 imp_j = aimag(delta_res(jb, ispin)%cc(ig))
1027 rep = real(delta_res_p(ig),
dp)
1028 imp = aimag(delta_res_p(ig))
1029 fmat(jb, ib) = fmat(jb, ib) + rep_j*rep + imp_j*imp
1033 CALL para_env%sum(fmat(1:nb, ib))
1035 fmat(ib, ib) = 1.0_dp
1038 fmat(ib, jb) = fmat(jb, ib)
1042 a(jb, jb) = weight(jb, ispin)*weight(jb, ispin)*fmat(jb, jb)
1044 a(jb, kb) = weight(jb, ispin)*weight(kb, ispin)*fmat(jb, kb)
1045 a(kb, jb) = weight(jb, ispin)*weight(kb, ispin)*fmat(kb, jb)
1050 CALL invert_matrix(a, b, inv_err)
1053 ALLOCATE (work(ib - 1, ib - 1), aval(ib - 1), av(ib - 1, ib - 1), au(ib - 1, ib - 1))
1055 ALLOCATE (iwork(8*(ib - 1)), work_dgesdd(1))
1057 CALL dgesdd(
'S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info)
1058 lwork = int(work_dgesdd(1))
1059 DEALLOCATE (work_dgesdd);
ALLOCATE (work_dgesdd(lwork))
1060 CALL dgesdd(
'S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info)
1064 IF (aval(kb) < 1.e-6_dp)
THEN
1067 aval(kb) = 1.0_dp/aval(kb)
1069 av(kb, :) = av(kb, :)*aval(kb)
1071 CALL dgemm(
'T',
'T', ib - 1, ib - 1, ib - 1, 1.0_dp, av, ib - 1, au, ib - 1, 0.0_dp, b, ib - 1)
1072 DEALLOCATE (iwork, aval, au, av, work_dgesdd, work)
1082 bq(jb, kb) = bq(jb, kb) - weight(jb, ispin)*weight(kkb, ispin)*b(jb, kkb)*fmat(kkb, kb)
1085 bq(jb, jb) = 1.0_dp + bq(jb, jb)
1088 IF (.NOT. skip_bq)
THEN
1091 ALLOCATE (tmp_z(nb))
1093 ALLOCATE (tmp_z(jb)%cc(ng))
1097 tmp(:) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1098 IF (.NOT. skip_bq)
THEN
1102 IF (weight(kb, ispin) >= (10.0_dp*wmax)) cycle
1104 tmp(ig) = tmp(ig) + bq(jb, kb)*z_vec(kb, ispin)%cc(ig)
1112 tmp(ig) = tmp(ig) + weight(kb, ispin)*weight(jb, ispin)*b(jb, kb)*u_vec(kb, ispin)%cc(ig)
1117 IF (skip_bq .OR. (weight(jb, ispin) >= (10._dp*wmax)))
THEN
1118 z_vec(jb, ispin)%cc(:) = tmp(:)
1121 tmp_z(jb)%cc(:) = tmp(:)
1125 IF (.NOT. skip_bq)
THEN
1127 IF (weight(jb, ispin) < (10._dp*wmax)) z_vec(jb, ispin)%cc(:) = tmp_z(jb)%cc(:)
1128 DEALLOCATE (tmp_z(jb)%cc)
1134 rho_g(ispin)%array(:) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1138 rep_j = real(delta_res(jb, ispin)%cc(ig),
dp)
1139 imp_j = aimag(delta_res(jb, ispin)%cc(ig))
1140 rep = real(res_rho_p(ig),
dp)
1141 imp = aimag(res_rho_p(ig))
1142 norm = norm + rep_j*rep + imp_j*imp
1144 CALL para_env%sum(norm)
1147 rho_g(ispin)%array(ig) = rho_g(ispin)%array(ig) - norm*z_vec(jb, ispin)%cc(ig)*sqt_vol
1153 f_mix = alpha*mixing_store%kerker_factor(ig)
1154 rho_g(ispin)%array(ig) = rho_g(ispin)%array(ig) + &
1155 mixing_store%rhoin_buffer(ib, ispin)%cc(ig) + f_mix*res_rho(ig)
1156 mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1159 mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
1160 mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
1161 mixing_store%last_res(ispin)%cc(:) = res_rho(:)
1165 IF (.NOT. only_kerker)
THEN
1166 DEALLOCATE (a, b, bq)
1167 DEALLOCATE (delta_res_p, tmp)
1169 DEALLOCATE (res_rho, res_rho_p)
1171 CALL timestop(handle)
1173 END SUBROUTINE broyden_mixing_new
1188 SUBROUTINE multisecant_mixing(mixing_store, rho, para_env)
1190 TYPE(mixing_storage_type),
POINTER :: mixing_store
1191 TYPE(qs_rho_type),
POINTER :: rho
1192 TYPE(mp_para_env_type),
POINTER :: para_env
1194 CHARACTER(len=*),
PARAMETER :: routinen =
'multisecant_mixing'
1195 COMPLEX(KIND=dp),
PARAMETER :: cmone = (-1.0_dp, 0.0_dp), &
1196 cone = (1.0_dp, 0.0_dp), &
1197 czero = (0.0_dp, 0.0_dp)
1199 COMPLEX(dp) :: saa, yaa
1200 COMPLEX(dp),
ALLOCATABLE,
DIMENSION(:) :: gn_global, pgn, res_matrix_global, &
1202 COMPLEX(dp),
ALLOCATABLE,
DIMENSION(:, :) :: a_matrix, res_matrix, sa, step_matrix, ya
1203 COMPLEX(dp),
DIMENSION(:),
POINTER :: gn
1204 INTEGER :: handle, ib, ib_next, ib_prev, ig, &
1205 ig_global, iig, ispin, jb, kb, nb, &
1206 nbuffer, ng, ng_global, nspin
1207 LOGICAL :: use_zgemm, use_zgemm_rev
1208 REAL(
dp) :: alpha, f_mix, gn_norm, gn_norm_old, inv_err, n_low, n_up, pgn_norm, prec, &
1209 r_step, reg_par, sigma_max, sigma_tilde, step_size
1210 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: norm_res, norm_res_low, norm_res_up
1211 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: b_matrix, binv_matrix
1212 REAL(
dp),
DIMENSION(:),
POINTER :: g2
1213 REAL(
dp),
SAVE :: sigma_old = 1.0_dp
1214 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rho_g
1216 cpassert(
ASSOCIATED(mixing_store))
1218 CALL timeset(routinen, handle)
1223 use_zgemm_rev = .true.
1228 nspin =
SIZE(rho_g, 1)
1230 cpassert(rho_g(1)%pw_grid%ngpts < huge(ng_global))
1231 ng_global = int(rho_g(1)%pw_grid%ngpts)
1232 ng =
SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
1233 alpha = mixing_store%alpha
1235 sigma_max = mixing_store%sigma_max
1236 reg_par = mixing_store%reg_par
1237 r_step = mixing_store%r_step
1238 nbuffer = mixing_store%nbuffer
1241 nb = min(mixing_store%ncall, nbuffer - 1)
1242 ib =
modulo(mixing_store%ncall, nbuffer) + 1
1243 IF (mixing_store%ncall > 0)
THEN
1244 ib_prev =
modulo(mixing_store%ncall - 1, nbuffer) + 1
1248 mixing_store%ncall = mixing_store%ncall + 1
1249 ib_next =
modulo(mixing_store%ncall, nbuffer) + 1
1253 gn => mixing_store%res_buffer(ib, ispin)%cc
1256 gn(ig) = (rho_g(ispin)%array(ig) - mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1257 gn_norm = gn_norm + &
1258 REAL(gn(ig),
dp)*
REAL(gn(ig), dp) + aimag(gn(ig))*aimag(gn(ig))
1260 CALL para_env%sum(gn_norm)
1261 gn_norm = sqrt(gn_norm)
1262 mixing_store%norm_res_buffer(ib, ispin) = gn_norm
1269 f_mix = alpha*mixing_store%kerker_factor(ig)
1270 rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(1, ispin)%cc(ig) + &
1271 f_mix*mixing_store%res_buffer(1, ispin)%cc(ig)
1272 mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1275 CALL timestop(handle)
1281 ALLOCATE (step_matrix(ng, nb))
1283 ALLOCATE (res_matrix(ng, nb))
1285 ALLOCATE (a_matrix(nb, ng_global))
1287 ALLOCATE (norm_res(nb))
1288 ALLOCATE (norm_res_low(nb))
1289 ALLOCATE (norm_res_up(nb))
1292 ALLOCATE (b_matrix(nb, nb))
1294 ALLOCATE (binv_matrix(nb, nb))
1296 ALLOCATE (gn_global(ng_global))
1297 ALLOCATE (res_matrix_global(ng_global))
1299 ALLOCATE (sa(ng, ng_global))
1300 ALLOCATE (ya(ng, ng_global))
1302 IF (use_zgemm_rev)
THEN
1303 ALLOCATE (tmp_vec(nb))
1310 gn => mixing_store%res_buffer(ib, ispin)%cc
1311 gn_global = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1313 ig_global = mixing_store%ig_global_index(ig)
1314 gn_global(ig_global) = gn(ig)
1316 CALL para_env%sum(gn_global)
1323 g2 => rho_g(1)%pw_grid%gsq
1328 ELSEIF (jb > ib)
THEN
1333 norm_res(kb) = 0.0_dp
1334 norm_res_low(kb) = 0.0_dp
1335 norm_res_up(kb) = 0.0_dp
1340 step_matrix(ig, kb) = prec*(mixing_store%rhoin_buffer(jb, ispin)%cc(ig) - &
1341 mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
1342 res_matrix(ig, kb) = (mixing_store%res_buffer(jb, ispin)%cc(ig) - &
1343 mixing_store%res_buffer(ib, ispin)%cc(ig))
1344 norm_res(kb) = norm_res(kb) + real(res_matrix(ig, kb),
dp)*real(res_matrix(ig, kb),
dp) + &
1345 aimag(res_matrix(ig, kb))*aimag(res_matrix(ig, kb))
1346 IF (g2(ig) < 4.0_dp)
THEN
1347 norm_res_low(kb) = norm_res_low(kb) + &
1348 REAL(res_matrix(ig, kb),
dp)*
REAL(res_matrix(ig, kb), dp) + &
1349 aimag(res_matrix(ig, kb))*aimag(res_matrix(ig, kb))
1351 norm_res_up(kb) = norm_res_up(kb) + &
1352 REAL(res_matrix(ig, kb),
dp)*
REAL(res_matrix(ig, kb), dp) + &
1353 aimag(res_matrix(ig, kb))*aimag(res_matrix(ig, kb))
1355 res_matrix(ig, kb) = prec*res_matrix(ig, kb)
1360 CALL para_env%sum(norm_res)
1361 CALL para_env%sum(norm_res_up)
1362 CALL para_env%sum(norm_res_low)
1363 norm_res(1:nb) = 1.0_dp/sqrt(norm_res(1:nb))
1367 n_low = n_low + norm_res_low(jb)/norm_res(jb)
1368 n_up = n_up + norm_res_up(jb)/norm_res(jb)
1371 IF (g2(ig) > 4.0_dp)
THEN
1372 step_matrix(ig, 1:nb) = step_matrix(ig, 1:nb)*sqrt(n_low/n_up)
1373 res_matrix(ig, 1:nb) = res_matrix(ig, 1:nb)*sqrt(n_low/n_up)
1377 step_matrix(1:ng, kb) = step_matrix(1:ng, kb)*norm_res(kb)
1378 res_matrix(1:ng, kb) = res_matrix(1:ng, kb)*norm_res(kb)
1385 b_matrix(kb, jb) = 0.0_dp
1389 b_matrix(kb, jb) = b_matrix(kb, jb) + real(res_matrix(ig, kb)*res_matrix(ig, jb),
dp)
1394 CALL para_env%sum(b_matrix)
1396 b_matrix(jb, jb) = b_matrix(jb, jb) + reg_par
1399 CALL invert_matrix(b_matrix, binv_matrix, inv_err)
1402 a_matrix = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1404 res_matrix_global = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1406 ig_global = mixing_store%ig_global_index(ig)
1407 res_matrix_global(ig_global) = res_matrix(ig, kb)
1409 CALL para_env%sum(res_matrix_global)
1413 ig_global = mixing_store%ig_global_index(ig)
1414 a_matrix(jb, ig_global) = a_matrix(jb, ig_global) + &
1415 binv_matrix(jb, kb)*res_matrix_global(ig_global)
1419 CALL para_env%sum(a_matrix)
1422 gn => mixing_store%res_buffer(ib, ispin)%cc
1427 CALL zgemm(
"N",
"N", ng, ng_global, nb, cmone, step_matrix(1, 1), ng, &
1428 a_matrix(1, 1), nb, czero, sa(1, 1), ng)
1429 CALL zgemm(
"N",
"N", ng, ng_global, nb, cmone, res_matrix(1, 1), ng, &
1430 a_matrix(1, 1), nb, czero, ya(1, 1), ng)
1432 ig_global = mixing_store%ig_global_index(ig)
1433 ya(ig, ig_global) = ya(ig, ig_global) + cmplx(1.0_dp, 0.0_dp, kind=
dp)
1436 CALL zgemv(
"N", ng, ng_global, cone, sa(1, 1), &
1437 ng, gn_global(1), 1, czero, pgn(1), 1)
1438 CALL zgemv(
"N", ng, ng_global, cone, ya(1, 1), &
1439 ng, gn_global(1), 1, czero, ugn(1), 1)
1442 pgn_norm = pgn_norm + real(pgn(ig),
dp)*real(pgn(ig),
dp) + &
1443 aimag(pgn(ig))*aimag(pgn(ig))
1445 CALL para_env%sum(pgn_norm)
1446 ELSEIF (use_zgemm_rev)
THEN
1448 CALL zgemv(
"N", nb, ng_global, cone, a_matrix(1, 1), &
1449 nb, gn_global(1), 1, czero, tmp_vec(1), 1)
1451 CALL zgemv(
"N", ng, nb, cmone, step_matrix(1, 1), ng, &
1452 tmp_vec(1), 1, czero, pgn(1), 1)
1454 CALL zgemv(
"N", ng, nb, cmone, res_matrix(1, 1), ng, &
1455 tmp_vec(1), 1, czero, ugn(1), 1)
1458 pgn_norm = pgn_norm + real(pgn(ig),
dp)*real(pgn(ig),
dp) + &
1459 aimag(pgn(ig))*aimag(pgn(ig))
1460 ugn(ig) = ugn(ig) + gn(ig)
1462 CALL para_env%sum(pgn_norm)
1466 pgn(ig) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1467 ugn(ig) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1468 ig_global = mixing_store%ig_global_index(ig)
1469 DO iig = 1, ng_global
1470 saa = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1471 yaa = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1473 IF (ig_global == iig) yaa = cmplx(1.0_dp, 0.0_dp, kind=
dp)
1476 saa = saa - step_matrix(ig, jb)*a_matrix(jb, iig)
1477 yaa = yaa - res_matrix(ig, jb)*a_matrix(jb, iig)
1479 pgn(ig) = pgn(ig) + saa*gn_global(iig)
1480 ugn(ig) = ugn(ig) + yaa*gn_global(iig)
1484 pgn_norm = pgn_norm + real(pgn(ig),
dp)*real(pgn(ig),
dp) + &
1485 aimag(pgn(ig))*aimag(pgn(ig))
1487 CALL para_env%sum(pgn_norm)
1490 gn_norm = mixing_store%norm_res_buffer(ib, ispin)
1491 gn_norm_old = mixing_store%norm_res_buffer(ib_prev, ispin)
1492 IF (ib_prev /= 0)
THEN
1493 sigma_tilde = sigma_old*max(0.5_dp, min(2.0_dp, gn_norm_old/gn_norm))
1495 sigma_tilde = 0.5_dp
1497 sigma_tilde = 0.1_dp
1499 step_size = min(sigma_tilde, r_step*pgn_norm/gn_norm, sigma_max)
1500 sigma_old = step_size
1504 prec = mixing_store%kerker_factor(ig)
1505 rho_g(ispin)%array(ig) = mixing_store%rhoin_buffer(ib, ispin)%cc(ig) &
1506 - prec*step_size*ugn(ig) + prec*pgn(ig)
1507 mixing_store%rhoin_buffer(ib_next, ispin)%cc(ig) = rho_g(ispin)%array(ig)
1513 DEALLOCATE (step_matrix, res_matrix)
1514 DEALLOCATE (norm_res)
1515 DEALLOCATE (norm_res_up)
1516 DEALLOCATE (norm_res_low)
1517 DEALLOCATE (a_matrix, b_matrix, binv_matrix)
1518 DEALLOCATE (ugn, pgn)
1522 IF (use_zgemm_rev)
THEN
1523 DEALLOCATE (tmp_vec)
1525 DEALLOCATE (gn_global, res_matrix_global)
1527 CALL timestop(handle)
1529 END SUBROUTINE multisecant_mixing
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Defines the basic variable types.
integer, parameter, public dp
Collection of simple mathematical functions and subroutines.
Interface to the message passing library MPI.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
module that contains the definitions of the scf types
integer, parameter, public broyden_mixing_nr
integer, parameter, public multisecant_mixing_nr
integer, parameter, public pulay_mixing_nr
integer, parameter, public gspace_mixing_nr
integer, parameter, public broyden_mixing_new_nr
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
Driver for the g-space mixing, calls the proper routine given the requested method.
superstucture that hold various representations of the density and keeps track of which ones are vali...
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...