From ddd02c5d4c9c6606d842092b6628ec2d0dfb1543 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Thu, 21 Apr 2022 17:38:14 +0200 Subject: [PATCH 01/34] Added functionality for determining the an environment label for displaying SCF solver settings. --- python/mrchem/helpers.py | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py index 8eadda3cb..c0bc4f1f0 100644 --- a/python/mrchem/helpers.py +++ b/python/mrchem/helpers.py @@ -195,10 +195,39 @@ def write_scf_solver(user_dict, wf_dict): if start_prec < 0.0: start_prec = final_prec + # TODO: Improve label determination. Make it easier to extent for future env additions + # TODO: Move determination to its own function to clean up this one + # TODO: Control precision of x,y,z strengths in label + # Determine the appropriate `Environment` printout in SCF + # The relevant algorithms are implicit solvation and finite external fields + env = user_dict['Environment'] + ext = user_dict['ExternalFields'] + + # Keep track of active environments + has_solvent = env['run_environment'] + has_external_fields = len(ext['electric_field']) > 0 + + # If no external fields, then the list will be empty + # Need to catch the exception and store placeholders + try: + x, y, z = ext['electric_field'] + except ValueError: + x, y, z = None, None, None # Useless placeholders + + # Labels to aggregate + label_sol = 'PCM' + label_ext = f'Electric field ({x}, {y}, {z})' + + # Aggregate labels + labels = [label_sol, label_ext] + envs = [has_solvent, has_external_fields] + label_env = ' ; '.join([label for label, has_env in zip(labels, envs) if has_env]) + scf_dict = user_dict["SCF"] solver_dict = { "method": wf_dict["method_name"], "relativity": wf_dict["relativity_name"], + "environment": label_env, "kain": scf_dict["kain"], "max_iter": scf_dict["max_iter"], "rotation": scf_dict["rotation"], @@ -211,6 +240,7 @@ def write_scf_solver(user_dict, wf_dict): "orbital_thrs": scf_dict["orbital_thrs"], "helmholtz_prec": user_dict["Precisions"]["helmholtz_prec"] } + return solver_dict @@ -446,5 +476,3 @@ def parse_wf_method(user_dict): "dft_funcs": dft_funcs } return wf_dict - - From c5092f07eee587a28ffaac9414eb02641fb060ce Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 21 Apr 2022 16:22:26 +0200 Subject: [PATCH 02/34] Pass unique_ptr by value to ReactionOperator This is to emphasize transfer of ownership by explicit std::move in the constructor call. --- src/driver.cpp | 2 +- src/qmoperators/two_electron/ReactionOperator.h | 4 ++-- src/qmoperators/two_electron/ReactionPotential.cpp | 2 +- src/qmoperators/two_electron/ReactionPotential.h | 2 +- tests/solventeffect/reaction_operator.cpp | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/driver.cpp b/src/driver.cpp index d766b677b..f5be254a3 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -1034,7 +1034,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild Permittivity dielectric_func(*cavity_r, eps_in_r, eps_out_r, formulation); auto scrf_p = std::make_unique(dielectric_func, nuclei, P_r, D_r, poisson_prec, hist_r, max_iter, accelerate_Vr, convergence_criterion, algorithm, density_type); - auto V_R = std::make_shared(scrf_p, Phi_p); + auto V_R = std::make_shared(std::move(scrf_p), Phi_p); F.getReactionOperator() = V_R; } /////////////////////////////////////////////////////////// diff --git a/src/qmoperators/two_electron/ReactionOperator.h b/src/qmoperators/two_electron/ReactionOperator.h index 58f12a2e7..71e2eb95d 100644 --- a/src/qmoperators/two_electron/ReactionOperator.h +++ b/src/qmoperators/two_electron/ReactionOperator.h @@ -42,8 +42,8 @@ class SCRF; class ReactionOperator final : public RankZeroOperator { public: - ReactionOperator(std::unique_ptr &scrf_p, std::shared_ptr Phi_p) { - potential = std::make_shared(scrf_p, Phi_p); + ReactionOperator(std::unique_ptr scrf_p, std::shared_ptr Phi_p) { + potential = std::make_shared(std::move(scrf_p), Phi_p); // Invoke operator= to assign *this operator RankZeroOperator &V = (*this); V = potential; diff --git a/src/qmoperators/two_electron/ReactionPotential.cpp b/src/qmoperators/two_electron/ReactionPotential.cpp index da723d7b3..ea2fbedb2 100644 --- a/src/qmoperators/two_electron/ReactionPotential.cpp +++ b/src/qmoperators/two_electron/ReactionPotential.cpp @@ -32,7 +32,7 @@ using OrbitalVector_p = std::shared_ptr; namespace mrchem { -ReactionPotential::ReactionPotential(SCRF_p &scrf_p, OrbitalVector_p Phi_p) +ReactionPotential::ReactionPotential(SCRF_p scrf_p, OrbitalVector_p Phi_p) : QMPotential(1, false) , helper(std::move(scrf_p)) , Phi(Phi_p) {} diff --git a/src/qmoperators/two_electron/ReactionPotential.h b/src/qmoperators/two_electron/ReactionPotential.h index 4af12b514..7aade0014 100644 --- a/src/qmoperators/two_electron/ReactionPotential.h +++ b/src/qmoperators/two_electron/ReactionPotential.h @@ -44,7 +44,7 @@ class ReactionPotential final : public QMPotential { /** @brief Initializes the ReactionPotential class. * @param scrf_p A SCRF instance which contains the parameters needed to compute the ReactionPotential. * @param Phi_p A pointer to a vector which contains the orbitals optimized in the SCF procedure. */ - ReactionPotential(std::unique_ptr &scrf_p, std::shared_ptr Phi_p); + ReactionPotential(std::unique_ptr scrf_p, std::shared_ptr Phi_p); ~ReactionPotential() override { free(NUMBER::Total); } SCRF *getHelper() { return this->helper.get(); } diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp index 60e553d77..23dc9d234 100644 --- a/tests/solventeffect/reaction_operator.cpp +++ b/tests/solventeffect/reaction_operator.cpp @@ -90,7 +90,7 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") { int kain = 4; auto scrf_p = std::make_unique(dielectric_func, molecule, P_p, D_p, prec, kain, 100, true, "dynamic", "scrf", "total"); - auto Reo = std::make_shared(scrf_p, Phi_p); + auto Reo = std::make_shared(std::move(scrf_p), Phi_p); Reo->setTesting(); Reo->setup(prec); double total_energy = Reo->getTotalEnergy(); From 92a6b53ba214ad1b80e6459ccb2563fe8bb23c56 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Thu, 21 Apr 2022 19:41:06 +0200 Subject: [PATCH 03/34] Added option for right-aligning text printing, and for adding a spacing buffer when printing long text or scalar names --- src/utils/print_utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h index 25127a510..86e4142ed 100644 --- a/src/utils/print_utils.h +++ b/src/utils/print_utils.h @@ -37,7 +37,7 @@ void headline(int level, const std::string &txt); void text(int level, const std::string &txt, const std::string &val); void json(int level, const nlohmann::json &j, bool ralign = false); void coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p = -1, bool s = false); -void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false); +void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false, int txtBuffer = 0); void vector(int level, const std::string &txt, const DoubleVector &val, int p = -1, bool s = false); void matrix(int level, const std::string &txt, const DoubleMatrix &val, int p = -1, bool s = false); void qmfunction(int level, const std::string &txt, const QMFunction &func, mrcpp::Timer &t); From 185e58d2ce48d3c5cc348afcbcc5565873a2cb7d Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Thu, 21 Apr 2022 19:42:53 +0200 Subject: [PATCH 04/34] Print information concerning the environment. An aggregated env label is printed in SCF solver settings, and a larger table is printed of some key SCRF settings. --- src/driver.cpp | 4 ++++ src/environment/SCRF.cpp | 18 ++++++++++++++++++ src/environment/SCRF.h | 1 + src/scf_solver/GroundStateSolver.cpp | 1 + src/scf_solver/SCFSolver.h | 2 ++ 5 files changed, 26 insertions(+) diff --git a/src/driver.cpp b/src/driver.cpp index f5be254a3..77ecb9c50 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -258,6 +258,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) { auto kain = json_scf["scf_solver"]["kain"]; auto method = json_scf["scf_solver"]["method"]; auto relativity = json_scf["scf_solver"]["relativity"]; + auto environment = json_scf["scf_solver"]["environment"]; auto max_iter = json_scf["scf_solver"]["max_iter"]; auto rotation = json_scf["scf_solver"]["rotation"]; auto localize = json_scf["scf_solver"]["localize"]; @@ -275,6 +276,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) { solver.setLocalize(localize); solver.setMethodName(method); solver.setRelativityName(relativity); + solver.setEnvironmentName(environment); solver.setCheckpoint(checkpoint); solver.setCheckpointFile(file_chk); solver.setMaxIterations(max_iter); @@ -1036,6 +1038,8 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto scrf_p = std::make_unique(dielectric_func, nuclei, P_r, D_r, poisson_prec, hist_r, max_iter, accelerate_Vr, convergence_criterion, algorithm, density_type); auto V_R = std::make_shared(std::move(scrf_p), Phi_p); F.getReactionOperator() = V_R; + + V_R->getHelper()->printParameters(); } /////////////////////////////////////////////////////////// //////////////////// XC Operator ////////////////////// diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index fd1fab2d8..50e716a30 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -26,6 +26,7 @@ #include "SCRF.h" #include +#include #include "chemistry/PhysicalConstants.h" #include "chemistry/chemistry_utils.h" @@ -290,4 +291,21 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { gamma_np1.free(NUMBER::Real); } +void SCRF::printParameters() { + int buffer = 3; // Extra space for long names + mrcpp::print::header(0, "Self-Consistent Reaction Field"); + print_utils::scalar(0, "Dielectric constant (inside) ", epsilon.getEpsIn(), "", 5, true, buffer); + print_utils::scalar(0, "Dielectric constant (outside)", epsilon.getEpsOut(), "", 5, true, buffer); + print_utils::scalar(0, "Max number of micro-ierations", max_iter, "", 1, false, buffer); + print_utils::text(0, "Accelerate with KAIN", (accelerate_Vr) ? "true" : "false", true, buffer); + print_utils::text(0, "Algorithm", algorithm, true, buffer); + print_utils::text(0, "Density type", density_type, true, buffer); + print_utils::text(0, "Convergence criterion", convergence_criterion, true, buffer); + mrcpp::print::separator(0, '='); + + // A couple of new lines to have a pleasing spacing to next section + println(0, ""); + println(0, ""); +} + } // namespace mrchem diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h index 6f1d20a4d..ec807f752 100644 --- a/src/environment/SCRF.h +++ b/src/environment/SCRF.h @@ -60,6 +60,7 @@ class SCRF final { QMFunction &getCurrentDifferenceGamma() { return this->dgamma_n; } void updateMOResidual(double const err_t) { this->mo_residual = err_t; } + void printParameters(); friend class ReactionPotential; diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp index a211ddf05..e04e5d347 100644 --- a/src/scf_solver/GroundStateSolver.cpp +++ b/src/scf_solver/GroundStateSolver.cpp @@ -171,6 +171,7 @@ void GroundStateSolver::printParameters(const std::string &calculation) const { print_utils::text(0, "Calculation ", calculation); print_utils::text(0, "Method ", this->methodName); print_utils::text(0, "Relativity ", this->relativityName); + print_utils::text(0, "Environment ", this->environmentName); print_utils::text(0, "Checkpointing ", (this->checkpoint) ? "On" : "Off"); print_utils::text(0, "Max iterations ", o_iter.str()); print_utils::text(0, "KAIN solver ", o_kain.str()); diff --git a/src/scf_solver/SCFSolver.h b/src/scf_solver/SCFSolver.h index a875980fa..b6e8fe4c7 100644 --- a/src/scf_solver/SCFSolver.h +++ b/src/scf_solver/SCFSolver.h @@ -58,6 +58,7 @@ class SCFSolver { void setMaxIterations(int iter) { this->maxIter = iter; } void setMethodName(const std::string &name) { this->methodName = name; } void setRelativityName(const std::string &name) { this->relativityName = name; } + void setEnvironmentName(const std::string &name) { this->environmentName = name; } protected: int history{0}; ///< Maximum length of KAIN history @@ -69,6 +70,7 @@ class SCFSolver { double orbPrec[3]{-1.0, -1.0, -1.0}; ///< Dynamic precision: [current_prec, start_prec, end_prec] std::string methodName; ///< Name of electronic structure method to appear in output std::string relativityName{"None"}; ///< Name of ZORA method + std::string environmentName{"None"}; ///< Aggregated name for external environment std::vector error; ///< Convergence orbital error std::vector property; ///< Convergence property error From c1db194341cfb08048014b3e438ecbb96d53ccc7 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 25 Apr 2022 09:59:42 +0200 Subject: [PATCH 05/34] Fixes suggested by stig: Typo, integer print, and newlines. --- src/environment/SCRF.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 50e716a30..4a2fe4107 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -296,16 +296,12 @@ void SCRF::printParameters() { mrcpp::print::header(0, "Self-Consistent Reaction Field"); print_utils::scalar(0, "Dielectric constant (inside) ", epsilon.getEpsIn(), "", 5, true, buffer); print_utils::scalar(0, "Dielectric constant (outside)", epsilon.getEpsOut(), "", 5, true, buffer); - print_utils::scalar(0, "Max number of micro-ierations", max_iter, "", 1, false, buffer); + print_utils::scalar(0, "Max number of micro-iterations", max_iter, "", 0, false, buffer); print_utils::text(0, "Accelerate with KAIN", (accelerate_Vr) ? "true" : "false", true, buffer); print_utils::text(0, "Algorithm", algorithm, true, buffer); print_utils::text(0, "Density type", density_type, true, buffer); print_utils::text(0, "Convergence criterion", convergence_criterion, true, buffer); - mrcpp::print::separator(0, '='); - - // A couple of new lines to have a pleasing spacing to next section - println(0, ""); - println(0, ""); + mrcpp::print::separator(0, '=', 2); } } // namespace mrchem From de65cc3e7dc2758ce92a88e0a3e12c96d3069346 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 25 Apr 2022 11:41:11 +0200 Subject: [PATCH 06/34] Added a method for printing json objects that automatically adjusts the space needed for its keys, with option for right and left-aligning the values. --- src/utils/print_utils.cpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index d2bb985fb..334ad9a2f 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -147,6 +147,30 @@ void print_utils::json(int level, const nlohmann::json &j, bool ralign) { } } +void print_utils::json(int level, const nlohmann::json &j, bool ralign) { + // Determine longest name + int w = 0; + for (const auto &item : j.items()) { + if (item.key().size() > w) w = item.key().size(); + } + + // Print + for (const auto &item : j.items()) { + std::string key = item.key(); + std::stringstream o_val; + o_val << item.value(); + std::string val = o_val.str(); + + // Remove quotes from val and print + val.erase(std::remove(val.begin(), val.end(), '\"'), val.end()); + + // If right-align, determine how much to shift the vals + int shift = (ralign) ? Printer::getWidth() - w - val.size() - 3 : 0; + + std::printf("%-*s%-s%-s%-s\n", w, key.c_str(), " : ", std::string(shift, ' ').c_str(), val.c_str()); + } +} + void print_utils::coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p, bool s) { if (p < 0) p = Printer::getPrecision(); int w0 = Printer::getWidth() - 2; From 35eeb3883d06ce3293e4062c6f45dc9dd9fa415d Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 25 Apr 2022 11:41:29 +0200 Subject: [PATCH 07/34] Use the json print method for printing SCRF settings. --- src/environment/SCRF.cpp | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 4a2fe4107..7c5ba9acf 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -36,6 +36,8 @@ #include "scf_solver/KAIN.h" #include "utils/print_utils.h" +#include + using mrcpp::Printer; using mrcpp::Timer; @@ -292,15 +294,18 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { } void SCRF::printParameters() { - int buffer = 3; // Extra space for long names + nlohmann::json data = { + {"Dielectric constant (inside)", epsilon.getEpsIn()}, + {"Dielectric constant (outside)", epsilon.getEpsOut()}, + {"Max. number of micro-iterations", max_iter}, + {"Accelerate with KAIN", (accelerate_Vr) ? "true" : "false"}, + {"Algorithm", algorithm}, + {"Density type", density_type}, + {"Convergence criterion", convergence_criterion}, + }; + mrcpp::print::header(0, "Self-Consistent Reaction Field"); - print_utils::scalar(0, "Dielectric constant (inside) ", epsilon.getEpsIn(), "", 5, true, buffer); - print_utils::scalar(0, "Dielectric constant (outside)", epsilon.getEpsOut(), "", 5, true, buffer); - print_utils::scalar(0, "Max number of micro-iterations", max_iter, "", 0, false, buffer); - print_utils::text(0, "Accelerate with KAIN", (accelerate_Vr) ? "true" : "false", true, buffer); - print_utils::text(0, "Algorithm", algorithm, true, buffer); - print_utils::text(0, "Density type", density_type, true, buffer); - print_utils::text(0, "Convergence criterion", convergence_criterion, true, buffer); + print_utils::json(0, data, true); mrcpp::print::separator(0, '=', 2); } From 679b10edd5e85bf0f75f5230a390c6e85afccc64 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 25 Apr 2022 11:47:23 +0200 Subject: [PATCH 08/34] Undo modifications to printing text and scalars. --- src/utils/print_utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h index 86e4142ed..25127a510 100644 --- a/src/utils/print_utils.h +++ b/src/utils/print_utils.h @@ -37,7 +37,7 @@ void headline(int level, const std::string &txt); void text(int level, const std::string &txt, const std::string &val); void json(int level, const nlohmann::json &j, bool ralign = false); void coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p = -1, bool s = false); -void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false, int txtBuffer = 0); +void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false); void vector(int level, const std::string &txt, const DoubleVector &val, int p = -1, bool s = false); void matrix(int level, const std::string &txt, const DoubleMatrix &val, int p = -1, bool s = false); void qmfunction(int level, const std::string &txt, const QMFunction &func, mrcpp::Timer &t); From 8df34de0abe1dde5c8159ac5c584be9e3c147977 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 25 Apr 2022 12:28:23 +0200 Subject: [PATCH 09/34] Remove ralign traces in text, and add a check for too long names in json print. --- src/utils/print_utils.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 334ad9a2f..ac0ef2ebe 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -167,6 +167,9 @@ void print_utils::json(int level, const nlohmann::json &j, bool ralign) { // If right-align, determine how much to shift the vals int shift = (ralign) ? Printer::getWidth() - w - val.size() - 3 : 0; + // Avoid runtime errors due to negative shifts caused by very long names + if (shift < 0) shift = 0; + std::printf("%-*s%-s%-s%-s\n", w, key.c_str(), " : ", std::string(shift, ' ').c_str(), val.c_str()); } } From f3aac7621667c4055ea00a06d2ea506135646552 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 25 Apr 2022 14:42:11 +0200 Subject: [PATCH 10/34] Print information about the solvation cavity if SCRF is requested. --- src/chemistry/Molecule.cpp | 48 ++++++++++++++++++++++++++++++++++++++ src/chemistry/Molecule.h | 1 + src/driver.cpp | 1 + src/environment/Cavity.h | 1 + 4 files changed, 51 insertions(+) diff --git a/src/chemistry/Molecule.cpp b/src/chemistry/Molecule.cpp index d51db7374..8473c9691 100644 --- a/src/chemistry/Molecule.cpp +++ b/src/chemistry/Molecule.cpp @@ -208,6 +208,54 @@ void Molecule::printGeometry() const { mrcpp::print::separator(0, '=', 2); } +/** @brief Pretty output of solvation cavity spheres */ +void Molecule::printCavity() { + // Collect relevant quantities + Cavity &cavity = getCavity(); + std::vector> coords = cavity.getCoordinates(); + std::vector radii = cavity.getRadii(); + + // Set widths + auto w0 = Printer::getWidth() - 1; + auto w1 = 5; + auto w2 = 12; + auto w3 = 2 * w0 / 9; + auto w4 = w0 - w1 - w2 - 3 * w3; + + // Build table column headers + std::stringstream o_head; + o_head << std::setw(w1) << "N"; + o_head << std::setw(w2) << "Radius"; + o_head << std::string(w4 - 1, ' ') << ':'; + o_head << std::setw(w3) << "x"; + o_head << std::setw(w3) << "y"; + o_head << std::setw(w3) << "z"; + + // Print + mrcpp::print::header(0, "Solvation Cavity"); + print_utils::scalar(0, "Cavity width", cavity.getWidth(), "", 6); + mrcpp::print::separator(0, '-'); + println(0, o_head.str()); + mrcpp::print::separator(0, '-'); + for (int i = 0; i < coords.size(); i++) { + mrcpp::Coord<3> coord = coords[i]; + double r = radii[i]; + double x = coord[0]; + double y = coord[1]; + double z = coord[2]; + + std::stringstream o_coord; + o_coord << std::setw(w1) << i; + o_coord << std::setw(w2) << std::setprecision(6) << std::fixed << r; + o_coord << std::string(w4 - 1, ' ') << ':'; + o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << x; + o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << y; + o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << z; + println(0, o_coord.str()); + } + mrcpp::print::separator(0, '=', 2); +} + void Molecule::printEnergies(const std::string &txt) const { energy.print(txt); epsilon.print(txt); diff --git a/src/chemistry/Molecule.h b/src/chemistry/Molecule.h index 19e4727da..a778a5575 100644 --- a/src/chemistry/Molecule.h +++ b/src/chemistry/Molecule.h @@ -106,6 +106,7 @@ class Molecule final { void printGeometry() const; void printEnergies(const std::string &txt) const; void printProperties() const; + void printCavity(); void initPerturbedOrbitals(bool dynamic); void initCavity(std::vector> &coords, std::vector &R, double slope); diff --git a/src/driver.cpp b/src/driver.cpp index 77ecb9c50..1a387f9e6 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -1039,6 +1039,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto V_R = std::make_shared(std::move(scrf_p), Phi_p); F.getReactionOperator() = V_R; + mol.printCavity(); V_R->getHelper()->printParameters(); } /////////////////////////////////////////////////////////// diff --git a/src/environment/Cavity.h b/src/environment/Cavity.h index e55888ee5..165febb62 100644 --- a/src/environment/Cavity.h +++ b/src/environment/Cavity.h @@ -50,6 +50,7 @@ class Cavity final : public mrcpp::RepresentableFunction<3> { auto getGradVector() const { return this->gradvector; } std::vector> getCoordinates() const { return centers; } //!< Returns #centers. std::vector getRadii() const { return radii; } //!< Returns #radii. + double getWidth() { return width; } //!< Returns #width. friend class Permittivity; protected: From bd1620fdf36e6533adcea229d4a8055abfabe4f4 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Tue, 26 Apr 2022 11:07:13 +0200 Subject: [PATCH 11/34] Shortened some names in SCRF printout. --- src/environment/SCRF.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 7c5ba9acf..1285a0522 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -295,9 +295,9 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { void SCRF::printParameters() { nlohmann::json data = { - {"Dielectric constant (inside)", epsilon.getEpsIn()}, - {"Dielectric constant (outside)", epsilon.getEpsOut()}, - {"Max. number of micro-iterations", max_iter}, + {"Dielectric const. (inside)", epsilon.getEpsIn()}, + {"Dielectric const. (outside)", epsilon.getEpsOut()}, + {"Max. iterations", max_iter}, {"Accelerate with KAIN", (accelerate_Vr) ? "true" : "false"}, {"Algorithm", algorithm}, {"Density type", density_type}, From d9b923a4fd53b482dd5b6d015ecd1fd0decd2077 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Tue, 26 Apr 2022 11:17:48 +0200 Subject: [PATCH 12/34] Removed print statements of solvation and SCRF. --- src/driver.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/driver.cpp b/src/driver.cpp index 1a387f9e6..642aeafc2 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -1038,9 +1038,6 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto scrf_p = std::make_unique(dielectric_func, nuclei, P_r, D_r, poisson_prec, hist_r, max_iter, accelerate_Vr, convergence_criterion, algorithm, density_type); auto V_R = std::make_shared(std::move(scrf_p), Phi_p); F.getReactionOperator() = V_R; - - mol.printCavity(); - V_R->getHelper()->printParameters(); } /////////////////////////////////////////////////////////// //////////////////// XC Operator ////////////////////// From 10ce99690dc6121fd36b02d3efdd2f7652ed7132 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Tue, 26 Apr 2022 11:19:17 +0200 Subject: [PATCH 13/34] Modified json printer to have one empty space at start and stop, and to align with other print sections if possible. --- src/utils/print_utils.cpp | 36 +++++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index ac0ef2ebe..04216cd1c 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -149,28 +149,46 @@ void print_utils::json(int level, const nlohmann::json &j, bool ralign) { void print_utils::json(int level, const nlohmann::json &j, bool ralign) { // Determine longest name - int w = 0; + // This will be used for the spacing left of : + // if the name is too long to be aligned with + // the other sections + int lshift = 0; for (const auto &item : j.items()) { - if (item.key().size() > w) w = item.key().size(); + if (item.key().size() > lshift) lshift = item.key().size(); } - // Print + // Loop over json items for (const auto &item : j.items()) { + // Extract key and value from json std::string key = item.key(); std::stringstream o_val; o_val << item.value(); std::string val = o_val.str(); - // Remove quotes from val and print + // Remove quotes from val val.erase(std::remove(val.begin(), val.end(), '\"'), val.end()); - // If right-align, determine how much to shift the vals - int shift = (ralign) ? Printer::getWidth() - w - val.size() - 3 : 0; + // Standard shift to align all colons + int w0 = Printer::getWidth() - 2; // Defines the printable width + int w1 = w0 * 2 / 9; // Space dedicated to the json key + int w2 = w0 - 3 * w1; + int w3 = w2 - (val.size() + 1); + + // Some paddings for book keeping + int frontEndPadding = 2; // Empty space at beginning and end + int colonPadding = 3; // Two empty spaces around a single colon - // Avoid runtime errors due to negative shifts caused by very long names - if (shift < 0) shift = 0; + // Use standard spacing if longest name fits + if (w3 > lshift) lshift = w3 + frontEndPadding; - std::printf("%-*s%-s%-s%-s\n", w, key.c_str(), " : ", std::string(shift, ' ').c_str(), val.c_str()); + // Calculate the shift needed for right-aligning + int rshift = (ralign) ? Printer::getWidth() - lshift - val.size() - frontEndPadding - colonPadding : 0; + + // Check that rshift is not negative (caused by very long names) + if (rshift < 0) rshift = 0; + + // Print line + std::printf(" %-*s%-s%-s%-s \n", lshift, key.c_str(), " : ", std::string(rshift, ' ').c_str(), val.c_str()); } } From cab5dc9a323080b06e759d91e401a847b01437d7 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Sun, 1 May 2022 09:15:10 +0200 Subject: [PATCH 14/34] Remove duplicated json printer, uncomment unnecessary line. --- src/utils/print_utils.cpp | 47 +-------------------------------------- 1 file changed, 1 insertion(+), 46 deletions(-) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 04216cd1c..982a74102 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -127,7 +127,7 @@ void print_utils::json(int level, const nlohmann::json &j, bool ralign) { int w0 = Printer::getWidth() - 2; // Defines the printable width int w1 = w0 * 2 / 9; // Space dedicated to the json key int w2 = w0 - 3 * w1; - int w3 = w2 - (val.size() + 1); + // int w3 = w2 - (val.size() + 1); // Some paddings for book keeping int frontEndPadding = 2; // Empty space at beginning and end @@ -147,51 +147,6 @@ void print_utils::json(int level, const nlohmann::json &j, bool ralign) { } } -void print_utils::json(int level, const nlohmann::json &j, bool ralign) { - // Determine longest name - // This will be used for the spacing left of : - // if the name is too long to be aligned with - // the other sections - int lshift = 0; - for (const auto &item : j.items()) { - if (item.key().size() > lshift) lshift = item.key().size(); - } - - // Loop over json items - for (const auto &item : j.items()) { - // Extract key and value from json - std::string key = item.key(); - std::stringstream o_val; - o_val << item.value(); - std::string val = o_val.str(); - - // Remove quotes from val - val.erase(std::remove(val.begin(), val.end(), '\"'), val.end()); - - // Standard shift to align all colons - int w0 = Printer::getWidth() - 2; // Defines the printable width - int w1 = w0 * 2 / 9; // Space dedicated to the json key - int w2 = w0 - 3 * w1; - int w3 = w2 - (val.size() + 1); - - // Some paddings for book keeping - int frontEndPadding = 2; // Empty space at beginning and end - int colonPadding = 3; // Two empty spaces around a single colon - - // Use standard spacing if longest name fits - if (w3 > lshift) lshift = w3 + frontEndPadding; - - // Calculate the shift needed for right-aligning - int rshift = (ralign) ? Printer::getWidth() - lshift - val.size() - frontEndPadding - colonPadding : 0; - - // Check that rshift is not negative (caused by very long names) - if (rshift < 0) rshift = 0; - - // Print line - std::printf(" %-*s%-s%-s%-s \n", lshift, key.c_str(), " : ", std::string(rshift, ' ').c_str(), val.c_str()); - } -} - void print_utils::coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p, bool s) { if (p < 0) p = Printer::getPrecision(); int w0 = Printer::getWidth() - 2; From 928f66758bafc93734babec7556f1896e2cde124 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Mon, 2 May 2022 16:26:28 +0200 Subject: [PATCH 15/34] Print information about the cavity if solvation is requested. --- src/driver.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/driver.cpp b/src/driver.cpp index 642aeafc2..3df43c31c 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -141,6 +141,7 @@ void driver::init_molecule(const json &json_mol, Molecule &mol) { auto xyz = coord["xyz"]; nuclei.push_back(atom, xyz); } + mol.printGeometry(); if (json_mol.contains("cavity")) { auto json_cavity = json_mol["cavity"]; @@ -153,8 +154,8 @@ void driver::init_molecule(const json &json_mol, Molecule &mol) { auto width = json_cavity["width"]; mol.initCavity(coords, radii, width); + mol.printCavity(); } - mol.printGeometry(); } void driver::init_properties(const json &json_prop, Molecule &mol) { From 74ea843fed3ff12f00959e965713ce0f9c2c53e8 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Tue, 10 May 2022 08:10:31 +0200 Subject: [PATCH 16/34] Use tilde separators when printing physical constants. --- src/chemistry/PhysicalConstants.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/chemistry/PhysicalConstants.cpp b/src/chemistry/PhysicalConstants.cpp index 91ef82c7e..59923cf41 100644 --- a/src/chemistry/PhysicalConstants.cpp +++ b/src/chemistry/PhysicalConstants.cpp @@ -58,9 +58,9 @@ PhysicalConstants &PhysicalConstants::Initialize(const json &constants) { /** @brief Pretty print physical constants */ void PhysicalConstants::Print() { - mrcpp::print::header(0, "Physical Constants"); + mrcpp::print::separator(0, '~'); print_utils::json(0, constants_, true); - mrcpp::print::separator(0, '=', 2); + mrcpp::print::separator(0, '~', 2); } } // namespace mrchem From 17a4dcca0b77f377bdaaccf6c79003f6154e94f1 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Tue, 10 May 2022 10:08:11 +0200 Subject: [PATCH 17/34] Print SCRF settings when solvent is requested. --- src/environment/SCRF.cpp | 12 +++++++----- src/scf_solver/GroundStateSolver.cpp | 4 ++++ src/utils/print_utils.cpp | 6 ++++++ src/utils/print_utils.h | 1 + 4 files changed, 18 insertions(+), 5 deletions(-) diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 1285a0522..60310e204 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -295,8 +295,8 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { void SCRF::printParameters() { nlohmann::json data = { - {"Dielectric const. (inside)", epsilon.getEpsIn()}, - {"Dielectric const. (outside)", epsilon.getEpsOut()}, + {"Dielec. const. (in)", epsilon.getEpsIn()}, + {"Dielec. const. (out)", epsilon.getEpsOut()}, {"Max. iterations", max_iter}, {"Accelerate with KAIN", (accelerate_Vr) ? "true" : "false"}, {"Algorithm", algorithm}, @@ -304,9 +304,11 @@ void SCRF::printParameters() { {"Convergence criterion", convergence_criterion}, }; - mrcpp::print::header(0, "Self-Consistent Reaction Field"); - print_utils::json(0, data, true); - mrcpp::print::separator(0, '=', 2); + mrcpp::print::separator(0, '~'); + print_utils::centeredText(0, "Self-Consistent Reaction Field"); + mrcpp::print::separator(0, '~'); + print_utils::json(0, data, false); + mrcpp::print::separator(0, '~', 2); } } // namespace mrchem diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp index e04e5d347..0f87f7c1b 100644 --- a/src/scf_solver/GroundStateSolver.cpp +++ b/src/scf_solver/GroundStateSolver.cpp @@ -167,6 +167,8 @@ void GroundStateSolver::printParameters(const std::string &calculation) const { o_prec_0 << std::setprecision(5) << std::scientific << this->orbPrec[1]; o_prec_1 << std::setprecision(5) << std::scientific << this->orbPrec[2]; + mrcpp::print::separator(0, '~'); + print_utils::centeredText(0, "Self-Consistent Field"); mrcpp::print::separator(0, '~'); print_utils::text(0, "Calculation ", calculation); print_utils::text(0, "Method ", this->methodName); @@ -218,6 +220,8 @@ void GroundStateSolver::reset() { */ json GroundStateSolver::optimize(Molecule &mol, FockBuilder &F) { printParameters("Optimize ground state orbitals"); + if (F.getReactionOperator() != nullptr) F.getReactionOperator()->getHelper()->printParameters(); + Timer t_tot; json json_out; diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 982a74102..21bba14ce 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -72,6 +72,12 @@ std::string print_utils::dbl_to_str(double d, int p, bool sci) { } return o_dbl.str(); } +/** @brief Prints text centered according to current print width */ / void print_utils::centeredText(int level, const std::string &txt) { + int pwidth = Printer::getWidth(); + int txt_width = txt.size(); + int pre_spaces = (pwidth - txt_width) / 2; + println(level, std::string(pre_spaces, ' ') << txt); +} void print_utils::headline(int level, const std::string &txt) { auto pwidth = Printer::getWidth(); diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h index 25127a510..df62fa7ca 100644 --- a/src/utils/print_utils.h +++ b/src/utils/print_utils.h @@ -33,6 +33,7 @@ namespace mrchem { namespace print_utils { +void centeredText(int level, const std::string &txt); void headline(int level, const std::string &txt); void text(int level, const std::string &txt, const std::string &val); void json(int level, const nlohmann::json &j, bool ralign = false); From e36d7c07111236dc57eb05d6e80eb6afc891ded4 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Tue, 10 May 2022 10:24:30 +0200 Subject: [PATCH 18/34] Fix typo in docstring --- src/utils/print_utils.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 21bba14ce..954af12a1 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -72,7 +72,8 @@ std::string print_utils::dbl_to_str(double d, int p, bool sci) { } return o_dbl.str(); } -/** @brief Prints text centered according to current print width */ / void print_utils::centeredText(int level, const std::string &txt) { +/** @brief Prints text centered according to current print width */ +void print_utils::centeredText(int level, const std::string &txt) { int pwidth = Printer::getWidth(); int txt_width = txt.size(); int pre_spaces = (pwidth - txt_width) / 2; From 6774e91ac2015cc2a32f5a65ae79cf622fb6eb0c Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:01:20 +0200 Subject: [PATCH 19/34] Added a tilde headline to the printed table. --- src/chemistry/PhysicalConstants.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/chemistry/PhysicalConstants.cpp b/src/chemistry/PhysicalConstants.cpp index 59923cf41..afca4c3aa 100644 --- a/src/chemistry/PhysicalConstants.cpp +++ b/src/chemistry/PhysicalConstants.cpp @@ -58,6 +58,8 @@ PhysicalConstants &PhysicalConstants::Initialize(const json &constants) { /** @brief Pretty print physical constants */ void PhysicalConstants::Print() { + mrcpp::print::separator(0, '~'); + print_utils::centeredText(0, "Physical Constants"); mrcpp::print::separator(0, '~'); print_utils::json(0, constants_, true); mrcpp::print::separator(0, '~', 2); From 40f1743645bc282a2198df4801fd7a15cebef566 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:03:07 +0200 Subject: [PATCH 20/34] Add separate headlines for initial guess and and SCF, and remove printCavity statement. --- src/driver.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/driver.cpp b/src/driver.cpp index 3df43c31c..284009368 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -154,7 +154,7 @@ void driver::init_molecule(const json &json_mol, Molecule &mol) { auto width = json_cavity["width"]; mol.initCavity(coords, radii, width); - mol.printCavity(); + // mol.printCavity(); } } @@ -223,7 +223,7 @@ void driver::init_properties(const json &json_prop, Molecule &mol) { * This function expects the "scf_calculation" subsection of the input. */ json driver::scf::run(const json &json_scf, Molecule &mol) { - print_utils::headline(0, "Computing Ground State Wavefunction"); + // print_utils::headline(0, "Computing Ground State Wavefunction"); json json_out = {{"success", true}}; if (json_scf.contains("properties")) driver::init_properties(json_scf["properties"], mol); @@ -241,6 +241,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) { /////////////////////////////////////////////////////////// /////////////// Setting Up Initial Guess ////////////// /////////////////////////////////////////////////////////// + print_utils::headline(0, "Computing Initial Guess Wavefunction"); const auto &json_guess = json_scf["initial_guess"]; if (scf::guess_orbitals(json_guess, mol)) { scf::guess_energy(json_guess, mol, F); @@ -256,6 +257,8 @@ json driver::scf::run(const json &json_scf, Molecule &mol) { // Run GroundStateSolver if present in input JSON if (json_scf.contains("scf_solver")) { + print_utils::headline(0, "Computing Ground State Wavefunction"); + auto kain = json_scf["scf_solver"]["kain"]; auto method = json_scf["scf_solver"]["method"]; auto relativity = json_scf["scf_solver"]["relativity"]; @@ -1039,6 +1042,9 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto scrf_p = std::make_unique(dielectric_func, nuclei, P_r, D_r, poisson_prec, hist_r, max_iter, accelerate_Vr, convergence_criterion, algorithm, density_type); auto V_R = std::make_shared(std::move(scrf_p), Phi_p); F.getReactionOperator() = V_R; + + // dielectric_func.printParameters(); + // V_R->getHelper()->printParameters(); } /////////////////////////////////////////////////////////// //////////////////// XC Operator ////////////////////// From 27f9d636cb6926874ec059ec2ad102d57606f025 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:08:08 +0200 Subject: [PATCH 21/34] Added option for right-aligning text prints. --- src/utils/print_utils.cpp | 8 ++++++-- src/utils/print_utils.h | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 954af12a1..c2071f125 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -98,14 +98,18 @@ void print_utils::headline(int level, const std::string &txt) { mrcpp::print::separator(level, ' '); } -void print_utils::text(int level, const std::string &txt, const std::string &val) { +void print_utils::text(int level, const std::string &txt, const std::string &val, bool ralign) { int w0 = Printer::getWidth() - 2; int w1 = w0 * 2 / 9; int w2 = w0 - 3 * w1; int w3 = w2 - (txt.size() + 1); + // Right-aligning + int shift; + shift = (ralign) ? w0 - w2 - val.size() - 1 : 0; + std::stringstream o; - o << " " << txt << std::string(w3, ' ') << ": " << val; + o << " " << txt << std::string(w3, ' ') << ": " << std::string(shift, ' ') << val; println(level, o.str()); } diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h index df62fa7ca..28f76f4b6 100644 --- a/src/utils/print_utils.h +++ b/src/utils/print_utils.h @@ -35,7 +35,7 @@ namespace print_utils { void centeredText(int level, const std::string &txt); void headline(int level, const std::string &txt); -void text(int level, const std::string &txt, const std::string &val); +void text(int level, const std::string &txt, const std::string &val, bool ralign = false); void json(int level, const nlohmann::json &j, bool ralign = false); void coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p = -1, bool s = false); void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false); From b4fc267f8cd67d4ec2e364b89d504b35586d64c5 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:08:27 +0200 Subject: [PATCH 22/34] Added option for right-aligning text prints. --- src/utils/print_utils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index c2071f125..1f90a0c1b 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -109,7 +109,7 @@ void print_utils::text(int level, const std::string &txt, const std::string &val shift = (ralign) ? w0 - w2 - val.size() - 1 : 0; std::stringstream o; - o << " " << txt << std::string(w3, ' ') << ": " << std::string(shift, ' ') << val; + o << " " << txt << std::string(w3, ' ') << ": " << std::string(shift, ' ') << val; println(level, o.str()); } From f003bec5e098c0c005ca3b7578c416b361d797c9 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:09:09 +0200 Subject: [PATCH 23/34] Added getter for the permittivity. --- src/environment/SCRF.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h index ec807f752..83ecf0c72 100644 --- a/src/environment/SCRF.h +++ b/src/environment/SCRF.h @@ -59,6 +59,8 @@ class SCRF final { QMFunction &getPreviousGamma() { return this->gamma_nm1; } QMFunction &getCurrentDifferenceGamma() { return this->dgamma_n; } + Permittivity &getPermittivity() { return this->epsilon; } + void updateMOResidual(double const err_t) { this->mo_residual = err_t; } void printParameters(); From 6133f6a3738f24532d463cdd525103cd6f36576a Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:10:06 +0200 Subject: [PATCH 24/34] Updated the parameter printout. --- src/environment/SCRF.cpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 60310e204..284c09aaa 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -294,15 +294,11 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { } void SCRF::printParameters() { - nlohmann::json data = { - {"Dielec. const. (in)", epsilon.getEpsIn()}, - {"Dielec. const. (out)", epsilon.getEpsOut()}, - {"Max. iterations", max_iter}, - {"Accelerate with KAIN", (accelerate_Vr) ? "true" : "false"}, - {"Algorithm", algorithm}, - {"Density type", density_type}, - {"Convergence criterion", convergence_criterion}, - }; + nlohmann::json data = {{"Max iterations", max_iter}, + {"Accelerate with KAIN", (accelerate_Vr) ? "true" : "false"}, + {"Algorithm", algorithm}, + {"Density type", density_type}, + {"Convergence criterion", convergence_criterion}}; mrcpp::print::separator(0, '~'); print_utils::centeredText(0, "Self-Consistent Reaction Field"); From 0cc4f01c8989cf2be472b1f6d2c13dbdd9bf5516 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:10:43 +0200 Subject: [PATCH 25/34] Added functionality for printint permittivity and cavity parameters. --- src/environment/Permittivity.cpp | 50 ++++++++++++++++++++++++++++++++ src/environment/Permittivity.h | 11 +++++++ 2 files changed, 61 insertions(+) diff --git a/src/environment/Permittivity.cpp b/src/environment/Permittivity.cpp index 33cf9e38d..e22ae6f0a 100644 --- a/src/environment/Permittivity.cpp +++ b/src/environment/Permittivity.cpp @@ -44,4 +44,54 @@ double Permittivity::evalf(const mrcpp::Coord<3> &r) const { } } +void Permittivity::printParameters() { + // Collect relevant quantities + Cavity cavity = getCavity(); + std::vector> coords = cavity.getCoordinates(); + std::vector radii = cavity.getRadii(); + + // Set widths + auto w0 = mrcpp::Printer::getWidth() - 1; + auto w1 = 5; + auto w2 = 12; + auto w3 = 2 * w0 / 9; + auto w4 = w0 - w1 - w2 - 3 * w3; + + // Build table column headers + std::stringstream o_head; + o_head << std::setw(w1) << "N"; + o_head << std::setw(w2) << "Radius"; + o_head << std::string(w4 - 1, ' ') << ':'; + o_head << std::setw(w3) << "x"; + o_head << std::setw(w3) << "y"; + o_head << std::setw(w3) << "z"; + + // Print + mrcpp::print::header(0, "Solvation Cavity"); + print_utils::scalar(0, "Cavity width", cavity.getWidth(), "", 6); + print_utils::scalar(0, "Dielec. Const. (in)", getEpsIn(), "", 6); + print_utils::scalar(0, "Dielec. Const. (out)", getEpsOut(), "", 6); + print_utils::text(0, "Formulation", getFormulation(), true); + mrcpp::print::separator(0, '-'); + println(0, o_head.str()); + mrcpp::print::separator(0, '-'); + for (int i = 0; i < coords.size(); i++) { + mrcpp::Coord<3> coord = coords[i]; + double r = radii[i]; + double x = coord[0]; + double y = coord[1]; + double z = coord[2]; + + std::stringstream o_coord; + o_coord << std::setw(w1) << i; + o_coord << std::setw(w2) << std::setprecision(6) << std::fixed << r; + o_coord << std::string(w4 - 1, ' ') << ':'; + o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << x; + o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << y; + o_coord << std::setw(w3) << std::setprecision(6) << std::fixed << z; + println(0, o_coord.str()); + } + mrcpp::print::separator(0, '=', 2); +} + } // namespace mrchem diff --git a/src/environment/Permittivity.h b/src/environment/Permittivity.h index a7adac7db..b2f72d5b3 100644 --- a/src/environment/Permittivity.h +++ b/src/environment/Permittivity.h @@ -26,7 +26,9 @@ #pragma once #include "Cavity.h" +#include "utils/print_utils.h" #include +#include namespace mrchem { /** @class Permittivity @@ -83,6 +85,15 @@ class Permittivity final : public mrcpp::RepresentableFunction<3> { /** @brief Returns the value of #epsilon_out. */ auto getEpsOut() const { return this->epsilon_out; } + /** @brief Returns the cavity */ + Cavity getCavity() const { return this->cavity; } + + /** @brief Returns the formulation */ + std::string getFormulation() const { return this->formulation; } + + /** @brief Print parameters */ + void printParameters(); + private: bool inverse = false; //!< State of #evalf double epsilon_in; //!< Dielectric constant describing the permittivity of free space. From 966d15655c3b675480f7f8f15de9fd3c8f8e5ffe Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:11:35 +0200 Subject: [PATCH 26/34] Print SCRF and cavity/permittivity parameters at the start of an SCF run if solvent is requested. --- src/scf_solver/GroundStateSolver.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp index 0f87f7c1b..0f91515be 100644 --- a/src/scf_solver/GroundStateSolver.cpp +++ b/src/scf_solver/GroundStateSolver.cpp @@ -220,7 +220,12 @@ void GroundStateSolver::reset() { */ json GroundStateSolver::optimize(Molecule &mol, FockBuilder &F) { printParameters("Optimize ground state orbitals"); - if (F.getReactionOperator() != nullptr) F.getReactionOperator()->getHelper()->printParameters(); + + // Print solvation cavity and SCRF parameters + if (F.getReactionOperator() != nullptr) { + F.getReactionOperator()->getHelper()->printParameters(); + F.getReactionOperator()->getHelper()->getPermittivity().printParameters(); + } Timer t_tot; json json_out; From 3bb801b128c1cf86f9beacf64c37754ca2b3dc55 Mon Sep 17 00:00:00 2001 From: Anders Brakestad Date: Wed, 18 May 2022 12:15:56 +0200 Subject: [PATCH 27/34] Use capitalized True/False in bool printout. --- src/environment/SCRF.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 284c09aaa..1c2516509 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -295,7 +295,7 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { void SCRF::printParameters() { nlohmann::json data = {{"Max iterations", max_iter}, - {"Accelerate with KAIN", (accelerate_Vr) ? "true" : "false"}, + {"Accelerate with KAIN", (accelerate_Vr) ? "True" : "False"}, {"Algorithm", algorithm}, {"Density type", density_type}, {"Convergence criterion", convergence_criterion}}; From 55e13acd1cddf5dfe5349db77630c2bf2fa8d4fd Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 11:06:13 +0200 Subject: [PATCH 28/34] Add print_level escape to print_utils::json --- src/utils/print_utils.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index 1f90a0c1b..f6ea37cb0 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -114,6 +114,7 @@ void print_utils::text(int level, const std::string &txt, const std::string &val } void print_utils::json(int level, const nlohmann::json &j, bool ralign) { + if (level > mrcpp::Printer::getPrintLevel()) return; // Determine longest name // This will be used for the spacing left of : // if the name is too long to be aligned with From 834f952343d0eccdde1deb089d833b6bf8631004 Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 11:09:25 +0200 Subject: [PATCH 29/34] Add local print level to Accelerator --- src/scf_solver/Accelerator.cpp | 28 ++++++++++++++-------------- src/scf_solver/Accelerator.h | 2 ++ src/scf_solver/KAIN.cpp | 4 ++-- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/scf_solver/Accelerator.cpp b/src/scf_solver/Accelerator.cpp index 12aa165d6..25f20fdf2 100644 --- a/src/scf_solver/Accelerator.cpp +++ b/src/scf_solver/Accelerator.cpp @@ -107,7 +107,7 @@ void Accelerator::rotate(const ComplexMatrix &U, bool all) { F = U.adjoint() * F * U; dF = U.adjoint() * dF * U; } - mrcpp::print::time(2, "Rotating iterative subspace", t_tot); + mrcpp::print::time(this->pl + 2, "Rotating iterative subspace", t_tot); } /** @brief Update iterative history with the latest orbitals and updates @@ -137,7 +137,7 @@ void Accelerator::push_back(OrbitalVector &Phi, OrbitalVector &dPhi, ComplexMatr if (historyIsFull and this->dFock.size() > 0) this->dFock.pop_front(); if (not verifyOverlap(Phi)) { - println(2, " Clearing accelerator"); + println(this->pl + 2, " Clearing accelerator"); this->clear(); } @@ -146,7 +146,7 @@ void Accelerator::push_back(OrbitalVector &Phi, OrbitalVector &dPhi, ComplexMatr if (F != nullptr) this->fock.push_back(*F); if (dF != nullptr) this->dFock.push_back(*dF); - mrcpp::print::time(2, "Push back orbitals", t_tot); + mrcpp::print::time(this->pl + 2, "Push back orbitals", t_tot); } /** @brief Verify that the orbital overlap between the two last iterations is positive. @@ -170,7 +170,7 @@ bool Accelerator::verifyOverlap(OrbitalVector &Phi) { auto sqNorm = phi_i.squaredNorm(); auto overlap = orbital::dot(phi_i, last_i); if (std::abs(overlap) < 0.5 * sqNorm) { - mrcpp::print::value(2, "Overlap not verified ", std::abs(overlap)); + mrcpp::print::value(this->pl + 2, "Overlap not verified ", std::abs(overlap)); out(i) = 1; } } @@ -201,8 +201,8 @@ void Accelerator::accelerate(double prec, OrbitalVector &Phi, OrbitalVector &dPh if (this->maxHistory < 1) return; Timer t_tot; - auto plevel = Printer::getPrintLevel(); - mrcpp::print::header(2, "Iterative subspace accelerator"); + auto plevel = Printer::getPrintLevel(); // global print level + mrcpp::print::header(this->pl + 2, "Iterative subspace accelerator"); // Deep copy into history this->push_back(Phi, dPhi, F, dF); @@ -216,8 +216,8 @@ void Accelerator::accelerate(double prec, OrbitalVector &Phi, OrbitalVector &dPh clearLinearSystem(); } printSizeNodes(); - mrcpp::print::footer(2, t_tot, 2); - if (plevel == 1) mrcpp::print::time(1, "Iterative subspace accelerator", t_tot); + mrcpp::print::footer(this->pl + 2, t_tot, 2); + if (plevel == 1) mrcpp::print::time(this->pl + 1, "Iterative subspace accelerator", t_tot); } /** @brief Solve the linear system of equations @@ -232,7 +232,7 @@ void Accelerator::solveLinearSystem() { auto tmpC = this->A[n].colPivHouseholderQr().solve(this->b[n]).eval(); this->c.push_back(tmpC); } - mrcpp::print::time(2, "Solve linear system", t_tot); + mrcpp::print::time(this->pl + 2, "Solve linear system", t_tot); } /** @brief Copy orbitals at the given point in history into the given set. @@ -249,7 +249,7 @@ void Accelerator::copyOrbitals(OrbitalVector &Phi, int nHistory) { if (nHistory >= totHistory or nHistory < 0) MSG_ABORT("Requested orbitals unavailable"); int n = totHistory - 1 - nHistory; Phi = orbital::deep_copy(this->orbitals[n]); - mrcpp::print::time(2, "Copy orbitals", t_tot); + mrcpp::print::time(this->pl + 2, "Copy orbitals", t_tot); } /** @brief Copy orbital updates at the given point in history into the given set. @@ -265,7 +265,7 @@ void Accelerator::copyOrbitalUpdates(OrbitalVector &dPhi, int nHistory) { if (nHistory >= totHistory or nHistory < 0) MSG_ABORT("Requested orbitals unavailable"); int n = totHistory - 1 - nHistory; dPhi = orbital::deep_copy(this->dOrbitals[n]); - mrcpp::print::time(2, "Copy orbital updates", t_tot); + mrcpp::print::time(this->pl + 2, "Copy orbital updates", t_tot); } /** @brief Replaces the orbital set from a given point in history. @@ -355,9 +355,9 @@ void Accelerator::printSizeNodes() const { dn += orbital::get_n_nodes(orbs_i); dm += orbital::get_size_nodes(orbs_i); } - mrcpp::print::separator(2, '-'); - mrcpp::print::tree(2, "Orbital sizes", n, m, 0.0); - mrcpp::print::tree(2, "Update sizes", dn, dm, 0.0); + mrcpp::print::separator(this->pl + 2, '-'); + mrcpp::print::tree(this->pl + 2, "Orbital sizes", n, m, 0.0); + mrcpp::print::tree(this->pl + 2, "Update sizes", dn, dm, 0.0); } } // namespace mrchem diff --git a/src/scf_solver/Accelerator.h b/src/scf_solver/Accelerator.h index 23ba8bafc..17d0d85eb 100644 --- a/src/scf_solver/Accelerator.h +++ b/src/scf_solver/Accelerator.h @@ -58,6 +58,7 @@ class Accelerator { virtual ~Accelerator() { this->clear(); } void clear(); + void setLocalPrintLevel(int p) { this->pl = p; } void setMaxHistory(int max) { this->maxHistory = max; } void setMinHistory(int min) { this->minHistory = min; } @@ -79,6 +80,7 @@ class Accelerator { void printSizeNodes() const; protected: + int pl{0}; ///< Local print level, relative to global Printer int minHistory; ///< Accelerator is activated when history reaches this size int maxHistory; ///< Oldest iteration is discarded when history exceeds this size bool sepOrbitals; ///< Use separate subspace for each orbital diff --git a/src/scf_solver/KAIN.cpp b/src/scf_solver/KAIN.cpp index 24db10c8e..16d29a58b 100644 --- a/src/scf_solver/KAIN.cpp +++ b/src/scf_solver/KAIN.cpp @@ -120,7 +120,7 @@ void KAIN::setupLinearSystem() { } sortLinearSystem(A_matrices, b_vectors); - mrcpp::print::time(2, "Setup linear system", t_tot); + mrcpp::print::time(this->pl + 2, "Setup linear system", t_tot); } /** @brief Compute the next step for orbitals and orbital updates @@ -213,7 +213,7 @@ void KAIN::expandSolution(double prec, OrbitalVector &Phi, OrbitalVector &dPhi, // F is unchanged, dF is overwritten *dF = fockStep; } - mrcpp::print::time(2, "Expand solution", t_tot); + mrcpp::print::time(this->pl + 2, "Expand solution", t_tot); } } // namespace mrchem From 2bec4a54559c8844874ebd86d794c298648885b4 Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 11:10:02 +0200 Subject: [PATCH 30/34] Add average option to get_nodes and get_sizes --- src/qmfunctions/orbital_utils.cpp | 24 ++++++++++++++---------- src/qmfunctions/orbital_utils.h | 4 ++-- 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/qmfunctions/orbital_utils.cpp b/src/qmfunctions/orbital_utils.cpp index 97eff6b2d..6be524a51 100644 --- a/src/qmfunctions/orbital_utils.cpp +++ b/src/qmfunctions/orbital_utils.cpp @@ -1593,18 +1593,22 @@ int orbital::get_electron_number(const OrbitalVector &Phi, int spin) { return nElectrons; } -/** @brief Returns the total number of nodes in the vector */ -int orbital::get_n_nodes(const OrbitalVector &Phi) { - int nNodes = 0; - for (const auto &phi_i : Phi) nNodes += phi_i.getNNodes(NUMBER::Total); - return nNodes; +/** @brief Returns the total number of nodes in the vector, toggle to get average. */ +int orbital::get_n_nodes(const OrbitalVector &Phi, bool avg) { + long long totNodes = 0; + for (const auto &phi_i : Phi) totNodes += phi_i.getNNodes(NUMBER::Total); + if (avg) totNodes /= Phi.size(); + if (totNodes > INT_MAX) MSG_WARN("Integer overflow: " << totNodes); + return static_cast(totNodes); } -/** @brief Returns the size of the coefficients of all nodes in the vector in kBytes */ -int orbital::get_size_nodes(const OrbitalVector &Phi) { - int tot_size = 0; - for (const auto &phi_i : Phi) tot_size += phi_i.getSizeNodes(NUMBER::Total); - return tot_size; +/** @brief Returns the size of the coefficients of all nodes in the vector in kBytes, toggle to get average.*/ +int orbital::get_size_nodes(const OrbitalVector &Phi, bool avg) { + long long totSize = 0; + for (const auto &phi_i : Phi) totSize += phi_i.getSizeNodes(NUMBER::Total); + if (avg) totSize /= Phi.size(); + if (totSize > INT_MAX) MSG_WARN("Integer overflow: " << totSize); + return static_cast(totSize); } /** @brief Returns a vector containing the orbital spins */ diff --git a/src/qmfunctions/orbital_utils.h b/src/qmfunctions/orbital_utils.h index a10db9adf..9e0efe2d5 100644 --- a/src/qmfunctions/orbital_utils.h +++ b/src/qmfunctions/orbital_utils.h @@ -81,8 +81,8 @@ int size_beta(const OrbitalVector &Phi); int get_multiplicity(const OrbitalVector &Phi); int get_electron_number(const OrbitalVector &Phi, int spin = SPIN::Paired); int start_index(const OrbitalVector &Phi, int spin); -int get_n_nodes(const OrbitalVector &Phi); -int get_size_nodes(const OrbitalVector &Phi); +int get_n_nodes(const OrbitalVector &Phi, bool avg = false); +int get_size_nodes(const OrbitalVector &Phi, bool avg = false); bool orbital_vector_is_sane(const OrbitalVector &Phi); void set_spins(OrbitalVector &Phi, const IntVector &spins); From 4c0f9cfcdd21a33afd5e247e7b6e3c1d0b82951b Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 13:44:50 +0200 Subject: [PATCH 31/34] Name change print::centeredText -> centered_text --- src/chemistry/PhysicalConstants.cpp | 2 +- src/utils/print_utils.cpp | 2 +- src/utils/print_utils.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/chemistry/PhysicalConstants.cpp b/src/chemistry/PhysicalConstants.cpp index afca4c3aa..bf28b9299 100644 --- a/src/chemistry/PhysicalConstants.cpp +++ b/src/chemistry/PhysicalConstants.cpp @@ -59,7 +59,7 @@ PhysicalConstants &PhysicalConstants::Initialize(const json &constants) { /** @brief Pretty print physical constants */ void PhysicalConstants::Print() { mrcpp::print::separator(0, '~'); - print_utils::centeredText(0, "Physical Constants"); + print_utils::centered_text(0, "Physical Constants"); mrcpp::print::separator(0, '~'); print_utils::json(0, constants_, true); mrcpp::print::separator(0, '~', 2); diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp index f6ea37cb0..65e4e6626 100644 --- a/src/utils/print_utils.cpp +++ b/src/utils/print_utils.cpp @@ -73,7 +73,7 @@ std::string print_utils::dbl_to_str(double d, int p, bool sci) { return o_dbl.str(); } /** @brief Prints text centered according to current print width */ -void print_utils::centeredText(int level, const std::string &txt) { +void print_utils::centered_text(int level, const std::string &txt) { int pwidth = Printer::getWidth(); int txt_width = txt.size(); int pre_spaces = (pwidth - txt_width) / 2; diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h index 28f76f4b6..2b29cc020 100644 --- a/src/utils/print_utils.h +++ b/src/utils/print_utils.h @@ -33,7 +33,7 @@ namespace mrchem { namespace print_utils { -void centeredText(int level, const std::string &txt); +void centered_text(int level, const std::string &txt); void headline(int level, const std::string &txt); void text(int level, const std::string &txt, const std::string &val, bool ralign = false); void json(int level, const nlohmann::json &j, bool ralign = false); From 8cb9933e476dda0c35ef732a8213f927a00a8807 Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 16:58:14 +0200 Subject: [PATCH 32/34] Streamline printed output --- doc/users/user_inp.rst | 8 +- src/driver.cpp | 4 +- src/environment/SCRF.cpp | 103 ++++++++++++++---- src/environment/SCRF.h | 7 +- src/initial_guess/core.cpp | 4 +- src/initial_guess/sad.cpp | 46 ++++---- src/mrdft/MRDFT.cpp | 36 ++++-- src/parallel.cpp | 2 +- src/properties/SCFEnergy.h | 21 +++- src/qmoperators/one_electron/CMakeLists.txt | 1 + src/qmoperators/one_electron/ZoraOperator.cpp | 69 ++++++++++++ src/qmoperators/one_electron/ZoraOperator.h | 30 +---- .../two_electron/CoulombPotential.cpp | 13 ++- .../two_electron/CoulombPotentialD1.cpp | 4 +- .../two_electron/CoulombPotentialD2.cpp | 4 +- .../two_electron/ExchangePotential.cpp | 15 ++- .../two_electron/ExchangePotentialD1.cpp | 66 ++++++----- .../two_electron/ExchangePotentialD2.cpp | 6 +- src/qmoperators/two_electron/FockBuilder.cpp | 66 +++++++---- .../two_electron/ReactionOperator.h | 1 + .../two_electron/ReactionPotential.cpp | 14 +++ src/qmoperators/two_electron/XCPotential.cpp | 23 +++- .../two_electron/XCPotentialD1.cpp | 6 +- .../two_electron/XCPotentialD2.cpp | 12 +- src/scf_solver/GroundStateSolver.cpp | 8 -- src/scf_solver/SCFSolver.cpp | 6 +- src/tensor/RankZeroOperator.cpp | 4 +- 27 files changed, 396 insertions(+), 183 deletions(-) create mode 100644 src/qmoperators/one_electron/ZoraOperator.cpp diff --git a/doc/users/user_inp.rst b/doc/users/user_inp.rst index 190e0e662..7893191bf 100644 --- a/doc/users/user_inp.rst +++ b/doc/users/user_inp.rst @@ -100,8 +100,10 @@ Available print levels are: - ``print_level=0`` prints mainly properties - ``print_level=1`` adds timings for individual steps - ``print_level=2`` adds memory and timing information on ``OrbitalVector`` level -- ``print_level=3`` adds memory and timing information on ``Orbital`` level -- ``print_level>10`` adds a *lot* more output from deep within MRCPP +- ``print_level=3`` adds details for individual terms of the Fock operator +- ``print_level=4`` adds memory and timing information on ``Orbital`` level +- ``print_level>=5`` adds debug information at MRChem level +- ``print_level>=10`` adds debug information at MRCPP level MPI @@ -572,7 +574,7 @@ affect the convergence rate. For response calculations, the important convergence threshold is that of the orbitals, and by default this is set one order of magnitude higher than the -overall ``world_prec``. +overall ``world_prec``. .. note:: diff --git a/src/driver.cpp b/src/driver.cpp index 284009368..0fa626e81 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -1038,13 +1038,11 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto formulation = json_fock["reaction_operator"]["formulation"]; auto density_type = json_fock["reaction_operator"]["density_type"]; Permittivity dielectric_func(*cavity_r, eps_in_r, eps_out_r, formulation); + dielectric_func.printParameters(); auto scrf_p = std::make_unique(dielectric_func, nuclei, P_r, D_r, poisson_prec, hist_r, max_iter, accelerate_Vr, convergence_criterion, algorithm, density_type); auto V_R = std::make_shared(std::move(scrf_p), Phi_p); F.getReactionOperator() = V_R; - - // dielectric_func.printParameters(); - // V_R->getHelper()->printParameters(); } /////////////////////////////////////////////////////////// //////////////////// XC Operator ////////////////////// diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 1c2516509..81eec02d8 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -27,6 +27,7 @@ #include #include +#include #include "chemistry/PhysicalConstants.h" #include "chemistry/chemistry_utils.h" @@ -65,6 +66,7 @@ SCRF::SCRF(Permittivity e, , max_iter(max_iter) , history(kain_hist) , apply_prec(orb_prec) + , conv_thrs(1.0) , mo_residual(-1.0) , epsilon(e) , rho_nuc(false) @@ -90,6 +92,14 @@ void SCRF::clear() { this->rho_tot.free(NUMBER::Real); } +double SCRF::setConvergenceThreshold(double prec) { + // converge_thrs should be in the interval [prec, 1.0] + this->conv_thrs = prec; + bool dynamic = (this->convergence_criterion == "dynamic"); + if (dynamic and this->mo_residual > 10 * prec) this->conv_thrs = std::min(1.0, this->mo_residual); + return this->conv_thrs; +} + void SCRF::setDCavity() { mrcpp::FunctionTree<3> *dx_cavity = new mrcpp::FunctionTree<3>(*MRA); mrcpp::FunctionTree<3> *dy_cavity = new mrcpp::FunctionTree<3>(*MRA); @@ -101,6 +111,7 @@ void SCRF::setDCavity() { } void SCRF::computeDensities(OrbitalVector &Phi) { + Timer timer; resetQMFunction(this->rho_tot); Density rho_el(false); density::compute(this->apply_prec, rho_el, Phi, DensityType::Total); @@ -112,6 +123,7 @@ void SCRF::computeDensities(OrbitalVector &Phi) { } else { qmfunction::add(this->rho_tot, 1.0, rho_el, 1.0, this->rho_nuc, -1.0); // probably change this into a vector } + print_utils::qmfunction(3, "Vacuum density", this->rho_tot, timer); } void SCRF::computeGamma(QMFunction &potential, QMFunction &out_gamma) { @@ -162,18 +174,18 @@ void SCRF::accelerateConvergence(QMFunction &dfunc, QMFunction &func, KAIN &kain } void SCRF::nestedSCRF(QMFunction V_vac) { - print_utils::headline(2, "Calculating Reaction Potential"); - KAIN kain(this->history); + kain.setLocalPrintLevel(10); - // converge_thrs should be in the interval [apply_prec, 1.0] - double converge_thrs = this->apply_prec; - if (this->convergence_criterion == "dynamic" and this->mo_residual > 10 * this->apply_prec) converge_thrs = std::min(1.0, this->mo_residual); + mrcpp::print::separator(3, '-'); - double update = 10.0; - for (int iter = 1; update >= converge_thrs && iter <= max_iter; iter++) { + double update = 10.0, norm = 1.0; + int iter = 1; + while (update >= this->conv_thrs && iter <= max_iter) { + Timer t_iter; // solve the poisson equation QMFunction Vr_np1 = solvePoissonEquation(this->gamma_n); + norm = Vr_np1.norm(); // use a convergence accelerator resetQMFunction(this->dVr_n); @@ -204,11 +216,10 @@ void SCRF::nestedSCRF(QMFunction V_vac) { } updateCurrentGamma(gamma_np1); - - print_utils::text(2, "Update ", print_utils::dbl_to_str(update, 5, true)); - print_utils::text(2, "Microiteration ", std::to_string(iter)); + printConvergenceRow(iter, norm, update, t_iter.elapsed()); + iter++; } - println(2, " Converged Reaction Potential!"); + mrcpp::print::separator(3, '-'); this->dVr_n.real().clear(); this->dVr_n.real().setZero(); this->dgamma_n.real().clear(); @@ -216,15 +227,42 @@ void SCRF::nestedSCRF(QMFunction V_vac) { kain.clear(); } +void SCRF::printConvergenceRow(int i, double norm, double update, double time) const { + auto pprec = Printer::getPrecision(); + auto w0 = Printer::getWidth() - 1; + auto w1 = 9; + auto w2 = 2 * w0 / 9; + auto w3 = w0 - w1 - 2 * w2; + + std::string time_unit = " sec"; + if (time < 0.01) { + time = time * 1000.0; + time_unit = " ms"; + } else if (time > 60.0) { + time = time / 60.0; + time_unit = " min"; + } + + std::stringstream o_txt; + o_txt << " Iter " << std::setw(w1 - 6) << i; + o_txt << " : " << std::setw(w3 - 3) << std::setprecision(2 * pprec) << std::scientific << norm; + o_txt << std::setw(w2) << std::setprecision(pprec) << std::scientific << update; + o_txt << std::setw(w2 - 4) << std::setprecision(2) << std::fixed << time << time_unit; + println(3, o_txt.str()); +} + QMFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) { this->apply_prec = prec; computeDensities(*Phi); + Timer t_vac; QMFunction V_vac; V_vac.alloc(NUMBER::Real); mrcpp::apply(this->apply_prec, V_vac.real(), *poisson, this->rho_tot.real()); + print_utils::qmfunction(3, "Vacuum potential", V_vac, t_vac); // set up the zero-th iteration potential and gamma, so the first iteration gamma and potentials can be made + Timer t_gamma; if ((not this->Vr_n.hasReal()) or (not this->gamma_n.hasReal())) { QMFunction gamma_0; QMFunction V_tot; @@ -251,8 +289,11 @@ QMFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) { qmfunction::deep_copy(this->gamma_n, temp_gamma_n); temp_gamma_n.free(NUMBER::Real); } + print_utils::qmfunction(3, "Initial gamma", this->gamma_n, t_gamma); + Timer t_scrf; if (this->algorithm == "scrf") nestedSCRF(V_vac); + print_utils::qmfunction(3, "Reaction potential", this->Vr_n, t_scrf); return this->Vr_n; } @@ -293,18 +334,34 @@ void SCRF::updateCurrentGamma(QMFunction &gamma_np1) { gamma_np1.free(NUMBER::Real); } -void SCRF::printParameters() { - nlohmann::json data = {{"Max iterations", max_iter}, - {"Accelerate with KAIN", (accelerate_Vr) ? "True" : "False"}, - {"Algorithm", algorithm}, - {"Density type", density_type}, - {"Convergence criterion", convergence_criterion}}; - - mrcpp::print::separator(0, '~'); - print_utils::centeredText(0, "Self-Consistent Reaction Field"); - mrcpp::print::separator(0, '~'); - print_utils::json(0, data, false); - mrcpp::print::separator(0, '~', 2); +void SCRF::printParameters() const { + bool dynamic = (this->convergence_criterion == "dynamic"); + + std::stringstream o_iter; + if (this->max_iter > 0) { + o_iter << this->max_iter; + } else { + o_iter << "Off"; + } + + std::stringstream o_kain; + if (this->history > 0) { + o_kain << this->history; + } else { + o_kain << "Off"; + } + + nlohmann::json data = { + {"Method ", "SCRF"}, + {"Max iterations ", max_iter}, + {"KAIN solver ", o_kain.str()}, + {"Density type ", density_type}, + {"Dynamic threshold ", (dynamic) ? "On" : "Off"}, + }; + + mrcpp::print::separator(3, '~'); + print_utils::json(3, data, false); + mrcpp::print::separator(3, '~'); } } // namespace mrchem diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h index 83ecf0c72..8c006b396 100644 --- a/src/environment/SCRF.h +++ b/src/environment/SCRF.h @@ -51,6 +51,8 @@ class SCRF final { ~SCRF(); void UpdateExternalDensity(Density new_density) { this->rho_ext = new_density; } + double setConvergenceThreshold(double prec); + QMFunction &getCurrentReactionPotential() { return this->Vr_n; } QMFunction &getPreviousReactionPotential() { return this->Vr_nm1; } QMFunction &getCurrentDifferenceReactionPotential() { return this->dVr_n; } @@ -62,7 +64,6 @@ class SCRF final { Permittivity &getPermittivity() { return this->epsilon; } void updateMOResidual(double const err_t) { this->mo_residual = err_t; } - void printParameters(); friend class ReactionPotential; @@ -78,6 +79,7 @@ class SCRF final { int max_iter; int history; double apply_prec; + double conv_thrs; double mo_residual; Permittivity epsilon; @@ -117,5 +119,8 @@ class SCRF final { void resetQMFunction(QMFunction &function); void updateCurrentReactionPotential(QMFunction &Vr_np1); void updateCurrentGamma(QMFunction &gamma_np1); + + void printParameters() const; + void printConvergenceRow(int i, double norm, double update, double time) const; }; } // namespace mrchem diff --git a/src/initial_guess/core.cpp b/src/initial_guess/core.cpp index 80372b630..ddbd8f8bf 100644 --- a/src/initial_guess/core.cpp +++ b/src/initial_guess/core.cpp @@ -134,7 +134,7 @@ bool initial_guess::core::setup(OrbitalVector &Phi, double prec, const Nuclei &n V.clear(); p.clear(); - mrcpp::print::footer(2, t_lap, 2); + mrcpp::print::footer(2, t_tot, 2); if (plevel == 1) mrcpp::print::footer(1, t_tot, 2); return true; } @@ -258,10 +258,12 @@ void initial_guess::core::rotate_orbitals(OrbitalVector &Psi, double prec, Compl ComplexMatrix initial_guess::core::diagonalize(OrbitalVector &Phi, MomentumOperator &p, RankZeroOperator &V) { Timer t1; ComplexMatrix S_m12 = orbital::calc_lowdin_matrix(Phi); + mrcpp::print::separator(2, '-'); ComplexMatrix t_tilde = qmoperator::calc_kinetic_matrix(p, Phi, Phi); ComplexMatrix v_tilde = V(Phi, Phi); ComplexMatrix f_tilde = t_tilde + v_tilde; ComplexMatrix f = S_m12.adjoint() * f_tilde * S_m12; + mrcpp::print::separator(2, '-'); mrcpp::print::time(1, "Computing Fock matrix", t1); Timer t2; diff --git a/src/initial_guess/sad.cpp b/src/initial_guess/sad.cpp index ba6709e3b..9f41cf110 100644 --- a/src/initial_guess/sad.cpp +++ b/src/initial_guess/sad.cpp @@ -115,11 +115,11 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c initial_guess::core::project_ao(Psi, prec, nucs, zeta); if (plevel == 1) mrcpp::print::time(1, "Projecting Hydrogen AOs", t_lap); - mrcpp::print::header(2, "Building Fock operator"); + if (plevel == 2) mrcpp::print::header(2, "Building Fock operator"); t_lap.start(); p.setup(prec); V.setup(prec); - mrcpp::print::footer(2, t_lap, 2); + if (plevel == 2) mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Building Fock operator", t_lap); // Compute Fock matrix @@ -138,7 +138,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c p.clear(); V.clear(); - mrcpp::print::footer(2, t_lap, 2); + mrcpp::print::footer(2, t_tot, 2); if (plevel == 1) mrcpp::print::footer(1, t_tot, 2); return true; } @@ -194,12 +194,11 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c OrbitalVector Psi; initial_guess::gto::project_ao(Psi, prec, nucs); if (plevel == 1) mrcpp::print::time(1, "Projecting GTO AOs", t_lap); - - mrcpp::print::header(2, "Building Fock operator"); + if (plevel == 2) mrcpp::print::header(2, "Building Fock operator"); t_lap.start(); p.setup(prec); V.setup(prec); - mrcpp::print::footer(2, t_lap, 2); + if (plevel == 2) mrcpp::print::footer(2, t_lap, 2); if (plevel == 1) mrcpp::print::time(1, "Building Fock operator", t_lap); // Compute Fock matrix @@ -218,7 +217,7 @@ bool initial_guess::sad::setup(OrbitalVector &Phi, double prec, double screen, c p.clear(); V.clear(); - mrcpp::print::footer(2, t_lap, 2); + mrcpp::print::footer(2, t_tot, 2); if (plevel == 1) mrcpp::print::footer(1, t_tot, 2); return true; } @@ -257,11 +256,10 @@ void initial_guess::sad::project_atomic_densities(double prec, Density &rho_tot, rho_loc.alloc(NUMBER::Real); rho_loc.real().setZero(); - auto tot_nuc = 0.0; - auto tot_rho = 0.0; - Timer t_loc; - for (int k = 0; k < nucs.size(); k++) { + auto N_nucs = nucs.size(); + DoubleVector charges = DoubleVector::Zero(2 * N_nucs); + for (int k = 0; k < N_nucs; k++) { if (mpi::orb_rank != k % mpi::orb_size) continue; const std::string &sym = nucs[k].getElement().getSymbol(); @@ -273,24 +271,28 @@ void initial_guess::sad::project_atomic_densities(double prec, Density &rho_tot, rho_loc.add(1.0, rho_k); rho_loc.crop(crop_prec); - auto nuc_charge = nucs[k].getCharge(); - auto rho_charge = rho_k.integrate().real(); - tot_nuc += nuc_charge; - tot_rho += rho_charge; + charges[k] = nucs[k].getCharge(); + charges[N_nucs + k] = rho_k.integrate().real(); + } + t_loc.stop(); + + Timer t_com; + mpi::allreduce_vector(charges, mpi::comm_orb); + density::allreduce_density(prec, rho_tot, rho_loc); + t_com.stop(); + for (int k = 0; k < N_nucs; k++) { std::stringstream o_row; o_row << std::setw(w1) << k; - o_row << std::setw(w2) << sym; + o_row << std::setw(w2) << nucs[k].getElement().getSymbol(); o_row << std::setw(w4) << " "; - o_row << std::setw(w3) << print_utils::dbl_to_str(nuc_charge, 2 * pprec, false); - o_row << std::setw(w3) << print_utils::dbl_to_str(rho_charge, 2 * pprec, false); + o_row << std::setw(w3) << print_utils::dbl_to_str(charges[k], 2 * pprec, false); + o_row << std::setw(w3) << print_utils::dbl_to_str(charges[N_nucs + k], 2 * pprec, false); println(2, o_row.str()); } - t_loc.stop(); - Timer t_com; - density::allreduce_density(prec, rho_tot, rho_loc); - t_com.stop(); + auto tot_nuc = charges.segment(0, N_nucs).sum(); + auto tot_rho = charges.segment(N_nucs, N_nucs).sum(); std::stringstream o_row; o_row << " Total charge"; diff --git a/src/mrdft/MRDFT.cpp b/src/mrdft/MRDFT.cpp index 4a443f12c..0b7d03851 100644 --- a/src/mrdft/MRDFT.cpp +++ b/src/mrdft/MRDFT.cpp @@ -60,12 +60,22 @@ namespace mrdft { * out_vec[2] = v_xc_b (XC beta potential) */ mrcpp::FunctionTreeVector<3> MRDFT::evaluate(mrcpp::FunctionTreeVector<3> &inp) { - mrcpp::Timer timer; + mrcpp::Timer t_tot, t_pre; grid().unify(inp); functional().preprocess(inp); mrcpp::FunctionTreeVector<3> xcInpVec = functional().setupXCInput(); mrcpp::FunctionTreeVector<3> ctrInpVec = functional().setupCtrInput(); + int inpNodes = 0; + int inpSize = 0; + for (auto i = 0; i < xcInpVec.size(); i++) { + auto &f_i = mrcpp::get_func(xcInpVec, i); + inpNodes += f_i.getNNodes(); + inpSize += f_i.getSizeNodes(); + } + mrcpp::print::tree(3, "Preprocess input", inpNodes, inpSize, t_pre.elapsed()); + + mrcpp::Timer t_eval; int nCoefs = mrcpp::get_func(inp, 0).getEndFuncNode(0).getNCoefs(); int nOutCtr = functional().getCtrOutputLength(); int nFcs = functional().getXCOutputLength(); @@ -135,22 +145,34 @@ mrcpp::FunctionTreeVector<3> MRDFT::evaluate(mrcpp::FunctionTreeVector<3> &inp) */ // Reconstruct contracted output functions - int totNodes = 0; - int totSize = 0; + int ctrNodes = 0; + int ctrSize = 0; for (auto i = 0; i < ctrOutVec.size(); i++) { auto &f_i = mrcpp::get_func(ctrOutVec, i); f_i.mwTransform(mrcpp::BottomUp); f_i.calcSquareNorm(); // TODO? insert a crop - totNodes += f_i.getNNodes(); - totSize += f_i.getSizeNodes(); + ctrNodes += f_i.getNNodes(); + ctrSize += f_i.getSizeNodes(); } + mrcpp::print::tree(3, "Evaluate functional", ctrNodes, ctrSize, t_eval.elapsed()); + mrcpp::Timer t_post; auto potOutVec = functional().postprocess(ctrOutVec); mrcpp::clear(ctrOutVec, true); functional().clear(); - auto t = timer.elapsed(); - mrcpp::print::tree(2, "XC evaluate", totNodes, totSize, t); + + int outNodes = 0; + int outSize = 0; + for (auto i = 0; i < potOutVec.size(); i++) { + auto &f_i = mrcpp::get_func(potOutVec, i); + f_i.mwTransform(mrcpp::BottomUp); + f_i.calcSquareNorm(); + // TODO? insert a crop + outNodes += f_i.getNNodes(); + outSize += f_i.getSizeNodes(); + } + mrcpp::print::tree(3, "Postprocess potential", outNodes, outSize, t_post.elapsed()); return potOutVec; } diff --git a/src/parallel.cpp b/src/parallel.cpp index 223f1d74b..370b19521 100644 --- a/src/parallel.cpp +++ b/src/parallel.cpp @@ -196,7 +196,7 @@ void mpi::initialize() { void mpi::finalize() { #ifdef MRCHEM_HAS_MPI if (mpi::bank_size > 0 and mpi::grand_master()) { - println(3, " max data in bank " << dataBank.get_maxtotalsize() << " MB "); + println(4, " max data in bank " << dataBank.get_maxtotalsize() << " MB "); dataBank.close(); } MPI_Barrier(MPI_COMM_WORLD); // to ensure everybody got here diff --git a/src/properties/SCFEnergy.h b/src/properties/SCFEnergy.h index c82dd6b78..2a9a79b4e 100644 --- a/src/properties/SCFEnergy.h +++ b/src/properties/SCFEnergy.h @@ -74,6 +74,9 @@ class SCFEnergy final { auto E_kJ = E_au * PhysicalConstants::get("hartree2kjmol"); auto E_kcal = E_au * PhysicalConstants::get("hartree2kcalmol"); + bool has_ext = (std::abs(E_eext) > mrcpp::MachineZero) || (std::abs(E_next) > mrcpp::MachineZero); + bool has_react = (std::abs(Er_el) > mrcpp::MachineZero) || (std::abs(Er_nuc) > mrcpp::MachineZero); + auto pprec = 2 * mrcpp::Printer::getPrecision(); mrcpp::print::header(0, "Molecular Energy (" + id + ")"); print_utils::scalar(0, "Kinetic energy ", E_kin, "(au)", pprec, false); @@ -81,13 +84,19 @@ class SCFEnergy final { print_utils::scalar(0, "Coulomb energy ", E_ee, "(au)", pprec, false); print_utils::scalar(0, "Exchange energy ", E_x, "(au)", pprec, false); print_utils::scalar(0, "X-C energy ", E_xc, "(au)", pprec, false); - print_utils::scalar(0, "Ext. field (el) ", E_eext, "(au)", pprec, false); - mrcpp::print::separator(0, '-'); print_utils::scalar(0, "N-N energy ", E_nn, "(au)", pprec, false); - print_utils::scalar(0, "Ext. field (nuc) ", E_next, "(au)", pprec, false); - print_utils::scalar(0, "Tot. Reac. Energy", Er_tot, "(au)", pprec, false); - print_utils::scalar(0, "El. Reac. Energy ", Er_el, "(au)", pprec, false); - print_utils::scalar(0, "Nuc. Reac. Energy", Er_nuc, "(au)", pprec, false); + if (has_ext) { + mrcpp::print::separator(0, '-'); + print_utils::scalar(0, "External field (el) ", E_eext, "(au)", pprec, false); + print_utils::scalar(0, "External field (nuc) ", E_next, "(au)", pprec, false); + print_utils::scalar(0, "External field (tot) ", E_eext + E_next, "(au)", pprec, false); + } + if (has_react) { + mrcpp::print::separator(0, '-'); + print_utils::scalar(0, "Reaction energy (el) ", Er_el, "(au)", pprec, false); + print_utils::scalar(0, "Reaction energy (nuc) ", Er_nuc, "(au)", pprec, false); + print_utils::scalar(0, "Reaction energy (tot) ", Er_tot, "(au)", pprec, false); + } mrcpp::print::separator(0, '-'); print_utils::scalar(0, "Electronic energy", E_el, "(au)", pprec, false); print_utils::scalar(0, "Nuclear energy ", E_nuc, "(au)", pprec, false); diff --git a/src/qmoperators/one_electron/CMakeLists.txt b/src/qmoperators/one_electron/CMakeLists.txt index 2291d40b7..264297579 100644 --- a/src/qmoperators/one_electron/CMakeLists.txt +++ b/src/qmoperators/one_electron/CMakeLists.txt @@ -1,3 +1,4 @@ target_sources(mrchem PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/NuclearOperator.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ZoraOperator.cpp ) diff --git a/src/qmoperators/one_electron/ZoraOperator.cpp b/src/qmoperators/one_electron/ZoraOperator.cpp new file mode 100644 index 000000000..a31b071e3 --- /dev/null +++ b/src/qmoperators/one_electron/ZoraOperator.cpp @@ -0,0 +1,69 @@ +/* + * MRChem, a numerical real-space code for molecular electronic structure + * calculations within the self-consistent field (SCF) approximations of quantum + * chemistry (Hartree-Fock and Density Functional Theory). + * Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. + * + * This file is part of MRChem. + * + * MRChem is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRChem is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRChem. If not, see . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#include "ZoraOperator.h" + +#include +#include + +#include "qmfunctions/qmfunction_utils.h" +#include "qmoperators/QMPotential.h" +#include "utils/print_utils.h" + +using mrcpp::Printer; +using mrcpp::Timer; + +namespace mrchem { + +ZoraOperator::ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse) { + Timer timer; + double two_cc = 2.0 * c * c; + + auto k = std::make_shared(1); + qmfunction::deep_copy(*k, vz); + + if (k->hasImag()) MSG_ERROR("Inverse of complex function"); + if (k->hasReal()) { + mrcpp::refine_grid(k->real(), 1); + if (inverse) { + k->real().map([two_cc](double val) { return (two_cc - val) / two_cc; }); + } else { + k->real().map([two_cc](double val) { return two_cc / (two_cc - val); }); + } + k->real().crop(proj_prec); + } + + RankZeroOperator &kappa = (*this); + kappa = k; + if (inverse) { + kappa.name() = "kappa_m1"; + } else { + kappa.name() = "kappa"; + } + auto plevel = Printer::getPrintLevel(); + print_utils::qmfunction(2, "ZORA operator (" + kappa.name() + ")", *k, timer); +} + +} // namespace mrchem diff --git a/src/qmoperators/one_electron/ZoraOperator.h b/src/qmoperators/one_electron/ZoraOperator.h index d9abdbb73..0130b0c2e 100644 --- a/src/qmoperators/one_electron/ZoraOperator.h +++ b/src/qmoperators/one_electron/ZoraOperator.h @@ -25,39 +25,15 @@ #pragma once -#include "qmfunctions/qmfunction_utils.h" -#include "qmoperators/QMPotential.h" #include "tensor/RankZeroOperator.h" namespace mrchem { +class QMPotential; + class ZoraOperator final : public RankZeroOperator { public: - ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse = false) { - double two_cc = 2.0 * c * c; - - auto k = std::make_shared(1); - qmfunction::deep_copy(*k, vz); - - if (k->hasImag()) MSG_ERROR("Inverse of complex function"); - if (k->hasReal()) { - mrcpp::refine_grid(k->real(), 1); - if (inverse) { - k->real().map([two_cc](double val) { return (two_cc - val) / two_cc; }); - } else { - k->real().map([two_cc](double val) { return two_cc / (two_cc - val); }); - } - k->real().crop(proj_prec); - } - - RankZeroOperator &kappa = (*this); - kappa = k; - if (inverse) { - kappa.name() = "kappa_m1"; - } else { - kappa.name() = "kappa"; - } - } + ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse = false); }; } // namespace mrchem diff --git a/src/qmoperators/two_electron/CoulombPotential.cpp b/src/qmoperators/two_electron/CoulombPotential.cpp index 4799944c5..fb9f016db 100644 --- a/src/qmoperators/two_electron/CoulombPotential.cpp +++ b/src/qmoperators/two_electron/CoulombPotential.cpp @@ -74,6 +74,11 @@ CoulombPotential::CoulombPotential(PoissonOperator_p P, OrbitalVector_p Phi, boo void CoulombPotential::setup(double prec) { if (isSetup(prec)) return; setApplyPrec(prec); + Timer timer; + auto plevel = Printer::getPrintLevel(); + mrcpp::print::header(3, "Building Coulomb operator"); + mrcpp::print::value(3, "Precision", prec, "(rel)", 5); + mrcpp::print::separator(3, '-'); if (hasDensity()) { setupGlobalPotential(prec); } else if (mpi::numerically_exact) { @@ -86,6 +91,8 @@ void CoulombPotential::setup(double prec) { QMFunction V = setupLocalPotential(0.1 * prec); allreducePotential(0.1 * prec, V); } + if (plevel == 2) print_utils::qmfunction(2, "Coulomb operator", *this, timer); + mrcpp::print::footer(3, timer, 2); } /** @brief clear operator after application @@ -124,7 +131,7 @@ void CoulombPotential::setupGlobalPotential(double prec) { V.alloc(NUMBER::Real); if (need_to_apply) mrcpp::apply(abs_prec, V.real(), P, rho.real()); mpi::share_function(V, 0, 22445, mpi::comm_share); - print_utils::qmfunction(2, "Coulomb potential", V, timer); + print_utils::qmfunction(3, "Compute global potential", V, timer); } /** @brief compute Coulomb potential @@ -148,7 +155,7 @@ QMFunction CoulombPotential::setupLocalPotential(double prec) { QMFunction V(false); V.alloc(NUMBER::Real); mrcpp::apply(abs_prec, V.real(), P, rho.real()); - print_utils::qmfunction(2, "Coulomb potential", V, timer); + print_utils::qmfunction(3, "Compute local potential", V, timer); return V; } @@ -183,7 +190,7 @@ void CoulombPotential::allreducePotential(double prec, QMFunction &V_loc) { mrcpp::copy_grid(V_tot.real(), V_loc.real()); mrcpp::copy_func(V_tot.real(), V_loc.real()); } - print_utils::qmfunction(2, "Allreduce Coulomb", V_tot, t_com); + print_utils::qmfunction(3, "Allreduce potential", V_tot, t_com); } } // namespace mrchem diff --git a/src/qmoperators/two_electron/CoulombPotentialD1.cpp b/src/qmoperators/two_electron/CoulombPotentialD1.cpp index 9f2ce960d..4989adb33 100644 --- a/src/qmoperators/two_electron/CoulombPotentialD1.cpp +++ b/src/qmoperators/two_electron/CoulombPotentialD1.cpp @@ -53,7 +53,7 @@ void CoulombPotentialD1::setupGlobalDensity(double prec) { Timer timer; density::compute(prec, rho, Phi, DensityType::Total); - print_utils::qmfunction(2, "Coulomb density", rho, timer); + print_utils::qmfunction(3, "Compute global density", rho, timer); } /** @brief compute local electron density (only own MPI orbitals) @@ -71,7 +71,7 @@ void CoulombPotentialD1::setupLocalDensity(double prec) { Timer timer; density::compute_local(prec, rho, Phi, DensityType::Total); - print_utils::qmfunction(2, "Coulomb density", rho, timer); + print_utils::qmfunction(3, "Compute local density", rho, timer); } } // namespace mrchem diff --git a/src/qmoperators/two_electron/CoulombPotentialD2.cpp b/src/qmoperators/two_electron/CoulombPotentialD2.cpp index added2e14..ade5356a3 100644 --- a/src/qmoperators/two_electron/CoulombPotentialD2.cpp +++ b/src/qmoperators/two_electron/CoulombPotentialD2.cpp @@ -57,7 +57,7 @@ void CoulombPotentialD2::setupGlobalDensity(double prec) { Timer timer; density::compute(prec, rho, Phi, X, Y, DensityType::Total); - print_utils::qmfunction(2, "Coulomb density", rho, timer); + print_utils::qmfunction(3, "Compute global density", rho, timer); } void CoulombPotentialD2::setupLocalDensity(double prec) { @@ -73,7 +73,7 @@ void CoulombPotentialD2::setupLocalDensity(double prec) { Timer timer; density::compute_local(prec, rho, Phi, X, Y, DensityType::Total); - print_utils::qmfunction(2, "Coulomb density", rho, timer); + print_utils::qmfunction(3, "Compute local density", rho, timer); } } // namespace mrchem diff --git a/src/qmoperators/two_electron/ExchangePotential.cpp b/src/qmoperators/two_electron/ExchangePotential.cpp index d8747e2b7..1fbfcf7be 100644 --- a/src/qmoperators/two_electron/ExchangePotential.cpp +++ b/src/qmoperators/two_electron/ExchangePotential.cpp @@ -117,10 +117,23 @@ double ExchangePotential::getSpinFactor(Orbital phi_i, Orbital phi_j) const { * and optionally compute the internal constributions. */ void ExchangePotential::setup(double prec) { + Timer timer; + auto plevel = Printer::getPrintLevel(); + mrcpp::print::header(3, "Building Exchange operator"); + mrcpp::print::value(3, "Precision", prec, "(rel)", 5); + mrcpp::print::separator(3, '-'); if (mpi::world_size > 1 and mpi::bank_size < 1) MSG_ABORT("MPI bank required!"); setApplyPrec(prec); setupBank(); if (this->pre_compute) setupInternal(prec); + + if (plevel == 2) { + auto t = timer.elapsed(); + auto avgNodes = orbital::get_n_nodes(this->exchange, true); + auto avgSize = orbital::get_size_nodes(this->exchange, true); + mrcpp::print::tree(2, "Exchange operator (avg.)", avgNodes, avgSize, t); + } + mrcpp::print::footer(3, timer, 2); } /** @brief Clears the Exchange Operator @@ -215,7 +228,7 @@ void ExchangePotential::calcExchange_kij(double prec, Orbital phi_k, Orbital phi } timer_jji.stop(); - println(4, + println(5, " time " << (int)((float)timer_tot.elapsed() * 1000) << " ms " << " mult1:" << (int)((float)timer_ij.elapsed() * 1000) << " Pot:" << (int)((float)timer_p.elapsed() * 1000) << " mult2:" << (int)((float)timer_kij.elapsed() * 1000) << " " << (int)((float)timer_jji.elapsed() * 1000) << " Nnodes: " << N_i << " " << N_j << " " << N_ij << " " << N_p << " " << N_kij << " " << N_jji << " norms " << norm_ij << " " diff --git a/src/qmoperators/two_electron/ExchangePotentialD1.cpp b/src/qmoperators/two_electron/ExchangePotentialD1.cpp index 8a910ea4a..5c27c7768 100644 --- a/src/qmoperators/two_electron/ExchangePotentialD1.cpp +++ b/src/qmoperators/two_electron/ExchangePotentialD1.cpp @@ -64,7 +64,7 @@ void ExchangePotentialD1::setupBank() { if (mpi::my_orb(Phi[i])) PhiBank.put_orb(i, Phi[i]); } mpi::barrier(mpi::comm_orb); - mrcpp::print::time(4, "Setting up exchange bank", timer); + mrcpp::print::time(3, "Setting up exchange bank", timer); } /** @brief Clears rbital bank accounts. @@ -116,10 +116,8 @@ Orbital ExchangePotentialD1::apply(Orbital phi_p) { MSG_WARN("Not computing exchange contributions that are not mine"); return out_p; } - println(4, "On-the-fly exchange"); return calcExchange(phi_p); } else { - println(4, "Precomputed exchange"); return this->exchange[i]; } } @@ -144,7 +142,7 @@ QMOperatorVector ExchangePotentialD1::apply(QMOperator_p &O) { * The exchange potential is (pre)computed among the orbitals that define the operator */ void ExchangePotentialD1::setupInternal(double prec) { - Timer timerT, timerS(false), timerR(false), t_calc(false), t_add(false), t_get(false), t_wait(false); + Timer t_tot, t_snd(false), t_orb(false), t_calc(false), t_add(false), t_get(false), t_wait(false); setApplyPrec(prec); if (this->exchange.size() != 0) MSG_ERROR("Exchange not properly cleared"); @@ -157,7 +155,7 @@ void ExchangePotentialD1::setupInternal(double prec) { prec = mpi::numerically_exact ? -1.0 : prec; precf /= std::sqrt(1 * Phi.size()); // Initialize this->exchange and compute own diagonal elements - Timer timerD; + Timer t_diag; for (auto &phi_i : Phi) { Orbital ex_iii = phi_i.paramCopy(); t_calc.resume(); @@ -165,8 +163,9 @@ void ExchangePotentialD1::setupInternal(double prec) { t_calc.stop(); Ex.push_back(ex_iii); } - mrcpp::print::time(4, "Exchange time diagonal", timerD); + t_diag.stop(); + Timer t_offd; // We divide all the exchange contributions into a fixed number of tasks. // all "j" orbitals are fetched and stored, and used together with one "i" orbital // At the end of a task, the Exchange for j is summed up with the contributions @@ -250,25 +249,25 @@ void ExchangePotentialD1::setupInternal(double prec) { i0 = iorb; Orbital phi_i; - timerR.resume(); + t_orb.resume(); if (bank_size > 0) { PhiBank.get_orb(iorb, phi_i, 1); // fetch also own orbitals (simpler for clean up, and they are few) iorb_vec.push_back(phi_i); } else { iorb_vec.push_back(Phi[iorb]); } - timerR.stop(); + t_orb.stop(); } for (int j = 0; j < jtasks[task].size(); j++) { int jorb = jtasks[task][j]; Orbital phi_j; - timerR.resume(); + t_orb.resume(); if (bank_size > 0) PhiBank.get_orb(jorb, phi_j, 1); else phi_j = Phi[jorb]; - timerR.stop(); + t_orb.stop(); QMFunctionVector iijfunc_vec; ComplexVector coef_vec(N); for (int i = 0; i < iorb_vec.size(); i++) { @@ -284,7 +283,7 @@ void ExchangePotentialD1::setupInternal(double prec) { calcExchange_kij(precf, phi_i, phi_i, phi_j, ex_iij, &ex_jji); t_calc.stop(); if (ex_iij.norm() > prec) coef_vec[iijfunc_vec.size()] = j_fac; - timerS.resume(); + t_snd.resume(); if (bank_size > 0) { // store ex_jji if (ex_iij.norm() > prec) iijfunc_vec.push_back(ex_iij); @@ -295,7 +294,7 @@ void ExchangePotentialD1::setupInternal(double prec) { Ex[jorb].add(j_fac, ex_iij); } ex_jji.free(NUMBER::Total); - timerS.stop(); + t_snd.stop(); } Timer timerx; // fetch ready contributions to ex_j from others @@ -320,14 +319,14 @@ void ExchangePotentialD1::setupInternal(double prec) { // ex_j is sent to Bank if (ex_j.hasReal() or ex_j.hasImag()) { auto tT = timerx.elapsed(); - timerS.resume(); + t_snd.resume(); ex_j.crop(prec); if (ex_j.norm() > prec) { if (i0 < 0) MSG_ERROR("no exchange contributions?"); ExBank.put_orb(jorb + (i0 + N) * N, ex_j); tasksMaster.put_readytask(jorb, i0 + N); } - timerS.stop(); + t_snd.stop(); ex_j.free(NUMBER::Total); } else if (iijfunc_vec.size() > 0) { MSG_ERROR("Exchange exists but has no real and no Imag parts"); @@ -342,7 +341,6 @@ void ExchangePotentialD1::setupInternal(double prec) { mpi::barrier(mpi::comm_orb); t_wait.stop(); - IntVector sizes = IntVector::Zero(2 * N); for (int j = 0; j < N; j++) { if (not mpi::my_orb(Phi[j]) or bank_size == 0) continue; // fetch only own j std::vector iVec = tasksMaster.get_readytask(j, 1); @@ -385,26 +383,24 @@ void ExchangePotentialD1::setupInternal(double prec) { Ex[j].crop(prec); t_add.stop(); } - sizes[j] = Ex[j].getNNodes(NUMBER::Total); - sizes[j + N] = Ex[j].getSizeNodes(NUMBER::Total); } - mrcpp::print::time(4, "Time rcv orbitals", timerR); - mrcpp::print::time(4, "Time send exchanges", timerS); - mrcpp::print::time(4, "Time rcv exchanges", t_get); - mrcpp::print::time(4, "Time wait others finished", t_wait); - mrcpp::print::time(4, "Time add exchanges", t_add); - mrcpp::print::time(4, "Time calculate exchanges", t_calc); - - auto t = timerT.elapsed(); - mpi::allreduce_vector(sizes, mpi::comm_orb); - long long nsum = 0; - for (int j = 0; j < N; j++) nsum += sizes[j]; - long long msum = 0; - for (int j = 0; j < N; j++) msum += sizes[j + N]; - int n = nsum / N; - int m = msum / N; - mrcpp::print::tree(2, "HF exchange (av.)", n, m, t); -} // namespace mrchem + t_offd.stop(); + mrcpp::print::time(3, "Time receiving orbitals", t_orb); + mrcpp::print::time(3, "Time receiving exchanges", t_get); + mrcpp::print::time(3, "Time sending exchanges", t_snd); + mrcpp::print::time(3, "Time adding exchanges", t_add); + mrcpp::print::time(3, "Time waiting for others", t_wait); + mrcpp::print::time(3, "Time computing exchanges", t_calc); + mrcpp::print::separator(3, '-'); + mrcpp::print::time(3, "Time diagonal terms", t_diag); + mrcpp::print::time(3, "Time off-diagonal terms", t_offd); + mrcpp::print::separator(3, '-'); + + auto t = t_tot.elapsed() / N; + auto n = orbital::get_n_nodes(this->exchange, true); + auto m = orbital::get_size_nodes(this->exchange, true); + mrcpp::print::tree(3, "Average exchange term", n, m, t); +} /** @brief Computes the exchange potential on the fly * @@ -444,7 +440,7 @@ Orbital ExchangePotentialD1::calcExchange(Orbital phi_p) { Orbital ex_p = phi_p.paramCopy(); Eigen::Map coefs(coef_vec.data(), coef_vec.size()); qmfunction::linear_combination(ex_p, coefs, func_vec, prec); - print_utils::qmfunction(3, "Applied exchange", ex_p, timer); + print_utils::qmfunction(4, "Applied exchange", ex_p, timer); return ex_p; } diff --git a/src/qmoperators/two_electron/ExchangePotentialD2.cpp b/src/qmoperators/two_electron/ExchangePotentialD2.cpp index 4533ccbfa..169719967 100644 --- a/src/qmoperators/two_electron/ExchangePotentialD2.cpp +++ b/src/qmoperators/two_electron/ExchangePotentialD2.cpp @@ -76,7 +76,7 @@ void ExchangePotentialD2::setupBank() { if (mpi::my_orb(Y[i])) YBank.put_orb(i, Y[i]); } mpi::barrier(mpi::comm_orb); - mrcpp::print::time(4, "Setting up exchange bank", timer); + mrcpp::print::time(3, "Setting up exchange bank", timer); } /** @brief Clears the Exchange Operator @@ -144,7 +144,7 @@ Orbital ExchangePotentialD2::apply(Orbital phi_p) { Orbital out_p = phi_p.paramCopy(); Eigen::Map coefs(coef_vec.data(), coef_vec.size()); qmfunction::linear_combination(out_p, coefs, func_vec, prec); - print_utils::qmfunction(3, "Applied exchange", out_p, timer); + print_utils::qmfunction(4, "Applied exchange", out_p, timer); return out_p; } @@ -203,7 +203,7 @@ Orbital ExchangePotentialD2::dagger(Orbital phi_p) { Orbital ex_p = phi_p.paramCopy(); Eigen::Map coefs(coef_vec.data(), coef_vec.size()); qmfunction::linear_combination(ex_p, coefs, func_vec, prec); - print_utils::qmfunction(3, "Applied exchange", ex_p, timer); + print_utils::qmfunction(4, "Applied exchange", ex_p, timer); return ex_p; } diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp index 4e6b9cdd9..91a03a93b 100644 --- a/src/qmoperators/two_electron/FockBuilder.cpp +++ b/src/qmoperators/two_electron/FockBuilder.cpp @@ -75,15 +75,22 @@ void FockBuilder::build(double exx) { void FockBuilder::setup(double prec) { Timer t_tot; auto plevel = Printer::getPrintLevel(); - mrcpp::print::header(2, "Building Fock operator"); - mrcpp::print::value(2, "Precision", prec, "(rel)", 5); - mrcpp::print::separator(2, '-'); + if (plevel == 2) { + mrcpp::print::header(2, "Building Fock operator"); + mrcpp::print::value(2, "Precision", prec, "(rel)", 5); + mrcpp::print::separator(2, '-'); + } if (this->mom != nullptr) this->momentum().setup(prec); this->potential().setup(prec); this->perturbation().setup(prec); if (isZora()) { + Timer t_zora; auto c = getLightSpeed(); + mrcpp::print::header(3, "Building ZORA operators"); + mrcpp::print::value(3, "Precision", prec, "(rel)", 5); + mrcpp::print::value(3, "Light speed", c, "(au)", 5); + mrcpp::print::separator(3, '-'); auto vz = collectZoraBasePotential(); this->kappa = std::make_shared(*vz, c, prec, false); this->kappa_inv = std::make_shared(*vz, c, prec, true); @@ -91,10 +98,11 @@ void FockBuilder::setup(double prec) { this->kappa->setup(prec); this->kappa_inv->setup(prec); this->zora_base.setup(prec); + mrcpp::print::footer(3, t_zora, 2); }; t_tot.stop(); - mrcpp::print::footer(2, t_tot, 2); + if (plevel == 2) mrcpp::print::footer(2, t_tot, 2); if (plevel == 1) mrcpp::print::time(1, "Building Fock operator", t_tot); } @@ -203,12 +211,13 @@ ComplexMatrix FockBuilder::operator()(OrbitalVector &bra, OrbitalVector &ket) { } OrbitalVector FockBuilder::buildHelmholtzArgument(double prec, OrbitalVector Phi, ComplexMatrix F_mat, ComplexMatrix L_mat) { - Timer t_arg; + Timer t_tot; auto plevel = Printer::getPrintLevel(); mrcpp::print::header(2, "Computing Helmholtz argument"); + Timer t_rot; OrbitalVector Psi = orbital::rotate(Phi, L_mat - F_mat, prec); - mrcpp::print::time(2, "Rotating orbitals", t_arg); + mrcpp::print::time(2, "Rotating orbitals", t_rot); OrbitalVector out; if (isZora()) { @@ -218,8 +227,8 @@ OrbitalVector FockBuilder::buildHelmholtzArgument(double prec, OrbitalVector Phi } Psi.clear(); - mrcpp::print::footer(2, t_arg, 2); - if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_arg); + mrcpp::print::footer(2, t_tot, 2); + if (plevel == 1) mrcpp::print::time(1, "Computing Helmholtz argument", t_tot); return out; } @@ -239,17 +248,24 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita operOne.setup(prec); operThree.setup(prec); + // Compute OrbitalVectors + Timer t_1; + OrbitalVector termOne = operOne(Phi); + mrcpp::print::time(2, "Computing gradient term", t_1); + + Timer t_2; + OrbitalVector termTwo = V(Phi); + mrcpp::print::time(2, "Computing potential term", t_2); + // Compute transformed orbitals scaled by diagonal Fock elements + Timer t_3; OrbitalVector epsPhi = orbital::deep_copy(Phi); for (int i = 0; i < epsPhi.size(); i++) { if (not mpi::my_orb(epsPhi[i])) continue; epsPhi[i].rescale(eps[i] / two_cc); } - - // Compute OrbitalVectors - OrbitalVector termOne = operOne(Phi); - OrbitalVector termTwo = V(Phi); OrbitalVector termThree = operThree(epsPhi); + mrcpp::print::time(2, "Computing rescaled potential term", t_3); auto normsOne = orbital::get_norms(termOne); auto normsTwo = orbital::get_norms(termTwo); @@ -257,17 +273,23 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita auto normsPsi = orbital::get_norms(Psi); // Add up all the terms - OrbitalVector out = orbital::deep_copy(termOne); - for (int i = 0; i < out.size(); i++) { - if (not mpi::my_orb(out[i])) continue; - out[i].add(1.0, termTwo[i]); - out[i].add(1.0, termThree[i]); - out[i].add(1.0, Psi[i]); + Timer t_add; + OrbitalVector arg = orbital::deep_copy(termOne); + for (int i = 0; i < arg.size(); i++) { + if (not mpi::my_orb(arg[i])) continue; + arg[i].add(1.0, termTwo[i]); + arg[i].add(1.0, termThree[i]); + arg[i].add(1.0, Psi[i]); } + mrcpp::print::time(2, "Adding contributions", t_add); operThree.clear(); operOne.clear(); - return kappa_m1(out); + + Timer t_kappa; + auto out = kappa_m1(arg); + mrcpp::print::time(2, "Applying kappa inverse", t_kappa); + return out; } // Non-relativistic Helmholtz argument @@ -276,14 +298,18 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentNREL(OrbitalVector &Phi, Orbita RankZeroOperator &V = this->potential(); // Compute OrbitalVectors + Timer t_pot; OrbitalVector termOne = V(Phi); + mrcpp::print::time(2, "Computing potential term", t_pot); // Add up all the terms + Timer t_add; OrbitalVector out = orbital::deep_copy(termOne); for (int i = 0; i < out.size(); i++) { if (not mpi::my_orb(out[i])) continue; out[i].add(1.0, Psi[i]); }; + mrcpp::print::time(2, "Adding contributions", t_add); return out; } @@ -294,6 +320,7 @@ void FockBuilder::setZoraType(bool has_nuc, bool has_coul, bool has_xc) { } std::shared_ptr FockBuilder::collectZoraBasePotential() { + Timer timer; auto vz = std::make_shared(1, false); if (zora_has_nuc) { if (getNuclearOperator() != nullptr) { @@ -324,6 +351,7 @@ std::shared_ptr FockBuilder::collectZoraBasePotential() { MSG_ERROR("ZORA: XC requested but not available"); } } + print_utils::qmfunction(2, "ZORA operator (base)", *vz, timer); return vz; } diff --git a/src/qmoperators/two_electron/ReactionOperator.h b/src/qmoperators/two_electron/ReactionOperator.h index 71e2eb95d..dae73f5f2 100644 --- a/src/qmoperators/two_electron/ReactionOperator.h +++ b/src/qmoperators/two_electron/ReactionOperator.h @@ -47,6 +47,7 @@ class ReactionOperator final : public RankZeroOperator { // Invoke operator= to assign *this operator RankZeroOperator &V = (*this); V = potential; + V.name() = "V_r"; } ComplexDouble trace(OrbitalVector &Phi) { return RankZeroOperator::trace(Phi); } diff --git a/src/qmoperators/two_electron/ReactionPotential.cpp b/src/qmoperators/two_electron/ReactionPotential.cpp index ea2fbedb2..15f250693 100644 --- a/src/qmoperators/two_electron/ReactionPotential.cpp +++ b/src/qmoperators/two_electron/ReactionPotential.cpp @@ -23,6 +23,9 @@ * */ +#include +#include + #include "ReactionPotential.h" #include "qmfunctions/qmfunction_utils.h" @@ -43,8 +46,19 @@ void ReactionPotential::setup(double prec) { this->first_iteration = false; return; } + auto thrs = this->helper->setConvergenceThreshold(prec); + mrcpp::Timer timer; + auto plevel = mrcpp::Printer::getPrintLevel(); + mrcpp::print::separator(3, '='); + print_utils::centered_text(3, "Building Reaction operator"); + this->helper->printParameters(); + mrcpp::print::value(3, "Precision", prec, "(rel)", 5); + mrcpp::print::value(3, "Threshold", thrs, "(abs)", 5); + mrcpp::print::separator(3, '-'); auto potential = this->helper->setup(prec, this->Phi); qmfunction::deep_copy(*this, potential); + if (plevel == 2) print_utils::qmfunction(2, "Reaction operator", *this, timer); + mrcpp::print::footer(3, timer, 2); } void ReactionPotential::clear() { diff --git a/src/qmoperators/two_electron/XCPotential.cpp b/src/qmoperators/two_electron/XCPotential.cpp index e9be1e998..d57255f3d 100644 --- a/src/qmoperators/two_electron/XCPotential.cpp +++ b/src/qmoperators/two_electron/XCPotential.cpp @@ -23,9 +23,10 @@ * */ +#include +#include + #include "XCPotential.h" -#include "MRCPP/Printer" -#include "MRCPP/Timer" #include "parallel.h" #include "qmfunctions/Density.h" #include "qmfunctions/Orbital.h" @@ -56,6 +57,11 @@ namespace mrchem { void XCPotential::setup(double prec) { if (isSetup(prec)) return; setApplyPrec(prec); + Timer timer; + auto plevel = Printer::getPrintLevel(); + mrcpp::print::header(3, "Building XC operator"); + mrcpp::print::value(3, "Precision", prec, "(rel)", 5); + mrcpp::print::separator(3, '-'); if (this->mrdft == nullptr) MSG_ERROR("XCFunctional not initialized"); if (this->potentials.size() != 0) MSG_ERROR("Potential not properly cleared"); @@ -84,6 +90,19 @@ void XCPotential::setup(double prec) { } mrcpp::clear(xc_out, true); + + if (plevel == 2) { + int totNodes = 0; + int totSize = 0; + for (auto i = 0; i < this->potentials.size(); i++) { + auto &f_i = mrcpp::get_func(this->potentials, i); + totNodes += f_i.getNNodes(); + totSize += f_i.getSizeNodes(); + } + auto t = timer.elapsed(); + mrcpp::print::tree(2, "XC operator", totNodes, totSize, t); + } + mrcpp::print::footer(3, timer, 2); } /** @brief Clears all data in the XCPotential object */ diff --git a/src/qmoperators/two_electron/XCPotentialD1.cpp b/src/qmoperators/two_electron/XCPotentialD1.cpp index 14cf3f895..54a58995b 100644 --- a/src/qmoperators/two_electron/XCPotentialD1.cpp +++ b/src/qmoperators/two_electron/XCPotentialD1.cpp @@ -58,7 +58,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD1::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, DensityType::Total); } - print_utils::qmfunction(2, "XC rho", rho, timer); + print_utils::qmfunction(3, "Compute rho", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } } else { @@ -70,7 +70,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD1::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, DensityType::Alpha); } - print_utils::qmfunction(2, "XC rho (alpha)", rho, timer); + print_utils::qmfunction(3, "Compute rho (alpha)", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } { // Unperturbed beta density @@ -81,7 +81,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD1::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, DensityType::Beta); } - print_utils::qmfunction(2, "XC rho (beta)", rho, timer); + print_utils::qmfunction(3, "Compute rho (beta)", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } } diff --git a/src/qmoperators/two_electron/XCPotentialD2.cpp b/src/qmoperators/two_electron/XCPotentialD2.cpp index 62e0779f5..877f598a3 100644 --- a/src/qmoperators/two_electron/XCPotentialD2.cpp +++ b/src/qmoperators/two_electron/XCPotentialD2.cpp @@ -80,7 +80,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, DensityType::Total); } - print_utils::qmfunction(2, "XC rho_0", rho, timer); + print_utils::qmfunction(3, "Compute rho_0", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } { // Perturbed total density @@ -91,7 +91,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, *orbitals_x, *orbitals_y, DensityType::Total); } - print_utils::qmfunction(2, "XC rho_1", rho, timer); + print_utils::qmfunction(3, "Compute rho_1", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } } else { @@ -103,7 +103,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, DensityType::Alpha); } - print_utils::qmfunction(2, "XC rho_0 (alpha)", rho, timer); + print_utils::qmfunction(3, "Compute rho_0 (alpha)", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } { // Unperturbed beta density @@ -114,7 +114,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, DensityType::Beta); } - print_utils::qmfunction(2, "XC rho_0 (beta)", rho, timer); + print_utils::qmfunction(3, "Compute rho_0 (beta)", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } { // Perturbed alpha density @@ -125,7 +125,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, *orbitals_x, *orbitals_y, DensityType::Alpha); } - print_utils::qmfunction(2, "XC rho_1 (alpha)", rho, timer); + print_utils::qmfunction(3, "Compute rho_1 (alpha)", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } { // Perturbed beta density @@ -136,7 +136,7 @@ mrcpp::FunctionTreeVector<3> XCPotentialD2::setupDensities(double prec, mrcpp::F mrcpp::copy_grid(rho.real(), grid); density::compute(prec, rho, *orbitals, *orbitals_x, *orbitals_y, DensityType::Beta); } - print_utils::qmfunction(2, "XC rho_1 (beta)", rho, timer); + print_utils::qmfunction(3, "Compute rho_1 (beta)", rho, timer); dens_vec.push_back(std::make_tuple(1.0, &rho.real())); } } diff --git a/src/scf_solver/GroundStateSolver.cpp b/src/scf_solver/GroundStateSolver.cpp index 0f91515be..05b7a31a8 100644 --- a/src/scf_solver/GroundStateSolver.cpp +++ b/src/scf_solver/GroundStateSolver.cpp @@ -167,8 +167,6 @@ void GroundStateSolver::printParameters(const std::string &calculation) const { o_prec_0 << std::setprecision(5) << std::scientific << this->orbPrec[1]; o_prec_1 << std::setprecision(5) << std::scientific << this->orbPrec[2]; - mrcpp::print::separator(0, '~'); - print_utils::centeredText(0, "Self-Consistent Field"); mrcpp::print::separator(0, '~'); print_utils::text(0, "Calculation ", calculation); print_utils::text(0, "Method ", this->methodName); @@ -221,12 +219,6 @@ void GroundStateSolver::reset() { json GroundStateSolver::optimize(Molecule &mol, FockBuilder &F) { printParameters("Optimize ground state orbitals"); - // Print solvation cavity and SCRF parameters - if (F.getReactionOperator() != nullptr) { - F.getReactionOperator()->getHelper()->printParameters(); - F.getReactionOperator()->getHelper()->getPermittivity().printParameters(); - } - Timer t_tot; json json_out; diff --git a/src/scf_solver/SCFSolver.cpp b/src/scf_solver/SCFSolver.cpp index de2b16e4a..152129f7d 100644 --- a/src/scf_solver/SCFSolver.cpp +++ b/src/scf_solver/SCFSolver.cpp @@ -82,10 +82,10 @@ double SCFSolver::adjustPrecision(double error) { this->orbPrec[0] = std::max(this->orbPrec[0], this->orbPrec[2]); mrcpp::print::separator(2, '='); - mrcpp::print::value(1, "Current precision", this->orbPrec[0], "", 5); + mrcpp::print::value(1, "Current precision", this->orbPrec[0], "(rel)", 5); mrcpp::print::separator(1, '-'); - mrcpp::print::value(2, "Orbital threshold", this->orbThrs, "", 5); - mrcpp::print::value(2, "Property threshold", this->propThrs, "", 5); + mrcpp::print::value(2, "Orbital threshold", this->orbThrs, "(abs)", 5); + mrcpp::print::value(2, "Property threshold", this->propThrs, "(abs)", 5); mrcpp::print::separator(2, '=', 2); return this->orbPrec[0]; } diff --git a/src/tensor/RankZeroOperator.cpp b/src/tensor/RankZeroOperator.cpp index e1f89c099..656c3a7b3 100644 --- a/src/tensor/RankZeroOperator.cpp +++ b/src/tensor/RankZeroOperator.cpp @@ -292,7 +292,7 @@ OrbitalVector RankZeroOperator::operator()(OrbitalVector &inp) { out.push_back(out_i); std::stringstream o_name; o_name << O.name() << "|" << i << ">"; - print_utils::qmfunction(3, o_name.str(), out_i, t1); + print_utils::qmfunction(4, o_name.str(), out_i, t1); } return out; } @@ -312,7 +312,7 @@ OrbitalVector RankZeroOperator::dagger(OrbitalVector &inp) { out.push_back(out_i); std::stringstream o_name; o_name << O.name() << "^dagger|" << i << ">"; - print_utils::qmfunction(3, o_name.str(), out_i, t1); + print_utils::qmfunction(4, o_name.str(), out_i, t1); } return out; } From 3a42623ca0f142fff1ea32fda0f3d80ae0a41fec Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 15:37:08 +0200 Subject: [PATCH 33/34] Remove unnecessary Lowdin in core guess --- src/initial_guess/core.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/initial_guess/core.cpp b/src/initial_guess/core.cpp index ddbd8f8bf..3b2a37463 100644 --- a/src/initial_guess/core.cpp +++ b/src/initial_guess/core.cpp @@ -112,7 +112,6 @@ bool initial_guess::core::setup(OrbitalVector &Phi, double prec, const Nuclei &n t_lap.start(); OrbitalVector Psi; initial_guess::core::project_ao(Psi, prec, nucs, zeta); - ComplexMatrix S_m12 = orbital::calc_lowdin_matrix(Psi); if (plevel == 1) mrcpp::print::time(1, "Projecting Hydrogen AOs", t_lap); p.setup(prec); From bd52a5b824fdb6f7d0366ea50e1694ef489ce30b Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Thu, 7 Jul 2022 17:06:10 +0200 Subject: [PATCH 34/34] Update gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index c1b170cb2..cd73b1d7c 100644 --- a/.gitignore +++ b/.gitignore @@ -128,3 +128,6 @@ Pipfile.lock #Containers *.sif + +#Visual Studio Code +.vscode