Skip to content

Commit

Permalink
Remove BROYDEN_MIXING_NEW option (cp2k#3346)
Browse files Browse the repository at this point in the history
  • Loading branch information
@fstein93
fstein93 authored Apr 5, 2024
1 parent 3204c26 commit 2ed94b8
Show file tree
Hide file tree
Showing 7 changed files with 17 additions and 387 deletions.
5 changes: 1 addition & 4 deletions src/qs_charge_mixing.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,7 @@ MODULE qs_charge_mixing
USE kinds, ONLY: dp
USE mathlib, ONLY: get_pseudo_inverse_svd
USE message_passing, ONLY: mp_para_env_type
USE qs_density_mixing_types, ONLY: broyden_mixing_new_nr,&
broyden_mixing_nr,&
USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
gspace_mixing_nr,&
mixing_storage_type,&
multisecant_mixing_nr,&
Expand Down Expand Up @@ -89,8 +88,6 @@ SUBROUTINE charge_mixing(mixing_method, mixing_store, charges, para_env, iter_co
ELSEIF (mixing_method == broyden_mixing_nr) THEN
CALL broyden_mixing(mixing_store, charges, imin, nvec, ns, para_env)
mixing_store%iter_method = "Broy."
ELSEIF (mixing_method == broyden_mixing_new_nr) THEN
CPABORT("Broyden_mixing_new method not available for Charge Mixing")
ELSEIF (mixing_method == multisecant_mixing_nr) THEN
CPABORT("Multisecant_mixing method not available for Charge Mixing")
END IF
Expand Down
12 changes: 2 additions & 10 deletions src/qs_density_mixing_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
MODULE qs_density_mixing_types
USE ao_util, ONLY: exp_radius
USE input_constants, ONLY: broy_mix,&
broy_mix_new,&
direct_p_mix,&
gaussian,&
kerker_mix,&
Expand Down Expand Up @@ -46,7 +45,6 @@ MODULE qs_density_mixing_types
INTEGER, PARAMETER, PUBLIC :: no_mixing_nr = 0, direct_mixing_nr = 1, &
gspace_mixing_nr = 2, pulay_mixing_nr = 3, &
broyden_mixing_nr = 4, &
broyden_mixing_new_nr = 5, &
multisecant_mixing_nr = 6
PUBLIC :: cp_1d_z_p_type, mixing_storage_create, mixing_storage_type, mixing_storage_release, create_mixing_section

Expand Down Expand Up @@ -188,11 +186,6 @@ SUBROUTINE mixing_storage_create(mixing_store, mixing_section, mixing_method, ec
CASE (broyden_mixing_nr)
CALL section_vals_val_get(mixing_section, "BROY_W0", r_val=mixing_store%broy_w0)
mixing_store%bconst = 20.0_dp
CASE (broyden_mixing_new_nr)
CALL section_vals_val_get(mixing_section, "BROY_WREF", r_val=mixing_store%wc)
CALL section_vals_val_get(mixing_section, "BROY_WMAX", r_val=mixing_store%wmax)
mixing_store%bconst = 20.0_dp
mixing_store%p_metric_method = 1
CASE (multisecant_mixing_nr)
CALL section_vals_val_get(mixing_section, "REGULARIZATION", r_val=mixing_store%reg_par)
CALL section_vals_val_get(mixing_section, "MAX_STEP", r_val=mixing_store%sigma_max)
Expand Down Expand Up @@ -507,14 +500,13 @@ SUBROUTINE create_mixing_section(section, ls_scf)
"KERKER_MIXING", &
"PULAY_MIXING", &
"BROYDEN_MIXING", &
"BROYDEN_MIXING_NEW", &
"MULTISECANT_MIXING"), &
enum_i_vals=(/no_mix, direct_p_mix, kerker_mix, pulay_mix, broy_mix, &
broy_mix_new, multisec_mix/), &
multisec_mix/), &
enum_desc=s2a("No mixing is applied", &
"Direct mixing of new and old density matrices", &
"Mixing of the potential in reciprocal space using the Kerker damping", &
"Pulay mixing", "Broyden mixing", "Broyden mixing second version", &
"Pulay mixing", "Broyden mixing", &
"Multisecant scheme for mixing"))

CALL section_add_keyword(section, keyword)
Expand Down
310 changes: 1 addition & 309 deletions src/qs_gspace_mixing.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,9 +23,7 @@ MODULE qs_gspace_mixing
USE pw_pool_types, ONLY: pw_pool_type
USE pw_types, ONLY: pw_c1d_gs_type,&
pw_r3d_rs_type
USE qs_density_mixing_types, ONLY: broyden_mixing_new_nr,&
broyden_mixing_nr,&
cp_1d_z_p_type,&
USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
gspace_mixing_nr,&
mixing_storage_type,&
multisecant_mixing_nr,&
Expand Down Expand Up @@ -157,10 +155,6 @@ SUBROUTINE gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, ite
ELSEIF (mixing_method == broyden_mixing_nr) THEN
CALL broyden_mixing(qs_env, mixing_store, rho, para_env)
mixing_store%iter_method = "Broy."
ELSEIF (mixing_method == broyden_mixing_new_nr) THEN
CPASSERT(.NOT. gapw)
CALL broyden_mixing_new(mixing_store, rho, para_env)
mixing_store%iter_method = "Broy."
ELSEIF (mixing_method == multisecant_mixing_nr) THEN
CPASSERT(.NOT. gapw)
CALL multisecant_mixing(mixing_store, rho, para_env)
Expand Down Expand Up @@ -870,308 +864,6 @@ SUBROUTINE broyden_mixing(qs_env, mixing_store, rho, para_env)

END SUBROUTINE broyden_mixing

! **************************************************************************************************
!> \brief ...
!> \param mixing_store ...
!> \param rho ...
!> \param para_env ...
! **************************************************************************************************
SUBROUTINE broyden_mixing_new(mixing_store, rho, para_env)

TYPE(mixing_storage_type), POINTER :: mixing_store
TYPE(qs_rho_type), POINTER :: rho
TYPE(mp_para_env_type), POINTER :: para_env

CHARACTER(len=*), PARAMETER :: routineN = 'broyden_mixing_new'

COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: delta_res_p, res_rho, res_rho_p, tmp
INTEGER :: handle, ib, ibb, ig, info, ispin, jb, &
kb, kkb, lwork, nb, nbuffer, ng, nspin
INTEGER, ALLOCATABLE, DIMENSION(:) :: IWORK
LOGICAL :: only_kerker, skip_bq
REAL(dp) :: alpha, beta, delta, delta_p, delta_rhog, delta_rhog_p, f_mix, imp, imp_j, &
inv_err, norm, norm_ig, rep, rep_j, sqt_uvol, sqt_vol, uvol, vol, wc, wmax
REAL(dp), ALLOCATABLE, DIMENSION(:) :: aval, work_dgesdd
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, au, av, b, bq, work
REAL(dp), DIMENSION(:), POINTER :: p_metric
REAL(dp), DIMENSION(:, :), POINTER :: fmat, weight
TYPE(cp_1d_z_p_type), DIMENSION(:), POINTER :: tmp_z
TYPE(cp_1d_z_p_type), DIMENSION(:, :), POINTER :: delta_res, u_vec, z_vec
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g

CPASSERT(ASSOCIATED(mixing_store%rhoin_buffer))

CALL timeset(routineN, handle)

NULLIFY (delta_res, u_vec, z_vec)
NULLIFY (fmat, rho_g)
CALL qs_rho_get(rho, rho_g=rho_g)

nspin = SIZE(rho_g, 1)
ng = SIZE(mixing_store%rhoin_buffer(1, 1)%cc)
vol = rho_g(1)%pw_grid%vol
sqt_vol = SQRT(vol)
uvol = 1.0_dp/vol
sqt_uvol = SQRT(uvol)
alpha = mixing_store%alpha
beta = mixing_store%beta

wc = mixing_store%wc
wmax = mixing_store%wmax
nbuffer = mixing_store%nbuffer

mixing_store%ncall = mixing_store%ncall + 1
IF (mixing_store%ncall == 1) THEN
only_kerker = .TRUE.
ELSE
only_kerker = .FALSE.
nb = MIN(mixing_store%ncall - 1, nbuffer)
ib = MODULO(mixing_store%ncall - 2, nbuffer) + 1

ALLOCATE (a(nb, nb))
a = 0.0_dp
ALLOCATE (b(nb, nb))
b = 0.0_dp
ALLOCATE (bq(nb, nb))
bq = 0.0_dp

ALLOCATE (tmp(ng))
ALLOCATE (delta_res_p(ng))
END IF

ibb = 0

ALLOCATE (res_rho(ng))
ALLOCATE (res_rho_p(ng))

p_metric => mixing_store%p_metric
weight => mixing_store%weight
CPASSERT(ASSOCIATED(mixing_store%delta_res))
delta_res => mixing_store%delta_res
CPASSERT(ASSOCIATED(mixing_store%u_vec))
u_vec => mixing_store%u_vec
CPASSERT(ASSOCIATED(mixing_store%z_vec))
z_vec => mixing_store%z_vec

delta_rhog = 0.0_dp
delta_rhog_p = 0.0_dp

DO ispin = 1, nspin

fmat => mixing_store%fmat(:, :, ispin)

delta = 0.0_dp
delta_p = 0.0_dp
! Residual at this step R_i(G) (rho_out(G)-rho_in(G))
! Residual multiplied by the metrics RP_i(G) = (rho_out(G)-rho_in(G)) * P(G)
! Delta is the norm of the residual, measures how far we are from convergence
DO ig = 1, ng
res_rho(ig) = rho_g(ispin)%array(ig) - mixing_store%rhoin(ispin)%cc(ig)
res_rho_p(ig) = res_rho(ig)*p_metric(ig) !*sqt_uvol
norm_ig = REAL(res_rho(ig), dp)*REAL(res_rho(ig), dp) + AIMAG(res_rho(ig))*AIMAG(res_rho(ig))
delta = delta + norm_ig
delta_p = delta_p + norm_ig*p_metric(ig) !*p_metric(ig)
END DO
CALL para_env%sum(delta)
delta = SQRT(delta)
CALL para_env%sum(delta_p)
delta_p = SQRT(delta_p)
delta_rhog = delta_rhog + delta
delta_rhog_p = delta_rhog_p + delta_p

weight(ib, ispin) = 1.0_dp ! wc
IF (wc < 0.0_dp) weight(ib, ispin) = 0.01_dp*ABS(wc)/(delta_p*delta_p)
IF (weight(ib, ispin) == 0.0_dp) weight(ib, ispin) = 100.0_dp
IF (weight(ib, ispin) < 1.0_dp) weight(ib, ispin) = 1.0_dp

IF (only_kerker) THEN
! Simple Kerker damping : linear mixing rho(G) = rho_in(G) - alpha k(G)*(rho_out(G)-rho_in(G))
DO ig = 1, ng
f_mix = alpha*mixing_store%kerker_factor(ig)
rho_g(ispin)%array(ig) = mixing_store%rhoin(ispin)%cc(ig) + f_mix*res_rho(ig)
mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
END DO
ELSE
norm = 0.0_dp
! Difference of residuals DR_{i-1)} (G) = R_i(G) - R_{i-1}(G)
DO ig = 1, ng
delta_res(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
delta_res_p(ig) = p_metric(ig)*delta_res(ib, ispin)%cc(ig)
norm_ig = REAL(delta_res(ib, ispin)%cc(ig), dp)*REAL(delta_res_p(ig), dp) + &
AIMAG(delta_res(ib, ispin)%cc(ig))*AIMAG(delta_res_p(ig))
norm = norm + norm_ig
END DO
CALL para_env%sum(norm)
norm = 1._dp/SQRT(norm)
delta_res(ib, ispin)%cc(:) = delta_res(ib, ispin)%cc(:)*norm
delta_res_p(:) = delta_res_p(:)*norm

IF (para_env%is_source()) WRITE (*, *) 'norm ', norm
! Vector U_{i-1}(G) = Drho_{i-1} + k(G) * DR_{i-1}(G)
DO ig = 1, ng
tmp(ig) = (mixing_store%rhoin(ispin)%cc(ig) - &
mixing_store%rhoin_old(ispin)%cc(ig))*norm
u_vec(ib, ispin)%cc(ig) = (tmp(ig) + &
mixing_store%kerker_factor(ig)*delta_res(ib, ispin)%cc(ig))
END DO

DO jb = 1, nb
fmat(jb, ib) = 0.0_dp

DO ig = 1, ng
rep_j = REAL(delta_res(jb, ispin)%cc(ig), dp)
imp_j = AIMAG(delta_res(jb, ispin)%cc(ig))
! < DR_{j} | DR_{i-1} >
rep = REAL(delta_res_p(ig), dp)
imp = AIMAG(delta_res_p(ig))
fmat(jb, ib) = fmat(jb, ib) + rep_j*rep + imp_j*imp
END DO

END DO
CALL para_env%sum(fmat(1:nb, ib))

fmat(ib, ib) = 1.0_dp

DO jb = 1, nb
fmat(ib, jb) = fmat(jb, ib)
END DO

DO jb = 1, nb
a(jb, jb) = weight(jb, ispin)*weight(jb, ispin)*fmat(jb, jb)
DO kb = 1, jb - 1
a(jb, kb) = weight(jb, ispin)*weight(kb, ispin)*fmat(jb, kb)
a(kb, jb) = weight(jb, ispin)*weight(kb, ispin)*fmat(kb, jb)
END DO
END DO
IF (.TRUE.) THEN
b = 0.0_dp
CALL invert_matrix(a, b, inv_err)
ELSE
b = 0.0_dp
ALLOCATE (work(ib - 1, ib - 1), aval(ib - 1), av(ib - 1, ib - 1), au(ib - 1, ib - 1))
work(:, :) = a
ALLOCATE (iwork(8*(ib - 1)), work_dgesdd(1))
lwork = -1
CALL DGESDD('S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info)
lwork = INT(work_dgesdd(1))
DEALLOCATE (work_dgesdd); ALLOCATE (work_dgesdd(lwork))
CALL DGESDD('S', ib - 1, ib - 1, work, ib - 1, aval, au, ib - 1, av, ib - 1, work_dgesdd, lwork, iwork, info)
! construct the inverse
DO kb = 1, ib - 1
! invert SV
IF (aval(kb) < 1.E-6_dp) THEN
aval(kb) = 0.0_dp
ELSE
aval(kb) = 1.0_dp/aval(kb)
END IF
av(kb, :) = av(kb, :)*aval(kb)
END DO
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)
DEALLOCATE (iwork, aval, au, av, work_dgesdd, work)
END IF

bq = 0.0_dp
! Broyden second method requires also bq (NYI)
skip_bq = .TRUE.
DO jb = 1, nb
DO kb = 1, nb
bq(jb, kb) = 0.0_dp
DO kkb = 1, nb
bq(jb, kb) = bq(jb, kb) - weight(jb, ispin)*weight(kkb, ispin)*b(jb, kkb)*fmat(kkb, kb)
END DO
END DO
bq(jb, jb) = 1.0_dp + bq(jb, jb)
END DO

IF (.NOT. skip_bq) THEN
! in this case the old z_vec is needed
! a temporary array is needed to store the new one
ALLOCATE (tmp_z(nb))
DO jb = 1, nb
ALLOCATE (tmp_z(jb)%cc(ng))
END DO
END IF
DO jb = 1, nb
tmp(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
IF (.NOT. skip_bq) THEN
! sum_{kb} bq(jb,kb) * z_vec_{kb,iter-2}
! added only for small weights
DO kb = 1, nb
IF (weight(kb, ispin) >= (10.0_dp*wmax)) CYCLE
DO ig = 1, ng
tmp(ig) = tmp(ig) + bq(jb, kb)*z_vec(kb, ispin)%cc(ig)
END DO
END DO ! kb
END IF

! sum_{kb} w(jb)*w(kb)*b(jb,kb) * u_vec_{kb}
DO kb = 1, ib - 1
DO ig = 1, ng
tmp(ig) = tmp(ig) + weight(kb, ispin)*weight(jb, ispin)*b(jb, kb)*u_vec(kb, ispin)%cc(ig)
END DO
END DO

! store the new z_vec(jb)
IF (skip_bq .OR. (weight(jb, ispin) >= (10._dp*wmax))) THEN
z_vec(jb, ispin)%cc(:) = tmp(:)
ELSE
! temporary array: old z_vec may still be needed
tmp_z(jb)%cc(:) = tmp(:)
END IF
END DO !jb

IF (.NOT. skip_bq) THEN
DO jb = 1, ib - 1
IF (weight(jb, ispin) < (10._dp*wmax)) z_vec(jb, ispin)%cc(:) = tmp_z(jb)%cc(:)
DEALLOCATE (tmp_z(jb)%cc)
END DO
DEALLOCATE (tmp_z)
END IF

! Overwrite the density i reciprocal space
rho_g(ispin)%array(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
DO jb = 1, ib - 1
norm = 0.0_dp
DO ig = 1, ng
rep_j = REAL(delta_res(jb, ispin)%cc(ig), dp)
imp_j = AIMAG(delta_res(jb, ispin)%cc(ig))
rep = REAL(res_rho_p(ig), dp)
imp = AIMAG(res_rho_p(ig))
norm = norm + rep_j*rep + imp_j*imp
END DO
CALL para_env%sum(norm)
! Subtract |Z_jb)><DR_jb|P|R_{iter}>
DO ig = 1, ng
rho_g(ispin)%array(ig) = rho_g(ispin)%array(ig) - norm*z_vec(jb, ispin)%cc(ig)*sqt_vol

END DO
END DO

DO ig = 1, ng
f_mix = alpha*mixing_store%kerker_factor(ig)
rho_g(ispin)%array(ig) = rho_g(ispin)%array(ig) + &
mixing_store%rhoin_buffer(ib, ispin)%cc(ig) + f_mix*res_rho(ig)
mixing_store%rhoin_buffer(ibb, ispin)%cc(ig) = rho_g(ispin)%array(ig)
END DO

mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%array(ig)
mixing_store%last_res(ispin)%cc(:) = res_rho(:)
END IF ! ib

END DO ! ispin
IF (.NOT. only_kerker) THEN
DEALLOCATE (a, b, bq)
DEALLOCATE (delta_res_p, tmp)
END IF
DEALLOCATE (res_rho, res_rho_p)

CALL timestop(handle)

END SUBROUTINE broyden_mixing_new

! **************************************************************************************************
!> \brief Multisecant scheme to perform charge density Mixing
!> as preconditioner we use the Kerker damping factor
Expand Down
Loading

0 comments on commit 2ed94b8

Please sign in to comment.