LSD and G-space mixing (#1719)

* gspace mixing with lsd

* prettify

* GMIX_P is needed when GAPW and G-space mixing

* check that GMIX_P is active for mixing the gapw densities

* g-space mixing with spin polarized uses spin up and spin down densities
 rather than spin up density and spin up minus spin down density

* pretty

* changed tolerance

* reference changed
This commit is contained in:
marcella 2021-11-02 12:23:58 +01:00 committed by GitHub
parent e840aa2a81
commit e6ceced390
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 236 additions and 120 deletions

View file

@ -56,7 +56,7 @@ MODULE qs_density_mixing_types
TYPE mixing_storage_type
INTEGER :: ref_count
INTEGER :: ig_max, ncall, nbuffer, n_simple_mix, &
INTEGER :: ig_max, ncall, ncall_p(2), nbuffer, n_simple_mix, &
nskip_mixing, p_metric_method
INTEGER, POINTER, DIMENSION(:) :: ig_global_index
LOGICAL :: gmix_p
@ -70,8 +70,8 @@ MODULE qs_density_mixing_types
REAL(KIND=dp), DIMENSION(:), POINTER :: kerker_factor
REAL(KIND=dp), DIMENSION(:), POINTER :: special_metric
REAL(KIND=dp), DIMENSION(:, :), POINTER :: weight
REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_res_buffer, pulay_matrix
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: fmat, gmat, smat
REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_res_buffer
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: fmat, gmat, pulay_matrix, smat
!
INTEGER :: nat_local, max_shell
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: acharge
@ -121,6 +121,7 @@ CONTAINS
mixing_store%nbuffer = 0
mixing_store%n_simple_mix = 0
mixing_store%ncall = 0
mixing_store%ncall_p = 0
mixing_store%alpha = 1.0_dp
mixing_store%pulay_beta = 1.0_dp
mixing_store%beta = 1.0_dp

View file

@ -300,7 +300,8 @@ CONTAINS
ng, nspin
LOGICAL :: gapw
REAL(dp) :: alpha_kerker, alpha_pulay, beta, f_mix, &
inv_err, norm, norm_c_inv, vol
inv_err, norm, norm_c_inv, res_norm, &
vol
REAL(dp), ALLOCATABLE, DIMENSION(:) :: alpha_c, ev
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, c, c_inv, cpc_h_mix, cpc_s_mix
TYPE(dft_control_type), POINTER :: dft_control
@ -328,52 +329,86 @@ CONTAINS
natom = SIZE(rho_atom)
END IF
ib = MODULO(mixing_store%ncall, nbuffer) + 1
mixing_store%ncall = mixing_store%ncall + 1
nb = MIN(mixing_store%ncall, nbuffer)
ibb = MODULO(mixing_store%ncall, nbuffer) + 1
nb1 = nb + 1
ALLOCATE (a(nb1, nb1))
a = 0.0_dp
ALLOCATE (b(nb1, nb1))
b = 0.0_dp
ALLOCATE (c(nb, nb))
c = 0.0_dp
ALLOCATE (c_inv(nb, nb))
c_inv = 0.0_dp
ALLOCATE (alpha_c(nb))
alpha_c = 0.0_dp
norm_c_inv = 0.0_dp
ALLOCATE (ev(nb1))
ev = 0.0_dp
ALLOCATE (cc_mix(ng))
IF (nspin == 2) THEN
DO ig = 1, ng
rho_g(2)%pw%cc(ig) = rho_g(1)%pw%cc(ig) - rho_g(2)%pw%cc(ig)
mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
END DO
IF (gapw .AND. mixing_store%gmix_p) THEN
DO iatom = 1, natom
IF (mixing_store%paw(iatom)) THEN
rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
mixing_store%cpc_h_in(iatom, 2)%r_coef
mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
mixing_store%cpc_s_in(iatom, 2)%r_coef
END IF
END DO
END IF
ENDIF
DO ispin = 1, nspin
ib = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) + 1
nb = MIN(mixing_store%ncall_p(ispin), nbuffer)
ibb = MODULO(mixing_store%ncall_p(ispin), nbuffer) + 1
mixing_store%res_buffer(ib, ispin)%cc(:) = CMPLX(0._dp, 0._dp, KIND=dp)
norm = 0.0_dp
IF (nb == 1) mixing_store%rhoin_buffer(1, ispin)%cc = mixing_store%rhoin(ispin)%cc
res_norm = 0.0_dp
DO ig = 1, ng
f_mix = mixing_store%kerker_factor(ig)
mixing_store%res_buffer(ib, ispin)%cc(ig) = f_mix*(rho_g(ispin)%pw%cc(ig) - &
mixing_store%rhoin_buffer(ib, ispin)%cc(ig))
res_norm = res_norm + &
REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp)*REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))*AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
END DO
CALL mp_sum(res_norm, para_env%group)
res_norm = SQRT(res_norm)
IF (res_norm < 1.E-14_dp) THEN
mixing_store%ncall_p(ispin) = mixing_store%ncall_p(ispin) - 1
ELSE
nb1 = nb + 1
ALLOCATE (a(nb1, nb1))
a = 0.0_dp
ALLOCATE (b(nb1, nb1))
b = 0.0_dp
ALLOCATE (c(nb, nb))
c = 0.0_dp
ALLOCATE (c_inv(nb, nb))
c_inv = 0.0_dp
ALLOCATE (alpha_c(nb))
alpha_c = 0.0_dp
norm_c_inv = 0.0_dp
ALLOCATE (ev(nb1))
ev = 0.0_dp
ENDIF
DO jb = 1, nb
mixing_store%pulay_matrix(jb, ib) = 0.0_dp
mixing_store%pulay_matrix(jb, ib, ispin) = 0.0_dp
DO ig = 1, ng
f_mix = mixing_store%special_metric(ig)
mixing_store%pulay_matrix(jb, ib) = mixing_store%pulay_matrix(jb, ib) + &
f_mix*(REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp) &
*REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig)))
mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin) + &
f_mix*(REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp) &
*REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig)))
END DO
CALL mp_sum(mixing_store%pulay_matrix(jb, ib), para_env%group)
mixing_store%pulay_matrix(jb, ib) = mixing_store%pulay_matrix(jb, ib)*vol*vol
mixing_store%pulay_matrix(ib, jb) = mixing_store%pulay_matrix(jb, ib)
CALL mp_sum(mixing_store%pulay_matrix(jb, ib, ispin), para_env%group)
mixing_store%pulay_matrix(jb, ib, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)*vol*vol
mixing_store%pulay_matrix(ib, jb, ispin) = mixing_store%pulay_matrix(jb, ib, ispin)
END DO
IF (nb == 1) THEN
IF (nb == 1 .OR. res_norm < 1.E-14_dp) THEN
DO ig = 1, ng
f_mix = alpha_kerker*mixing_store%kerker_factor(ig)
cc_mix(ig) = rho_g(ispin)%pw%cc(ig) - &
@ -406,13 +441,14 @@ CONTAINS
END IF
END IF
ELSE
b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb)
b(1:nb, 1:nb) = mixing_store%pulay_matrix(1:nb, 1:nb, ispin)
c(1:nb, 1:nb) = b(1:nb, 1:nb)
b(nb1, 1:nb) = -1.0_dp
b(1:nb, nb1) = -1.0_dp
b(nb1, nb1) = 0.0_dp
CALL invert_matrix(c, c_inv, inv_err)
CALL invert_matrix(c, c_inv, inv_err, improve=.TRUE.)
alpha_c = 0.0_dp
DO i = 1, nb
DO jb = 1, nb
@ -479,13 +515,37 @@ CONTAINS
END DO
END IF
END IF ! nb==1
IF (res_norm > 1.E-14_dp) THEN
DEALLOCATE (a)
DEALLOCATE (b)
DEALLOCATE (ev)
DEALLOCATE (c, c_inv, alpha_c)
ENDIF
END DO ! ispin
DEALLOCATE (a)
DEALLOCATE (b)
DEALLOCATE (ev)
DEALLOCATE (cc_mix)
DEALLOCATE (c, c_inv, alpha_c)
IF (nspin == 2) THEN
DO ig = 1, ng
rho_g(2)%pw%cc(ig) = rho_g(1)%pw%cc(ig) - rho_g(2)%pw%cc(ig)
mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
END DO
IF (gapw .AND. mixing_store%gmix_p) THEN
DO iatom = 1, natom
IF (mixing_store%paw(iatom)) THEN
rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
mixing_store%cpc_h_in(iatom, 2)%r_coef
mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
mixing_store%cpc_s_in(iatom, 2)%r_coef
END IF
END DO
END IF
ENDIF
CALL timestop(handle)
@ -520,7 +580,7 @@ CONTAINS
LOGICAL :: gapw, only_kerker
REAL(dp) :: alpha, dcpc_h_res, dcpc_s_res, &
delta_norm, f_mix, inv_err, res_norm, &
valh, vals, w0
rho_norm, valh, vals, w0
REAL(dp), ALLOCATABLE, DIMENSION(:) :: c, g
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a, b
TYPE(dft_control_type), POINTER :: dft_control
@ -563,6 +623,25 @@ CONTAINS
natom = 0
END IF
IF (nspin == 2) THEN
DO ig = 1, ng
rho_g(2)%pw%cc(ig) = rho_g(1)%pw%cc(ig) - rho_g(2)%pw%cc(ig)
mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
END DO
IF (gapw .AND. mixing_store%gmix_p) THEN
DO iatom = 1, natom
IF (mixing_store%paw(iatom)) THEN
rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
mixing_store%cpc_h_in(iatom, 2)%r_coef
mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
mixing_store%cpc_s_in(iatom, 2)%r_coef
END IF
END DO
END IF
ENDIF
DO ispin = 1, nspin
res_rho = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
DO ig = 1, ng
@ -613,6 +692,7 @@ CONTAINS
delta_norm = 0.0_dp
res_norm = 0.0_dp
rho_norm = 0.0_dp
DO ig = 1, ng
mixing_store%res_buffer(ib, ispin)%cc(ig) = res_rho(ig) - mixing_store%last_res(ispin)%cc(ig)
mixing_store%last_res(ispin)%cc(ig) = res_rho(ig)
@ -624,6 +704,9 @@ CONTAINS
REAL(mixing_store%res_buffer(ib, ispin)%cc(ig), dp) + &
AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))* &
AIMAG(mixing_store%res_buffer(ib, ispin)%cc(ig))
rho_norm = rho_norm + &
REAL(rho_g(ispin)%pw%cc(ig), dp)*REAL(rho_g(ispin)%pw%cc(ig), dp) + &
AIMAG(rho_g(ispin)%pw%cc(ig))*AIMAG(rho_g(ispin)%pw%cc(ig))
END DO
DO ig = 1, ng
mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
@ -634,80 +717,87 @@ CONTAINS
delta_norm = SQRT(delta_norm)
CALL mp_sum(res_norm, para_env%group)
res_norm = SQRT(res_norm)
mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm
mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm
CALL mp_sum(rho_norm, para_env%group)
rho_norm = SQRT(rho_norm)
IF (mixing_store%gmix_p .AND. gapw) THEN
DO iatom = 1, natom
IF (mixing_store%paw(iatom)) THEN
n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
DO i = 1, n2
DO j = 1, n1
mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
(mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
IF (res_norm > 1.E-14_dp) THEN
mixing_store%res_buffer(ib, ispin)%cc(:) = mixing_store%res_buffer(ib, ispin)%cc(:)/delta_norm
mixing_store%drho_buffer(ib, ispin)%cc(:) = mixing_store%drho_buffer(ib, ispin)%cc(:)/delta_norm
mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
(mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
IF (mixing_store%gmix_p .AND. gapw) THEN
DO iatom = 1, natom
IF (mixing_store%paw(iatom)) THEN
n1 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 1)
n2 = SIZE(mixing_store%cpc_s_in(iatom, ispin)%r_coef, 2)
DO i = 1, n2
DO j = 1, n1
mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
(mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i) - &
mixing_store%cpc_h_old(iatom, ispin)%r_coef(j, i))/delta_norm
dcpc_h_res = ((rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)) - &
mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
mixing_store%cpc_h_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) - &
mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)
mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
alpha*dcpc_h_res + &
mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
alpha*dcpc_s_res + &
mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
(mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i) - &
mixing_store%cpc_s_old(iatom, ispin)%r_coef(j, i))/delta_norm
dcpc_s_res = ((rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)) - &
mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i))/delta_norm
mixing_store%cpc_s_lastres(iatom, ispin)%r_coef(j, i) = rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) - &
mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)
mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i) = &
alpha*dcpc_h_res + &
mixing_store%dcpc_h_in(ib, iatom, ispin)%r_coef(j, i)
mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i) = &
alpha*dcpc_s_res + &
mixing_store%dcpc_s_in(ib, iatom, ispin)%r_coef(j, i)
END DO
END DO
END DO
END IF
END DO
END IF
a(:, :) = 0.0_dp
DO ig = 1, ng
f_mix = alpha*mixing_store%kerker_factor(ig)
mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
mixing_store%drho_buffer(ib, ispin)%cc(ig)
END DO
DO jb = 1, nb
DO kb = jb, nb
DO ig = 1, mixing_store%ig_max !ng
a(kb, jb) = a(kb, jb) + mixing_store%p_metric(ig)*( &
REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)* &
REAL(mixing_store%res_buffer(kb, ispin)%cc(ig), dp) + &
AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
AIMAG(mixing_store%res_buffer(kb, ispin)%cc(ig)))
END IF
END DO
a(jb, kb) = a(kb, jb)
END DO
END DO
CALL mp_sum(a, para_env%group)
END IF
C = 0.0_dp
DO jb = 1, nb
a(jb, jb) = w0 + a(jb, jb)
DO ig = 1, mixing_store%ig_max !ng
c(jb) = c(jb) + mixing_store%p_metric(ig)*( &
REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)*REAL(res_rho(ig), dp) + &
AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))*AIMAG(res_rho(ig)))
a(:, :) = 0.0_dp
DO ig = 1, ng
f_mix = alpha*mixing_store%kerker_factor(ig)
mixing_store%drho_buffer(ib, ispin)%cc(ig) = &
f_mix*mixing_store%res_buffer(ib, ispin)%cc(ig) + &
mixing_store%drho_buffer(ib, ispin)%cc(ig)
END DO
END DO
CALL mp_sum(c, para_env%group)
CALL invert_matrix(a, b, inv_err)
DO jb = 1, nb
DO kb = jb, nb
DO ig = 1, mixing_store%ig_max !ng
a(kb, jb) = a(kb, jb) + mixing_store%p_metric(ig)*( &
REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)* &
REAL(mixing_store%res_buffer(kb, ispin)%cc(ig), dp) + &
AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))* &
AIMAG(mixing_store%res_buffer(kb, ispin)%cc(ig)))
END DO
a(jb, kb) = a(kb, jb)
END DO
END DO
CALL mp_sum(a, para_env%group)
CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1)
C = 0.0_dp
DO jb = 1, nb
a(jb, jb) = w0 + a(jb, jb)
DO ig = 1, mixing_store%ig_max !ng
c(jb) = c(jb) + mixing_store%p_metric(ig)*( &
REAL(mixing_store%res_buffer(jb, ispin)%cc(ig), dp)*REAL(res_rho(ig), dp) + &
AIMAG(mixing_store%res_buffer(jb, ispin)%cc(ig))*AIMAG(res_rho(ig)))
END DO
END DO
CALL mp_sum(c, para_env%group)
CALL invert_matrix(a, b, inv_err)
CALL dgemv('T', nb, nb, 1.0_dp, B, nb, C, 1, 0.0_dp, G, 1)
ELSE
G = 0.0_dp
ENDIF
DO ig = 1, ng
cc_mix = CMPLX(0.0_dp, 0.0_dp, kind=dp)
@ -715,8 +805,9 @@ CONTAINS
cc_mix = cc_mix - G(jb)*mixing_store%drho_buffer(jb, ispin)%cc(ig)
END DO
f_mix = alpha*mixing_store%kerker_factor(ig)
rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
f_mix*res_rho(ig) + cc_mix
IF (res_norm > 1.E-14_dp) rho_g(ispin)%pw%cc(ig) = mixing_store%rhoin(ispin)%cc(ig) + &
f_mix*res_rho(ig) + cc_mix
mixing_store%rhoin_old(ispin)%cc(ig) = mixing_store%rhoin(ispin)%cc(ig)
mixing_store%rhoin(ispin)%cc(ig) = rho_g(ispin)%pw%cc(ig)
END DO
@ -735,12 +826,14 @@ CONTAINS
valh = valh - G(jb)*mixing_store%dcpc_h_in(jb, iatom, ispin)%r_coef(j, i)
vals = vals - G(jb)*mixing_store%dcpc_s_in(jb, iatom, ispin)%r_coef(j, i)
END DO
rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + valh
rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + vals
IF (res_norm > 1.E-14_dp) THEN
rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) = &
alpha*rho_atom(iatom)%cpc_h(ispin)%r_coef(j, i) + &
mixing_store%cpc_h_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + valh
rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) = &
alpha*rho_atom(iatom)%cpc_s(ispin)%r_coef(j, i) + &
mixing_store%cpc_s_in(iatom, ispin)%r_coef(j, i)*(1._dp - alpha) + vals
END IF
END DO
END DO
@ -757,6 +850,25 @@ CONTAINS
END IF
END DO ! ispin
IF (nspin == 2) THEN
DO ig = 1, ng
rho_g(2)%pw%cc(ig) = rho_g(1)%pw%cc(ig) - rho_g(2)%pw%cc(ig)
mixing_store%rhoin(2)%cc(ig) = mixing_store%rhoin(1)%cc(ig) - mixing_store%rhoin(2)%cc(ig)
END DO
IF (gapw .AND. mixing_store%gmix_p) THEN
DO iatom = 1, natom
IF (mixing_store%paw(iatom)) THEN
rho_atom(iatom)%cpc_h(2)%r_coef = rho_atom(iatom)%cpc_h(1)%r_coef - rho_atom(iatom)%cpc_h(2)%r_coef
rho_atom(iatom)%cpc_s(2)%r_coef = rho_atom(iatom)%cpc_s(1)%r_coef - rho_atom(iatom)%cpc_s(2)%r_coef
mixing_store%cpc_h_in(iatom, 2)%r_coef = mixing_store%cpc_h_in(iatom, 1)%r_coef - &
mixing_store%cpc_h_in(iatom, 2)%r_coef
mixing_store%cpc_s_in(iatom, 2)%r_coef = mixing_store%cpc_s_in(iatom, 1)%r_coef - &
mixing_store%cpc_s_in(iatom, 2)%r_coef
END IF
END DO
END IF
ENDIF
DEALLOCATE (res_rho)

View file

@ -215,7 +215,7 @@ CONTAINS
END IF
IF (mixing_method == pulay_mixing_nr) THEN
IF (ASSOCIATED(mixing_store%pulay_matrix)) DEALLOCATE (mixing_store%pulay_matrix)
ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer))
ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
mixing_store%pulay_matrix = 0.0_dp
ELSEIF (mixing_method == broyden_mixing_nr) THEN
IF (ASSOCIATED(mixing_store%abroy)) DEALLOCATE (mixing_store%abroy)
@ -280,7 +280,7 @@ CONTAINS
! *** allocate pulay
IF (mixing_method == pulay_mixing_nr) THEN
IF (.NOT. ASSOCIATED(mixing_store%pulay_matrix)) THEN
ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer))
ALLOCATE (mixing_store%pulay_matrix(nbuffer, nbuffer, nspins))
END IF
IF (.NOT. ASSOCIATED(mixing_store%rhoin_buffer)) THEN

View file

@ -16,8 +16,8 @@ CH3SH_xastpfh_list.inp 1 5e-7 -
# density mixing
c8_pmix_gapw_all.inp 1 4e-04 -302.63376760933824
c8_pmix_gapw_all_xashh.inp 1 1e-04 -297.77088288250752
c8_broy_gapw_all.inp 1 1e-13 -303.76231051467545
c8_broy_gapw_all_xashh.inp 1 2e-11 -296.95066303779714
c8_broy_gapw_all.inp 1 1e-13 -303.71942702489218
c8_broy_gapw_all_xashh.inp 1 2e-11 -294.57321406651721
# different scf_env
CO_xastpfh_gsot.inp 1 8e-14 -79.09997113686542
# XAS TP choice of core occ from input

View file

@ -38,6 +38,7 @@
ALPHA 0.1
BETA 1.0
NBUFFER 8
GMIX_P
&END
ADDED_MOS 20 20
&END SCF

View file

@ -38,6 +38,7 @@
ALPHA 0.2
BETA 1.0
NBUFFER 8
GMIX_P
&END
ADDED_MOS 20 20
&END SCF
@ -59,6 +60,7 @@
ALPHA 0.1
BETA 1.5
NBUFFER 8
GMIX_P
&END MIXING
&END

View file

@ -8,7 +8,7 @@ si8_pmix.inp 1 7e-14 -
si8_kerker.inp 1 1e-13 -31.15224294859119
si8_pulay.inp 1 2e-13 -31.16058546019718
si8_broy.inp 1 7e-14 -31.14863007620936
si8_pulay_md.inp 1 3e-14 -31.16148751713639
si8_pulay_md.inp 1 3e-12 -31.16216174849309
si8_pulay_skip.inp 1 5e-14 -31.16145405092692
# cholesky methods
si8_pulay_reduce.inp 1 2e-13 -31.16058546019707

View file

@ -1,7 +1,7 @@
si8_lsd_broy_wc_rst.inp 23 2e-09 2111.2137083255
si8_lsd_broy_wc_list.inp 23 4e-10 992.1693045962
si8_lsd_broy_wc_list_rst.inp 23 4e-10 991.9475824044
si8_lsd_broy_fm0.2.inp 48 1.0E-14 -0.049499
si8_lsd_broy_fm0.2.inp 48 1.0E-14 -0.050162
#
c8_kerker.inp 1 1e-13 -44.11857062837854
c8_pmix.inp 1 2e-13 -43.56688833744322