Merge pull request #408 from Andersmb/solvent-printouts

Print information about the external environments

.gitignore

@@ -128,3 +128,6 @@ Pipfile.lock
 
 #Containers
 *.sif
+
+#Visual Studio Code
+.vscode

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::
 

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
-
-

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<mrcpp::Coord<3>> coords = cavity.getCoordinates();
+    std::vector<double> 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);

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<mrcpp::Coord<3>> &coords, std::vector<double> &R, double slope);

src/chemistry/PhysicalConstants.cpp

@@ -58,9 +58,11 @@ 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::centered_text(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

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) {
@@ -222,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);
@@ -240,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);
@@ -255,9 +257,12 @@ 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"];
+        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 +280,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);
@@ -1032,9 +1038,10 @@ 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<SCRF>(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<ReactionOperator>(scrf_p, Phi_p);
+        auto V_R = std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
         F.getReactionOperator() = V_R;
     }
     ///////////////////////////////////////////////////////////

src/environment/Cavity.h

@@ -50,6 +50,7 @@ class Cavity final : public mrcpp::RepresentableFunction<3> {
     auto getGradVector() const { return this->gradvector; }
     std::vector<mrcpp::Coord<3>> getCoordinates() const { return centers; } //!< Returns #centers.
     std::vector<double> getRadii() const { return radii; }                  //!< Returns #radii.
+    double getWidth() { return width; }                                     //!< Returns #width.
     friend class Permittivity;
 
 protected:

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<mrcpp::Coord<3>> coords = cavity.getCoordinates();
+    std::vector<double> 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

src/environment/Permittivity.h

@@ -26,7 +26,9 @@
 #pragma once
 
 #include "Cavity.h"
+#include "utils/print_utils.h"
 #include <MRCPP/MWFunctions>
+#include <MRCPP/Printer>
 
 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.