From 54e1aba2a8209d0e93b9a5176e6ebf33cbe925fe Mon Sep 17 00:00:00 2001 From: Gabriel Gerez Date: Mon, 21 Nov 2022 19:44:37 +0100 Subject: [PATCH] Add two state solver, seg faults --- src/chemistry/Molecule.cpp | 6 +- src/chemistry/Molecule.h | 8 +- src/driver.cpp | 62 +++-- src/qmfunctions/qmfunction_fwd.h | 2 + src/scf_solver/ExcitedStatesSolver.cpp | 321 +++++++++++------------- src/scf_solver/ExcitedStatesSolver.h | 10 +- src/scf_solver/LinearResponseSolver.cpp | 4 +- 7 files changed, 208 insertions(+), 205 deletions(-) diff --git a/src/chemistry/Molecule.cpp b/src/chemistry/Molecule.cpp index 4a50b63b4..7976f31ff 100644 --- a/src/chemistry/Molecule.cpp +++ b/src/chemistry/Molecule.cpp @@ -75,10 +75,10 @@ Molecule::Molecule(const std::vector &coord_str, int c, int m) void Molecule::initPerturbedOrbitals(bool dynamic) { if (dynamic) { - this->orbitals_x = std::make_shared(); - this->orbitals_y = std::make_shared(); + this->orbitals_x = std::make_shared(); + this->orbitals_y = std::make_shared(); } else { - this->orbitals_x = std::make_shared(); + this->orbitals_x = std::make_shared(); this->orbitals_y = this->orbitals_x; } } diff --git a/src/chemistry/Molecule.h b/src/chemistry/Molecule.h index 44f1ca950..c2a897ce8 100644 --- a/src/chemistry/Molecule.h +++ b/src/chemistry/Molecule.h @@ -98,8 +98,8 @@ class Molecule final { const auto &getFockMatrix() const { return this->fock_matrix; } auto getOrbitals_p() const { return this->orbitals_0; } - auto getOrbitalsX_p() const { return this->orbitals_x; } - auto getOrbitalsY_p() const { return this->orbitals_y; } + auto getOrbitalsX_p(int state) const { return ((*(this->orbitals_x))[state]); } + auto getOrbitalsY_p(int state) const { return ((*(this->orbitals_y))[state]); } auto getCavity_p() const { return this->cavity; } nlohmann::json json() const; @@ -134,8 +134,8 @@ class Molecule final { std::shared_ptr cavity{nullptr}; std::shared_ptr orbitals_0{std::make_shared()}; - std::shared_ptr orbitals_x{nullptr}; - std::shared_ptr orbitals_y{nullptr}; + std::shared_ptr orbitals_x{nullptr}; + std::shared_ptr orbitals_y{nullptr}; // Properties SCFEnergy energy{}; diff --git a/src/driver.cpp b/src/driver.cpp index 5bad5c6b8..d97dbc3c8 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -105,6 +105,7 @@ DerivativeOperator_p get_derivative(const std::string &name); template RankOneOperator get_operator(const std::string &name, const json &json_oper); template RankTwoOperator get_operator(const std::string &name, const json &json_oper); void build_fock_operator(const json &input, Molecule &mol, FockBuilder &F, int order); +void build_fock_operator(const json &input, Molecule &mol, FockBuilder &F, int order, int states); void init_properties(const json &json_prop, Molecule &mol); namespace scf { @@ -117,8 +118,11 @@ void plot_quantities(const json &input, Molecule &mol); namespace rsp { bool guess_orbitals(const json &input, Molecule &mol); +bool guess_orbitals(const json &json_guess, Molecule &mol, int state); void write_orbitals(const json &input, Molecule &mol, bool dynamic); +void write_orbitals(const json &input, Molecule &mol, bool dynamic, int state); void calc_properties(const json &input, Molecule &mol, int dir, double omega); +void calc_properties(const json &input, Molecule &mol, int dir, double omega, int state); } // namespace rsp } // namespace driver @@ -745,12 +749,18 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { auto dynamic = false; // json_rsp["dynamic"]; mol.initPerturbedOrbitals(dynamic); - FockBuilder F_1; - const auto &json_fock_1 = json_rsp["fock_operator"]; - driver::build_fock_operator(json_fock_1, mol, F_1, 1); + auto N_states = 1; + + std::vector F_1_vec; // make into a matrix for the eigenvalue problem + const auto &json_fock = json_rsp["fock_operator"]; + for (auto i = 0; i < N_states; i++) { + FockBuilder F_1; + driver::build_fock_operator(json_fock, mol, F_1, 1, i); + F_1_vec.push_back(F_1); + } const auto &json_pert = json_rsp["perturbation"]; - auto h_1 = driver::get_operator<3>(json_pert["operator"], json_pert); + // auto h_1 = driver::get_operator<3>(json_pert["operator"], json_pert); json_out["perturbation"] = json_pert["operator"]; // json_out["frequency"] = omega; json_out["components"] = {}; @@ -764,7 +774,7 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { /////////////////////////////////////////////////////////// const auto &json_guess = json_comp["initial_guess"]; - json_out["success"] = rsp::guess_orbitals(json_guess, mol); + for (auto i = 0; i < N_states; i++) { json_out["success"] = rsp::guess_orbitals(json_guess, mol, i); } /////////////////////////////////////////////////////////// ///////////// Optimizing Perturbed Orbitals //////////// @@ -785,7 +795,7 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { auto helmholtz_prec = json_comp["rsp_solver"]["helmholtz_prec"]; // LinearResponseSolver solver(dynamic); - ExcitedStatesSolver solver(dynamic); + ExcitedStatesSolver solver(dynamic, N_states); solver.setHistory(kain); solver.setMethodName(method); solver.setMaxIterations(max_iter); @@ -796,7 +806,7 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { solver.setThreshold(orbital_thrs, property_thrs); solver.setOrthPrec(orth_prec); - comp_out["rsp_solver"] = solver.optimize(mol, F_0, F_1); + comp_out["rsp_solver"] = solver.optimize(mol, F_0, F_1_vec); json_out["success"] = comp_out["rsp_solver"]["converged"]; /////////////////////////////////////////////////////////// @@ -814,8 +824,10 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { //} F_0.clear(); mpi::barrier(mpi::comm_orb); - mol.getOrbitalsX_p().reset(); // Release shared_ptr - mol.getOrbitalsY_p().reset(); // Release shared_ptr + for (auto i = 0; i < N_states; i++) { + mol.getOrbitalsX_p(0).reset(); // Release shared_ptr + mol.getOrbitalsY_p(0).reset(); // Release shared_ptr + } return json_out; } @@ -831,6 +843,10 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { * This function expects the "initial_guess" subsection of the input. */ bool driver::rsp::guess_orbitals(const json &json_guess, Molecule &mol) { + return driver::rsp::guess_orbitals(json_guess, mol, 0); +} + +bool driver::rsp::guess_orbitals(const json &json_guess, Molecule &mol, int state) { auto type = json_guess["type"]; auto prec = json_guess["prec"]; auto mw_xp = json_guess["file_x_p"]; @@ -849,8 +865,8 @@ bool driver::rsp::guess_orbitals(const json &json_guess, Molecule &mol) { auto cube_yb = json_guess["file_CUBE_y_b"]; auto &Phi = mol.getOrbitals(); - auto &X = mol.getOrbitalsX(); - auto &Y = mol.getOrbitalsY(); + auto &X = *(mol.getOrbitalsX())[state]; + auto &Y = *(mol.getOrbitalsY())[state]; auto success_x = false; X = orbital::param_copy(Phi); @@ -894,12 +910,16 @@ bool driver::rsp::guess_orbitals(const json &json_guess, Molecule &mol) { } void driver::rsp::write_orbitals(const json &json_orbs, Molecule &mol, bool dynamic) { - auto &X = mol.getOrbitalsX(); + driver::rsp::write_orbitals(json_orbs, mol, dynamic, 0); +} + +void driver::rsp::write_orbitals(const json &json_orbs, Molecule &mol, bool dynamic, int state) { + auto &X = *((mol.getOrbitalsX())[state]); orbital::save_orbitals(X, json_orbs["file_x_p"], SPIN::Paired); orbital::save_orbitals(X, json_orbs["file_x_a"], SPIN::Alpha); orbital::save_orbitals(X, json_orbs["file_x_b"], SPIN::Beta); if (dynamic) { - auto &Y = mol.getOrbitalsY(); + auto &Y = *((mol.getOrbitalsY())[state]); orbital::save_orbitals(Y, json_orbs["file_y_p"], SPIN::Paired); orbital::save_orbitals(Y, json_orbs["file_y_a"], SPIN::Alpha); orbital::save_orbitals(Y, json_orbs["file_y_b"], SPIN::Beta); @@ -912,13 +932,17 @@ void driver::rsp::write_orbitals(const json &json_orbs, Molecule &mol, bool dyna * input section, and will compute all properties which are present in this input. */ void driver::rsp::calc_properties(const json &json_prop, Molecule &mol, int dir, double omega) { + driver::rsp::calc_properties(json_prop, mol, dir, omega, 0); +} + +void driver::rsp::calc_properties(const json &json_prop, Molecule &mol, int dir, double omega, int state) { Timer t_tot, t_lap; auto plevel = Printer::getPrintLevel(); if (plevel == 1) mrcpp::print::header(1, "Computing linear response properties"); auto &Phi = mol.getOrbitals(); - auto &X = mol.getOrbitalsX(); - auto &Y = mol.getOrbitalsY(); + auto &X = *((mol.getOrbitalsX())[state]); + auto &Y = *((mol.getOrbitalsY())[state]); if (json_prop.contains("polarizability")) { t_lap.start(); @@ -986,10 +1010,14 @@ void driver::rsp::calc_properties(const json &json_prop, Molecule &mol, int dir, * perturbation order of the operators. */ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuilder &F, int order) { + driver::build_fock_operator(json_fock, mol, F, order, 0); +} + +void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuilder &F, int order, int state) { auto &nuclei = mol.getNuclei(); auto Phi_p = mol.getOrbitals_p(); - auto X_p = mol.getOrbitalsX_p(); - auto Y_p = mol.getOrbitalsY_p(); + auto X_p = mol.getOrbitalsX_p(state); + auto Y_p = mol.getOrbitalsY_p(state); /////////////////////////////////////////////////////////// /////////////// Momentum Operator ///////////////// diff --git a/src/qmfunctions/qmfunction_fwd.h b/src/qmfunctions/qmfunction_fwd.h index 49df4b92e..63dd4c7b7 100644 --- a/src/qmfunctions/qmfunction_fwd.h +++ b/src/qmfunctions/qmfunction_fwd.h @@ -23,6 +23,7 @@ * */ +#include #include #pragma once @@ -53,6 +54,7 @@ using QMFunctionVector = std::vector; class Orbital; using OrbitalChunk = std::vector>; using OrbitalVector = std::vector; +using NStatesVector = std::vector>; class Density; diff --git a/src/scf_solver/ExcitedStatesSolver.cpp b/src/scf_solver/ExcitedStatesSolver.cpp index ccb243c4f..aaa148b86 100644 --- a/src/scf_solver/ExcitedStatesSolver.cpp +++ b/src/scf_solver/ExcitedStatesSolver.cpp @@ -67,34 +67,42 @@ namespace mrchem { * Only one state is optimized for now. random guess might help to get more states out * */ -json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder &F_1) { +json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, std::vector &F_1_vec) { Timer t_tot; json json_out; + double err_o = 1.0; + double err_t = 1.0; // Setup KAIN accelerators KAIN kain_x(this->history); KAIN kain_y(this->history); OrbitalVector &Phi_0 = mol.getOrbitals(); - OrbitalVector &X_n = mol.getOrbitalsX(); - OrbitalVector &Y_n = mol.getOrbitalsY(); + + NStatesVector &X_n_vec = mol.getOrbitalsX(); + ComplexMatrix &F_mat_0 = mol.getFockMatrix(); RankZeroOperator V_0 = F_0.potential(); - RankZeroOperator V_1 = F_1.potential(); - - double err_o = 1.0; - double err_t = 1.0; + std::vector V_1_vec; + double orb_prec = adjustPrecision(err_o); + for (auto i = 0; i < n_states; i++) { // so much repeated code + RankZeroOperator V_1 = F_1_vec[i].potential(); + V_1_vec.push_back(V_1); + V_1.setup(orb_prec); + } this->error.push_back(err_t); - double orb_prec = adjustPrecision(err_o); + F_0.setup(orb_prec); - F_1.setup(orb_prec); DoubleVector errors_x = DoubleVector::Zero(Phi_0.size()); DoubleVector errors_y = DoubleVector::Zero(Phi_0.size()); - orbital::orthogonalize(this->orth_prec, X_n, Phi_0); - orbital::orthogonalize(this->orth_prec, X_n); + // orthogonalize all orbitals wrt. the ground state and each other. + for (auto i = 0; i < n_states; i++) { + orbital::orthogonalize(this->orth_prec, *(X_n_vec[i]), Phi_0); + orbital::orthogonalize(this->orth_prec, *(X_n_vec[i])); + } // orbital::orthogonalize(this->orth_prec, Y_n, Phi_0); // orbital::orthogonalize(this->orth_prec, Y_n); @@ -104,11 +112,20 @@ json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder json_out["cycles"] = {}; // compute initial omega - auto omega_n = computeOmega(Phi_0, X_n, F_0, V_1); - auto domega_n = omega_n; - this->property.push_back(omega_n); - printParameters(omega_n, F_1.perturbation().name()); // have to change a bit this here + std::vector omega_n; + std::vector domega_n; + omega_n.push_back(0.46844); // trial energies for now + omega_n.push_back(0.47844); + + domega_n.push_back(0.46844); + domega_n.push_back(0.47844); + + // auto omega_n = computeOmega(Phi_0, X_n, F_0, V_1, F_mat_0); + + this->property.push_back(omega_n[0]); + + printParameters(omega_n[0], F_1_vec[0].perturbation().name()); // have to change a bit this here if (plevel < 1) { printConvergenceHeader("State 1 energy"); if (plevel < 1) printConvergenceRow(0); @@ -123,167 +140,121 @@ json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder // Initialize SCF cycle Timer t_scf, t_lap; double orb_prec = adjustPrecision(err_o); - - ComplexMatrix F_mat_x = F_mat_0 + omega_n * ComplexMatrix::Identity(Phi_0.size(), Phi_0.size()); - ComplexMatrix F_mat_y = F_mat_0 - omega_n * ComplexMatrix::Identity(Phi_0.size(), Phi_0.size()); - - // Setup Helmholtz operators (fixed, based on unperturbed system) double helm_prec = getHelmholtzPrec(); - HelmholtzVector H_x(helm_prec, F_mat_x.real().diagonal()); - HelmholtzVector H_y(helm_prec, F_mat_y.real().diagonal()); - - auto dot_of_X = orbital::dot(X_n, X_n).sum(); - std::cout << __FILE__ << " " << __LINE__ << ": som of the dot product : " << dot_of_X << "\n"; - - if (dynamic and plevel == 1) mrcpp::print::separator(1, '-'); - - { // Iterate X orbitals - // Compute argument: psi_i = sum_i F_0*x_j + (1 - rho_0)V_1(phi_i) - Timer t_arg; - mrcpp::print::header(2, "Computing Helmholtz argument"); - t_lap.start(); - V_1.setup(orb_prec); - OrbitalVector Psi = V_1(Phi_0); - mrcpp::print::time(2, "Applying V_1", t_lap); - t_lap.start(); - orbital::orthogonalize(this->orth_prec, Psi, Phi_0); - mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); - - mrcpp::print::footer(2, t_arg, 2); - if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_arg); - - // Apply Helmholtz operators - OrbitalVector X_np1 = H_x.apply(V_0, X_n, Psi); - Psi.clear(); - // Projecting (1 - rho_0)X - mrcpp::print::header(2, "Projecting occupied space"); - t_lap.start(); - orbital::orthogonalize(this->orth_prec, X_np1, Phi_0); - orbital::orthogonalize(this->orth_prec, X_np1); - - mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); - mrcpp::print::footer(2, t_lap, 2); - - if (plevel == 1) mrcpp::print::time(1, "Projecting occupied space", t_lap); - // orbital::orthogonalize(this->orth_prec, X_np1); // X_n orbitals should be orthogonal wrt. each other - // before next iteration, maybe do this before kain acceleration - // Compute update and errors - OrbitalVector dX_n = orbital::add(1.0, X_np1, -1.0, X_n); - domega_n = updateOmega(X_n, X_np1); - errors_x = orbital::get_norms(dX_n); - - // Compute KAIN update: - kain_x.accelerate(orb_prec, X_n, dX_n); - - // Prepare for next iteration - X_n = orbital::add(1.0, X_n, 1.0, dX_n); - // orbital::orthogonalize(this->orth_prec, X_n); - - // Setup perturbed Fock operator (including V_1) - F_1.setup(orb_prec); // do the x orbitals being uodated change this? it obviously should, but check - - // Compute omega - mrcpp::print::header(2, "Computing frequency update"); - t_lap.start(); - auto omega_np1 = computeOmega(Phi_0, X_n, F_0, V_1); /* computeOmega(Phi_0, X_n, Y_n, F_0, F_1, F_mat_0); */ - domega_n = omega_np1 - omega_n; // maybe I should do this before normalization - omega_n += domega_n; - mrcpp::print::footer(2, t_lap, 2); - if (plevel == 1) mrcpp::print::time(1, "Computing frequency update", t_lap); - this->property.push_back(omega_n); - X_np1.clear(); - - // Save checkpoint file - if (this->checkpoint) orbital::save_orbitals(X_n, this->chkFileX); - } - - if (dynamic and plevel == 1) mrcpp::print::separator(1, '-'); - - if (dynamic) { // Iterate Y orbitals - // Compute argument: psi_i = F_0*y_i + (1 - rho_0)V_1.dagger(phi_i) - Timer t_arg; - mrcpp::print::header(2, "Computing Helmholtz argument"); - t_lap.start(); - OrbitalVector Psi = V_1.dagger(Phi_0); - mrcpp::print::time(2, "Applying V_1.dagger()", t_lap); - - t_lap.start(); - orbital::orthogonalize(this->orth_prec, Psi, Phi_0); - mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); - - mrcpp::print::footer(2, t_arg, 2); - if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_arg); - - t_lap.start(); - orbital::orthogonalize(this->orth_prec, Y_n, Phi_0); - mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); - - // Apply Helmholtz operators - OrbitalVector Y_np1 = H_y.apply(V_0, Y_n, Psi); // quite sure this might be wrong - Psi.clear(); - - // Projecting (1 - rho_0)X - mrcpp::print::header(2, "Projecting occupied space"); - t_lap.start(); - orbital::orthogonalize(this->orth_prec, Y_np1, Phi_0); - mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); - mrcpp::print::footer(2, t_lap, 2); - if (plevel == 1) mrcpp::print::time(1, "Projecting occupied space", t_lap); - - orbital::orthogonalize(this->orth_prec, Y_np1); // X_n orbitals should be orthogonal wrt. each other - // before next iteration, maybe do this before kain acceler - // Compute update and errors - OrbitalVector dY_n = orbital::add(1.0, Y_np1, -1.0, Y_n); - errors_y = orbital::get_norms(dY_n); - Y_np1.clear(); - - // Compute KAIN update: - kain_y.accelerate(orb_prec, Y_n, dY_n); - - // Prepare for next iteration - Y_n = orbital::add(1.0, Y_n, 1.0, dY_n); - - // Save checkpoint file - if (this->checkpoint) orbital::save_orbitals(Y_n, this->chkFileY); - } - // Compute errors - err_o = std::max(errors_x.maxCoeff(), errors_y.maxCoeff()); - err_t = std::sqrt(errors_x.dot(errors_x) + errors_y.dot(errors_y)); - json_cycle["mo_residual"] = err_t; - // Collect convergence data - this->error.push_back(err_t); - double err_w = domega_n; - converged = checkConvergence(err_o, err_w); - json_cycle["frequency"] = omega_n; - json_cycle["frequency_update"] = err_w; - - // Finalize SCF cycle - if (plevel < 1) printConvergenceRow(nIter + 1); - printOrbitals(orbital::get_norms(X_n), errors_x, X_n, 1); - if (dynamic) printOrbitals(orbital::get_norms(Y_n), errors_y, Y_n, 1, false); - mrcpp::print::separator(1, '-'); - printResidual(err_t, converged); - mrcpp::print::separator(2, '=', 2); - printProperty(); - printMemory(); - t_scf.stop(); - json_cycle["wall_time"] = t_scf.elapsed(); - mrcpp::print::footer(1, t_scf, 2, '#'); - mrcpp::print::separator(2, ' ', 2); - json_out["cycles"].push_back(json_cycle); - if (converged) { - mrcpp::print::header(2, "Computing frequency"); - t_lap.start(); - omega_n = computeOmega(Phi_0, X_n, F_0, V_1); /* computeOmega(Phi_0, X_n, Y_n, F_0, F_1, F_mat_0); */ - mrcpp::print::footer(2, t_lap, 2); - if (plevel == 1) mrcpp::print::time(1, "Computing frequency", t_lap); - this->property.push_back(omega_n); - F_1.clear(); - break; + std::vector H_x_vec; + for (auto i = 0; i < n_states; i++) { + ComplexMatrix F_mat_x = F_mat_0 + omega_n[0] * ComplexMatrix::Identity(Phi_0.size(), Phi_0.size()); + + // Setup Helmholtz operators (fixed, based on unperturbed system) + HelmholtzVector H_x(helm_prec, F_mat_x.real().diagonal()); + H_x_vec.push_back(H_x); + + auto dot_of_X = orbital::dot(*(X_n_vec[i]), *(X_n_vec[i])).sum(); + std::cout << __FILE__ << " " << __LINE__ << ": som of the dot product : " << dot_of_X << "\n"; + + OrbitalVector &X_n = *(X_n_vec[i]); + RankZeroOperator V_1 = F_1_vec[i].potential(); + + if (dynamic and plevel == 1) mrcpp::print::separator(1, '-'); + + { // Iterate X orbitals + // Compute argument: psi_i = sum_i F_0*x_j + (1 - rho_0)V_1(phi_i) + Timer t_arg; + mrcpp::print::header(2, "Computing Helmholtz argument"); + t_lap.start(); + V_1.setup(orb_prec); + OrbitalVector Psi = V_1(Phi_0); + mrcpp::print::time(2, "Applying V_1", t_lap); + t_lap.start(); + orbital::orthogonalize(this->orth_prec, Psi, Phi_0); + mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); + + mrcpp::print::footer(2, t_arg, 2); + if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_arg); + + // Apply Helmholtz operators + OrbitalVector X_np1 = H_x.apply(V_0, X_n, Psi); + Psi.clear(); + // Projecting (1 - rho_0)X + mrcpp::print::header(2, "Projecting occupied space"); + t_lap.start(); + orbital::orthogonalize(this->orth_prec, X_np1, Phi_0); + orbital::orthogonalize(this->orth_prec, X_np1); + + mrcpp::print::time(2, "Projecting (1 - rho_0)", t_lap); + mrcpp::print::footer(2, t_lap, 2); + + if (plevel == 1) mrcpp::print::time(1, "Projecting occupied space", t_lap); + + // Compute update and errors + OrbitalVector dX_n = orbital::add(1.0, X_np1, -1.0, X_n); + domega_n[i] = updateOmega(X_n, X_np1); + errors_x = orbital::get_norms(dX_n); + + // Compute KAIN update: + kain_x.accelerate(orb_prec, X_n, dX_n); + + // Prepare for next iteration + X_n = orbital::add(1.0, X_n, 1.0, dX_n); + // orbital::orthogonalize(this->orth_prec, X_n); + + // Setup perturbed Fock operator (including V_1) + V_1.clear(); + V_1.setup(orb_prec); // do the x orbitals being uodated change this? it obviously should, but check + + // Compute omega + mrcpp::print::header(2, "Computing frequency update"); + t_lap.start(); + auto omega_np1 = computeOmega(Phi_0, X_n, F_0, V_1, F_mat_0); /* computeOmega(Phi_0, X_n, Y_n, F_0, F_1, F_mat_0); */ + domega_n[i] = omega_np1 - omega_n[i]; // maybe I should do this before normalization + omega_n[i] += domega_n[i]; + mrcpp::print::footer(2, t_lap, 2); + if (plevel == 1) mrcpp::print::time(1, "Computing frequency update", t_lap); + this->property.push_back(omega_n[0]); + X_np1.clear(); + + // Save checkpoint file + if (this->checkpoint) orbital::save_orbitals(X_n, this->chkFileX); + } + + // Compute errors + err_o = std::max(errors_x.maxCoeff(), errors_y.maxCoeff()); + err_t = std::sqrt(errors_x.dot(errors_x) + errors_y.dot(errors_y)); + json_cycle["mo_residual"] = err_t; + // Collect convergence data + this->error.push_back(err_t); + double err_w = domega_n[0]; + converged = checkConvergence(err_o, err_w); + json_cycle["frequency"] = omega_n; + json_cycle["frequency_update"] = err_w; + + // Finalize SCF cycle + if (plevel < 1) printConvergenceRow(nIter + 1); + printOrbitals(orbital::get_norms(X_n), errors_x, X_n, 1); + mrcpp::print::separator(1, '-'); + printResidual(err_t, converged); + mrcpp::print::separator(2, '=', 2); + printProperty(); + printMemory(); + t_scf.stop(); + json_cycle["wall_time"] = t_scf.elapsed(); + mrcpp::print::footer(1, t_scf, 2, '#'); + mrcpp::print::separator(2, ' ', 2); + json_out["cycles"].push_back(json_cycle); + if (converged) { + mrcpp::print::header(2, "Computing frequency"); + t_lap.start(); + omega_n[i] = computeOmega(Phi_0, X_n, F_0, V_1, F_mat_0); /* computeOmega(Phi_0, X_n, Y_n, F_0, F_1, F_mat_0); */ + mrcpp::print::footer(2, t_lap, 2); + if (plevel == 1) mrcpp::print::time(1, "Computing frequency", t_lap); + this->property.push_back(omega_n[0]); + V_1.clear(); + break; + } + // Clear perturbed Fock operator + V_1.clear(); } - // Clear perturbed Fock operator - F_1.clear(); } // Compute property @@ -296,7 +267,7 @@ json ExcitedStatesSolver::optimize(Molecule &mol, FockBuilder &F_0, FockBuilder } /** @brief consider only diagonals of the A and S matrices, first for single state. maybe only ok for tda */ -double ExcitedStatesSolver::computeOmega(OrbitalVector &Phi, OrbitalVector &X, FockBuilder &F_0, RankZeroOperator &V_1) { +double ExcitedStatesSolver::computeOmega(OrbitalVector &Phi, OrbitalVector &X, FockBuilder &F_0, RankZeroOperator &V_1, ComplexMatrix &F_mat_0) { /* std::cout << "diagonal of F_mat_0 " << F_mat_0.diagonal() << "\n"; auto &p_op = F_0.momentum(); @@ -335,7 +306,7 @@ double ExcitedStatesSolver::computeOmega(OrbitalVector &Phi, OrbitalVector &X, F MSG_INFO("\ncomputing omega\n") auto xi_t_xi_vec = orbital::dot(X, X); std::cout << "vector of : " << xi_t_xi_vec << "\n"; - auto ei_vec = F_0(Phi, Phi).diagonal(); // could probably append X to Phi_0 for this and then remove them after the computation + auto ei_vec = F_mat_0.diagonal(); // could probably append X to Phi_0 for this and then remove them after the computation std::cout << "vector of e_i = : " << ei_vec << "\n"; auto sum_ei_xi_t_xi = ei_vec.dot(xi_t_xi_vec); std::cout << "sum_i e_i : " << sum_ei_xi_t_xi << "\n"; diff --git a/src/scf_solver/ExcitedStatesSolver.h b/src/scf_solver/ExcitedStatesSolver.h index cd84e1957..106fed373 100644 --- a/src/scf_solver/ExcitedStatesSolver.h +++ b/src/scf_solver/ExcitedStatesSolver.h @@ -40,11 +40,12 @@ class FockBuilder; class ExcitedStatesSolver final : public SCFSolver { public: - explicit ExcitedStatesSolver(bool dyn = false) - : dynamic(dyn) {} + explicit ExcitedStatesSolver(bool dyn = false, int N_states = 1) + : dynamic(dyn) + , n_states(N_states) {} ~ExcitedStatesSolver() override = default; - nlohmann::json optimize(Molecule &mol, FockBuilder &F_0, FockBuilder &F_1); + nlohmann::json optimize(Molecule &mol, FockBuilder &F_0, std::vector &F_1_vec); void setOrthPrec(double prec) { this->orth_prec = prec; } void setCheckpointFile(const std::string &file_x, const std::string &file_y) { this->chkFileX = file_x; @@ -53,11 +54,12 @@ class ExcitedStatesSolver final : public SCFSolver { protected: const bool dynamic; + const int n_states; double orth_prec{mrcpp::MachineZero}; std::string chkFileX; ///< Name of checkpoint file std::string chkFileY; ///< Name of checkpoint file - double computeOmega(OrbitalVector &Phi_0, OrbitalVector &X, FockBuilder &F_0, RankZeroOperator &V_1); + double computeOmega(OrbitalVector &Phi_0, OrbitalVector &X, FockBuilder &F_0, RankZeroOperator &V_1, ComplexMatrix &F_mat_0); double updateOmega(OrbitalVector &X_n, OrbitalVector &X_np1); void printProperty() const; void printParameters(double omega, const std::string &oper) const; diff --git a/src/scf_solver/LinearResponseSolver.cpp b/src/scf_solver/LinearResponseSolver.cpp index cf92b412f..ecfa91964 100644 --- a/src/scf_solver/LinearResponseSolver.cpp +++ b/src/scf_solver/LinearResponseSolver.cpp @@ -69,8 +69,8 @@ json LinearResponseSolver::optimize(double omega, Molecule &mol, FockBuilder &F_ KAIN kain_x(this->history); KAIN kain_y(this->history); OrbitalVector &Phi_0 = mol.getOrbitals(); - OrbitalVector &X_n = mol.getOrbitalsX(); - OrbitalVector &Y_n = mol.getOrbitalsY(); + OrbitalVector &X_n = *((mol.getOrbitalsX())[0]); + OrbitalVector &Y_n = *((mol.getOrbitalsY())[0]); ComplexMatrix &F_mat_0 = mol.getFockMatrix(); ComplexMatrix F_mat_x = F_mat_0 + omega * ComplexMatrix::Identity(Phi_0.size(), Phi_0.size()); ComplexMatrix F_mat_y = F_mat_0 - omega * ComplexMatrix::Identity(Phi_0.size(), Phi_0.size());