diff --git a/src/apps/molresponse/Plot_VTK.cc b/src/apps/molresponse/Plot_VTK.cc index 0818ab255d2..14faa7b50f2 100644 --- a/src/apps/molresponse/Plot_VTK.cc +++ b/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 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(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(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(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(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(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(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 box_lo{-L / 4, -L / 4, -L / 4}; - Vector box_hi{L / 4, L / 4, L / 4}; + auto L_box = L / 10; + Vector box_lo{-L_box, -L_box, -L_box}; + Vector 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(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(rho_i, field_name.c_str(), world, filename, box_lo, box_hi, points, true, false); + state_number++; }); plotvtk_end<3>(world, filename, true); } diff --git a/src/apps/molresponse/ResponseBase.cpp b/src/apps/molresponse/ResponseBase.cpp index cd6ee7cae29..8ea514d66a3 100644 --- a/src/apps/molresponse/ResponseBase.cpp +++ b/src/apps/molresponse/ResponseBase.cpp @@ -1294,6 +1294,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() +#define MADCHEM_HAS_STD_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); diff --git a/src/apps/molresponse/ResponseBase.hpp b/src/apps/molresponse/ResponseBase.hpp index 622eed4616e..6031a24ca60 100644 --- a/src/apps/molresponse/ResponseBase.hpp +++ b/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; diff --git a/src/apps/molresponse/response_parameters.h b/src/apps/molresponse/response_parameters.h index 01ba38c0eb8..f3cccc677e4 100644 --- a/src/apps/molresponse/response_parameters.h +++ b/src/apps/molresponse/response_parameters.h @@ -38,7 +38,7 @@ namespace madness { initialize>("plot_cell", std::vector(), "lo-hi plot cell (default is all space)"); initialize("plot_l", -1.0, "Controls the plotting box size"); - initialize("plot_pts", 51, "Controls number of points in plots"); + initialize("plot_pts", 81, "Controls number of points in plots"); initialize("plot_all_orbitals", false, "Turn on 2D plotting of response orbitals "); initialize("maxiter", 25, "maximum number of iterations"); initialize("dconv", 1.e-4, "recommended values: 1.e-4 < dconv < 1.e-8"); diff --git a/src/apps/molresponse/testing/CMakeLists.txt b/src/apps/molresponse/testing/CMakeLists.txt index c08818d9593..22bcf4dba43 100644 --- a/src/apps/molresponse/testing/CMakeLists.txt +++ b/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 diff --git a/src/apps/molresponse/testing/mad-freq_vtk_plots.cpp b/src/apps/molresponse/testing/mad-freq_vtk_plots.cpp new file mode 100644 index 00000000000..455eeecac6f --- /dev/null +++ b/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 +#include + +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; +} diff --git a/src/apps/molresponse/testing/runners.hpp b/src/apps/molresponse/testing/runners.hpp index b9151c91074..2c9cf29baf9 100644 --- a/src/apps/molresponse/testing/runners.hpp +++ b/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 { + // 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>(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 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 &freq, const std::array