diff --git a/src/apps/molresponse/ExcitedResponse.cpp b/src/apps/molresponse/ExcitedResponse.cpp index c3e1f466e87..168c5c88043 100644 --- a/src/apps/molresponse/ExcitedResponse.cpp +++ b/src/apps/molresponse/ExcitedResponse.cpp @@ -2161,7 +2161,7 @@ auto ExcitedResponse::update_response(World &world, X_space &Chi, XCOperator0 or first run where there should not be a problem // computed new_chi and res if (r_params.kain() && (iter > 0) && true) { - new_chi = kain_x_space_update(world, rotated_chi, new_res, kain_x_space); + new_chi = kain_x_space_update(world, rotated_chi, new_res, kain_x_space, 1e-8); } if (false) { x_space_step_restriction(world, rotated_chi, new_chi, compute_y, maxrotn); } diff --git a/src/apps/molresponse/FrequencyResponse.cpp b/src/apps/molresponse/FrequencyResponse.cpp index 5ff1bd52d11..527548962e9 100644 --- a/src/apps/molresponse/FrequencyResponse.cpp +++ b/src/apps/molresponse/FrequencyResponse.cpp @@ -69,16 +69,7 @@ void FrequencyResponse::iterate(World &world) { auto bsh_x_ops = make_bsh_operators_response(world, x_shifts, omega); std::vector bsh_y_ops; bsh_y_ops = (compute_y) ? make_bsh_operators_response(world, y_shifts, -omega) : bsh_x_ops; - auto max_rotation = .5; - if (thresh >= 1e-2) { - max_rotation = 2; - } else if (thresh >= 1e-4) { - max_rotation = 2 * x_residual_target; - } else if (thresh >= 1e-6) { - max_rotation = 2 * x_residual_target; - } else if (thresh >= 1e-7) { - max_rotation = .01; - } + auto max_rotation = .25 * x_residual_target + x_residual_target; PQ = generator(world, *this); PQ.truncate(); @@ -280,7 +271,7 @@ auto FrequencyResponse::update_response(World &world, X_space &chi, XCOperator= 0) {// & (iteration % 3 == 0)) { - new_chi = kain_x_space_update(world, chi, new_res, kain_x_space); + new_chi = kain_x_space_update(world, chi, new_res, kain_x_space, max_rotation); } inner_to_json(world, "x_update", response_context.inner(new_chi, new_chi), iter_function_data); diff --git a/src/apps/molresponse/ResponseBase.cpp b/src/apps/molresponse/ResponseBase.cpp index 625a3406c6a..106adf06813 100644 --- a/src/apps/molresponse/ResponseBase.cpp +++ b/src/apps/molresponse/ResponseBase.cpp @@ -1174,24 +1174,42 @@ auto ResponseBase::update_residual(World &world, const X_space &chi, const X_spa } auto ResponseBase::kain_x_space_update(World &world, const X_space &chi, const X_space &residual_chi, - response_solver &kain_x_space) -> X_space { + response_solver &kain_x_space, double max_rotation) -> X_space { if (r_params.print_level() >= 1) { molresponse::start_timer(world); } size_t m = chi.num_states(); size_t n = chi.num_orbitals(); X_space kain_update = chi.copy(); response_matrix update(m); + // compute the norm of the residuals + + + Tensor residual_norms(m); + bool compute_y = r_params.omega() != 0.0; if (compute_y) { auto x_vectors = to_response_matrix(chi); auto x_residuals = to_response_matrix(residual_chi); + // compute the norm of the residuals + for (const auto &b: chi.active) { residual_norms(b) = norm2(world, x_residuals[b]); }// / norm2(world, gx[b]); } + for (const auto &i: Chi.active) { - auto temp = kain_x_space[i].update(x_vectors[i], x_residuals[i]); - std::copy(temp.begin(), temp.begin() + n, kain_update.x[i].begin()); - std::copy(temp.begin() + n, temp.end(), kain_update.y[i].begin()); + if (residual_norms(i) > max_rotation) { + auto temp = kain_x_space[i].update(x_vectors[i], x_residuals[i]); + std::copy(temp.begin(), temp.begin() + n, kain_update.x[i].begin()); + std::copy(temp.begin() + n, temp.end(), kain_update.y[i].begin()); + } }; } else { - for (const auto &i: Chi.active) { kain_update.x[i] = kain_x_space[i].update(chi.x[i], residual_chi.x[i]); } + // first compute the residuals + for (const auto &b: chi.active) { + residual_norms(b) = norm2(world, residual_chi.x[b]); + }// / norm2(world, g_chi.x[b]); } + for (const auto &i: Chi.active) { + if (residual_norms(i) > max_rotation) { + kain_update.x[i] = kain_x_space[i].update(chi.x[i], residual_chi.x[i]); + } + } } if (r_params.print_level() >= 1) { molresponse::end_timer(world, "kain_x_update", "kain_x_update", iter_timing); } return kain_update; diff --git a/src/apps/molresponse/ResponseBase.hpp b/src/apps/molresponse/ResponseBase.hpp index 6031a24ca60..1cd2785320e 100644 --- a/src/apps/molresponse/ResponseBase.hpp +++ b/src/apps/molresponse/ResponseBase.hpp @@ -679,7 +679,7 @@ class ResponseBase { auto kain_x_space_update(World &world, const X_space &chi, const X_space &residual_chi, - response_solver &kain_x_space) -> X_space; + response_solver &kain_x_space, double max_rotation) -> X_space; void x_space_step_restriction(World &world, const X_space &old_Chi, X_space &temp, bool restrict_y, const double &max_bsh_rotation);