new mad-freq_plot_vtk executable

src/apps/molresponse/Plot_VTK.cc

@@ -65,7 +65,6 @@ namespace madness {
         std::filesystem::create_directories(vtk_dir);
 
         std::string geo_file;
-        const char *filename;
 
         Vector<long, 3> points{npt_plot, npt_plot, npt_plot};
         // Plot the whole box?
@@ -83,12 +82,13 @@ namespace madness {
 
         int orb_num = 0;
 
+        auto orbital_file = vtk_dir + "/" + "orbitals.vts";
+        plotvtk_begin<3>(world, orbital_file.c_str(), box_lo, box_hi, points, true);
+
         std::for_each(ground_orbs.begin(), ground_orbs.end(), [&](const auto &phi0_i) {
-            auto orb_file = vtk_dir + "/phi0_" + std::to_string(orb_num) + ".vts";
-            filename = orb_file.c_str();
-            plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true);
-            plotvtk_data<double, 3>(phi0_i, "ground_orbtial", world, filename, box_lo, box_hi, points, true, false);
-            plotvtk_end<3>(world, filename, true);
+            auto orb_name = "/phi0_" + std::to_string(orb_num);
+            plotvtk_data<double, 3>(phi0_i, orb_name.c_str(), world, orbital_file.c_str(), box_lo, box_hi, points, true,
+                                    false);
             orb_num++;
         });
 
@@ -99,26 +99,21 @@ namespace madness {
             // plot the x first
             auto orb_num = 0;
             std::for_each(xy.begin(), xy.begin() + num_orbitals, [&](const auto &xi) {
-                auto orb_file = vtk_dir + "/" + "x" + std::to_string(state_number) + std::to_string(orb_num) + ".vts";
-                filename = orb_file.c_str();
-                plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true);
-                plotvtk_data<double, 3>(xi, "response_x_orbitals", world, filename, box_lo, box_hi, points, true,
-                                        false);
-                plotvtk_end<3>(world, filename, true);
+                auto field_name = "x_orbital_" + std::to_string(state_number) + "_" + std::to_string(orb_num);
+                plotvtk_data<double, 3>(xi, field_name.c_str(), world, orbital_file.c_str(), box_lo, box_hi, points,
+                                        true, false);
                 orb_num++;
             });
             orb_num = 0;
             std::for_each(xy.begin() + num_orbitals, xy.end(), [&](const auto &yi) {
-                auto orb_file = vtk_dir + "/" + "y" + std::to_string(state_number) + std::to_string(orb_num) + ".vts";
-                filename = orb_file.c_str();
-                plotvtk_begin<3>(world, filename, box_lo, box_hi, points, true);
-                plotvtk_data<double, 3>(yi, "response_y_orbitals", world, filename, box_lo, box_hi, points, true,
-                                        false);
-                plotvtk_end<3>(world, filename, true);
+                auto field_name = "y_orbital_" + std::to_string(state_number) + "_" + std::to_string(orb_num);
+                plotvtk_data<double, 3>(yi, field_name.c_str(), world, orbital_file.c_str(), box_lo, box_hi, points,
+                                        true, false);
                 orb_num++;
             });
             state_number++;
         });
+        plotvtk_end<3>(world, orbital_file.c_str(), true);
     }
     void do_response_density_vtk_plots(World &world, int npt_plot, double L, const Molecule &molecule,
                                        const real_function_3d &ground_density,
@@ -162,8 +157,9 @@ namespace madness {
                                            const vector_real_function_3d &response_density) {
         // Stuff needed to plot
         //
-        Vector<double, 3> box_lo{-L / 4, -L / 4, -L / 4};
-        Vector<double, 3> box_hi{L / 4, L / 4, L / 4};
+        auto L_box = L / 10;
+        Vector<double, 3> box_lo{-L_box, -L_box, -L_box};
+        Vector<double, 3> box_hi{L_box , L_box, L_box};
 
 
         std::string vtk_dir = "vtk_plots";
@@ -183,9 +179,9 @@ namespace madness {
         //***********************************ground density plot
         int state_number = 0;
         std::for_each(response_density.begin(), response_density.end(), [&](const auto &rho_i) {
-          auto field_name = "response_" + std::to_string(state_number);
-          plotvtk_data<double, 3>(rho_i, field_name.c_str(), world, filename, box_lo, box_hi, points, true, false);
-          state_number++;
+            auto field_name = "response_" + std::to_string(state_number);
+            plotvtk_data<double, 3>(rho_i, field_name.c_str(), world, filename, box_lo, box_hi, points, true, false);
+            state_number++;
         });
         plotvtk_end<3>(world, filename, true);
     }

src/apps/molresponse/ResponseBase.cpp

@@ -1364,6 +1364,81 @@ void ResponseBase::function_data_to_json(json &j_mol_in, size_t iter, const Tens
     j_mol_in["protocol_data"][index]["iter_data"].push_back(j);
 }
 
+void ResponseBase::write_vtk(World &world) {
+    // Get start time
+    molresponse::start_timer(world);
+
+    // Plotting input orbitals
+    //if (r_params.plot_initial()) { PlotGroundDensityVTK(world, *this); }
+    const auto protocol = r_params.protocol();
+    if (world.rank() == 0) {
+        print("Response State Calculation for the following protocols");
+        print("Protocol: ", protocol);
+    }
+    bool first_protocol = true;
+    // access the last protocol only
+    auto iter_thresh = protocol.back();
+    // We set the protocol and function defaults here for the given threshold of
+    set_protocol(world, iter_thresh);
+    if (world.rank() == 0) { print("Successfully set protocol"); }
+    // protocol
+    if (first_protocol) {
+        if (r_params.restart()) {
+            if (world.rank() == 0) { print("   Restarting from file:", r_params.restart_file()); }
+            load(world, r_params.restart_file());
+            first_protocol = false;
+        } else {
+            this->initialize(world);
+            if (world.rank() == 0) { print("Successfully initialized "); }
+        }
+        check_k(world, iter_thresh, FunctionDefaults<3>::get_k());
+        if (world.rank() == 0) { print("Successfully check K first initialization "); }
+        first_protocol = false;
+    } else {
+        check_k(world, iter_thresh, FunctionDefaults<3>::get_k());
+        if (world.rank() == 0) { print("Successfully check K not first initialization "); }
+    }
+    // Now actually ready to iterate...
+    // At this point we should know if calc converged maybe add a flag to response.json which states if it has
+    // check norms of vectors ground and response
+    auto x_norms = Chi.norm2s();
+    if (world.rank() == 0) {
+        print("x_norms");
+        print(x_norms);
+    }
+    auto ground_orbitals_norms = norm2s(world, ground_orbitals);
+    if (world.rank() == 0) {
+        print("rho_norms");
+        print(ground_orbitals_norms);
+    }
+    auto response_densities = make_density(world, Chi);
+    auto rho_norms = norm2s(world, response_densities);
+    if (world.rank() == 0) {
+        print("rho_norms");
+        print(rho_norms);
+    }
+
+#if defined(__has_include)
+#if __has_include(<filesystem>)
+#define MADCHEM_HAS_STD_FILESYSTEM
+// <filesystem> is not reliably usable on Linux with gcc < 9
+#if defined(__GNUC__)
+#if __GNUC__ >= 7 && __GNUC__ < 9
+#undef MADCHEM_HAS_STD_FILESYSTEM
+#endif
+#endif
+#if defined(MADCHEM_HAS_STD_FILESYSTEM)
+    auto r_matrix = to_response_matrix(Chi);
+    do_response_orbital_vtk_plots(world, r_params.plot_pts(), r_params.L(), molecule, ground_orbitals, r_matrix);
+    do_response_density_vtk_plots_new(world, r_params.plot_pts(), r_params.L(), molecule, ground_density,
+                                      response_densities);
+#endif
+#endif
+#endif
+
+
+    // Plot the response function if desired
+}
 void ResponseBase::solve(World &world) {
     // Get start time
     molresponse::start_timer(world);

src/apps/molresponse/ResponseBase.hpp

@@ -555,6 +555,8 @@ class ResponseBase {
     }
 
 
+    void write_vtk(World &world);
+
 protected:
     // Given molecule returns the nuclear potential of the molecule
     ResponseParameters r_params;

src/apps/molresponse/response_parameters.h

@@ -38,7 +38,7 @@ namespace madness {
             initialize<std::vector<double>>("plot_cell", std::vector<double>(),
                                             "lo-hi plot cell (default is all space)");
             initialize<double>("plot_l", -1.0, "Controls the plotting box size");
-            initialize<size_t>("plot_pts", 51, "Controls number of points in plots");
+            initialize<size_t>("plot_pts", 81, "Controls number of points in plots");
             initialize<bool>("plot_all_orbitals", false, "Turn on 2D plotting of response orbitals ");
             initialize<size_t>("maxiter", 25, "maximum number of iterations");
             initialize<double>("dconv", 1.e-4, "recommended values: 1.e-4 < dconv < 1.e-8");

src/apps/molresponse/testing/CMakeLists.txt

@@ -3,6 +3,7 @@
 
 set(MY_EXECUTABLES
         mad-freq
+        mad-freq_vtk_plots
         mad-excited
         frequency_calc
         excited_state_calc

src/apps/molresponse/testing/mad-freq_vtk_plots.cpp

@@ -0,0 +1,80 @@
+//
+// Created by adrianhurtado on 1/1/22.
+//
+#include "ResponseExceptions.hpp"
+#include "madness/tensor/tensor_json.hpp"
+#include "madness/world/worldmem.h"
+#include "response_functions.h"
+#include "runners.hpp"
+
+#if defined(HAVE_SYS_TYPES_H) && defined(HAVE_SYS_STAT_H) && defined(HAVE_UNISTD_H)
+
+#include <sys/stat.h>
+#include <unistd.h>
+
+static inline int file_exists(const char *input_name) {
+    struct stat buffer {};
+    size_t rc = stat(input_name, &buffer);
+    return (rc == 0);
+}
+
+#endif
+
+using path = std::filesystem::path;
+
+using namespace madness;
+
+
+auto main(int argc, char *argv[]) -> int {
+
+    madness::initialize(argc, argv);
+
+    {
+        World world(SafeMPI::COMM_WORLD);
+        startup(world, argc, argv, true);
+
+        try {
+            sleep(5);
+            std::cout.precision(6);
+            if (argc != 6) {
+                std::cout << "Wrong number of inputs" << std::endl;
+                return 1;
+            }
+            const std::string molecule_name{argv[1]};
+            const std::string xc{argv[2]};
+            const std::string op{argv[3]};
+            const std::string precision{argv[4]};
+            const std::string static_calc{argv[5]};
+            if (precision != "high" && precision != "low" && precision != "super") {
+                if (world.rank() == 0) { std::cout << "Set precision to low high super" << std::endl; }
+                return 1;
+            }
+            auto schema = runSchema(world, xc);
+            auto m_schema = moldftSchema(world, molecule_name, xc, schema);
+            auto f_schema = frequencySchema(world, schema, m_schema, op, static_calc == "true");
+            write_VTK_outputs(world, f_schema, precision);
+        } catch (const SafeMPI::Exception &e) {
+            print(e.what());
+            error("caught an MPI exception");
+        } catch (const madness::MadnessException &e) {
+            print(e);
+            error("caught a MADNESS exception");
+        } catch (const madness::TensorException &e) {
+            print(e.what());
+            error("caught a Tensor exception");
+        } catch (const nlohmann::detail::exception &e) {
+            print(e.what());
+            error("Caught JSON exception");
+        } catch (const std::filesystem::filesystem_error &ex) {
+            std::cerr << ex.what() << "\n";
+        } catch (const std::exception &e) {
+            print(e.what());
+            error("caught an STL exception");
+        } catch (...) { error("caught unhandled exception"); }
+        // Nearly all memory will be freed at this point
+        print_stats(world);
+        if (world.rank() == 0) { print("Finished All Frequencies"); }
+    }
+    finalize();
+    return 0;
+}

src/apps/molresponse/testing/runners.hpp

@@ -696,6 +696,56 @@ auto RunResponse(World &world, const std::string &filename, double frequency, co
     return {save_path, calc.j_molresponse["converged"]};
 }
 
+/**
+ *
+ * @param world
+ * @param filename
+ * @param frequency
+ * @param property
+ * @param xc
+ * @param moldft_path
+ * @param restart_path
+ * @return
+ */
+auto WriteVTKOutputs(World &world, const std::string &filename, double frequency, const std::string &property,
+                     const std::string &xc, const std::filesystem::path &moldft_path,
+                     std::filesystem::path restart_path, const std::string &precision)
+        -> std::pair<std::filesystem::path, bool> {
+    // Set the response parameters
+    ResponseParameters r_params{};
+    set_frequency_response_parameters(world, r_params, property, xc, frequency, precision);
+    auto save_path =
+            set_frequency_path_and_restart(world, r_params, property, frequency, xc, moldft_path, restart_path, true);
+    if (world.rank() == 0) { molresponse::write_response_input(r_params, filename); }
+    // if rbase exists and converged I just return save path and true
+    auto calc_params = initialize_calc_params(world, std::string(filename));
+    RHS_Generator rhs_generator;
+    if (property == "dipole") {
+        rhs_generator = dipole_generator;
+    } else {
+        rhs_generator = nuclear_generator;
+    }
+    FunctionDefaults<3>::set_pmap(pmapT(new LevelPmap<Key<3>>(world)));
+    FrequencyResponse calc(world, calc_params, frequency, rhs_generator);
+    if (world.rank() == 0) {
+        print("\n\n");
+        print(" MADNESS Time-Dependent Density Functional Theory Response "
+              "Program");
+        print(" ----------------------------------------------------------\n");
+        print("\n");
+        calc_params.molecule.print();
+        print("\n");
+        calc_params.response_parameters.print("response");
+        // put the response parameters in a j_molrespone json object
+        calc_params.response_parameters.to_json(calc.j_molresponse);
+    }
+    calc.write_vtk(world);
+    world.gop.fence();
+    // set protocol to the first
+    //calc.time_data.print_data();
+    return {save_path, true};
+}
+
 /***
  * sets the run path based on the run type set by r_params
  * creates the run directory and sets current directory to the run data
@@ -819,6 +869,41 @@ void runFrequencyTests(World &world, const frequencySchema &schema, const std::s
     }
 }
 
+/**
+ * Takes in the moldft path where moldft restart file exists
+ * runs a response calculations for given property at given frequencies.
+ *
+ *
+ * @param world
+ * @param moldft_path
+ * @param frequencies
+ * @param xc
+ * @param property
+ */
+void write_VTK_outputs(World &world, const frequencySchema &schema, const std::string &high_prec) {
+    std::filesystem::current_path(schema.moldft_path);
+    // add a restart path
+    auto restart_path =
+            addPath(schema.moldft_path, "/" + schema.op + "_0-000000.00000/restart_" + schema.op + "_0-000000.00000");
+    std::pair<std::filesystem::path, bool> success{schema.moldft_path, false};
+    bool first = true;
+    for (const auto &freq: schema.freq) {
+        if (world.rank() == 0) { print(success.second); }
+        std::filesystem::current_path(schema.moldft_path);
+        if (first) {
+            first = false;
+        } else if (success.second) {
+            // if the previous run succeeded then set the restart path
+            restart_path = success.first;
+            if (world.rank() == 0) { print("restart_path", restart_path); }
+        } else {
+            throw Response_Convergence_Error{};
+        }
+        success = WriteVTKOutputs(world, "response.in", freq, schema.op, schema.xc, schema.moldft_path, restart_path,
+                                  high_prec);
+        if (world.rank() == 0) { print("Frequency ", freq, " completed"); }
+    }
+}
 
 // set up a function that creates a beta_json with the fields defined  below.  in each field there will
 // a vector of values.
@@ -852,9 +937,6 @@ void append_to_beta_json(const std::array<double, 3> &freq, const std::array<dou
     // capitalize the direction
 
 
-
-
-
     for (int i = 0; i < 18; i++) {
         beta_json["A-freq"].push_back(freq[0]);
         beta_json["B-freq"].push_back(freq[1]);