diff --git a/src/qs_density_mixing_types.F b/src/qs_density_mixing_types.F index dffc3a9bbb..c863197379 100644 --- a/src/qs_density_mixing_types.F +++ b/src/qs_density_mixing_types.F @@ -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 diff --git a/src/qs_gspace_mixing.F b/src/qs_gspace_mixing.F index 3efeaadea6..6c72a93901 100644 --- a/src/qs_gspace_mixing.F +++ b/src/qs_gspace_mixing.F @@ -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) diff --git a/src/qs_mixing_utils.F b/src/qs_mixing_utils.F index 3dba26cad8..1f5148c9a4 100644 --- a/src/qs_mixing_utils.F +++ b/src/qs_mixing_utils.F @@ -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 diff --git a/tests/QS/regtest-gapw-1/TEST_FILES b/tests/QS/regtest-gapw-1/TEST_FILES index 7addb5755c..dade6159b7 100644 --- a/tests/QS/regtest-gapw-1/TEST_FILES +++ b/tests/QS/regtest-gapw-1/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-gapw-1/c8_broy_gapw_all.inp b/tests/QS/regtest-gapw-1/c8_broy_gapw_all.inp index 2663f9e685..6161a94fac 100644 --- a/tests/QS/regtest-gapw-1/c8_broy_gapw_all.inp +++ b/tests/QS/regtest-gapw-1/c8_broy_gapw_all.inp @@ -38,6 +38,7 @@ ALPHA 0.1 BETA 1.0 NBUFFER 8 + GMIX_P &END ADDED_MOS 20 20 &END SCF diff --git a/tests/QS/regtest-gapw-1/c8_broy_gapw_all_xashh.inp b/tests/QS/regtest-gapw-1/c8_broy_gapw_all_xashh.inp index bc567e2fff..0ed3958bb3 100644 --- a/tests/QS/regtest-gapw-1/c8_broy_gapw_all_xashh.inp +++ b/tests/QS/regtest-gapw-1/c8_broy_gapw_all_xashh.inp @@ -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 diff --git a/tests/QS/regtest-gpw-5/TEST_FILES b/tests/QS/regtest-gpw-5/TEST_FILES index 17402360bd..96fa296959 100644 --- a/tests/QS/regtest-gpw-5/TEST_FILES +++ b/tests/QS/regtest-gpw-5/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-gpw-6-3/TEST_FILES b/tests/QS/regtest-gpw-6-3/TEST_FILES index 53e92be07e..8a1f10d0ec 100644 --- a/tests/QS/regtest-gpw-6-3/TEST_FILES +++ b/tests/QS/regtest-gpw-6-3/TEST_FILES @@ -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