From fb275a72d463f786a46d423fc6b7f68f544bcbec Mon Sep 17 00:00:00 2001 From: Adrian Hurtado Date: Sun, 16 Jul 2023 20:33:59 -0400 Subject: [PATCH] The start of the new executable --- src/apps/molresponse/ResponseExceptions.hpp | 7 + .../molresponse/testing/mad-quadratic.cpp | 79 ++++ src/apps/molresponse/testing/runners.hpp | 378 ++++++++++++------ 3 files changed, 334 insertions(+), 130 deletions(-) create mode 100644 src/apps/molresponse/testing/mad-quadratic.cpp diff --git a/src/apps/molresponse/ResponseExceptions.hpp b/src/apps/molresponse/ResponseExceptions.hpp index 635e57a8a46..cb4156819af 100644 --- a/src/apps/molresponse/ResponseExceptions.hpp +++ b/src/apps/molresponse/ResponseExceptions.hpp @@ -34,4 +34,11 @@ class Response_Convergence_Error : public MadException { __LINE__, __FUNCTION__, __FILE__) {} }; +class Quad_Response_Error : public MadException { +public: + explicit Quad_Response_Error() + : MadException("Previous frequency response calculation did not converge ", nullptr, 25, + __LINE__, __FUNCTION__, __FILE__) {} +}; + #endif// MADNESS_RESPONSEEXCEPTIONS_HPP diff --git a/src/apps/molresponse/testing/mad-quadratic.cpp b/src/apps/molresponse/testing/mad-quadratic.cpp new file mode 100644 index 00000000000..26ac83b1d16 --- /dev/null +++ b/src/apps/molresponse/testing/mad-quadratic.cpp @@ -0,0 +1,79 @@ +// +// 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 != 5) { + 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]}; + + + auto schema = runSchema(world, xc); + auto m_schema = moldftSchema(world, molecule_name, xc, schema); + auto f_schema = frequencySchema(world, schema, m_schema, op); + + + runQuadraticResponse(world, f_schema, precision); + world.gop.fence(); + } 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); + } + finalize(); + return 0; +} diff --git a/src/apps/molresponse/testing/runners.hpp b/src/apps/molresponse/testing/runners.hpp index 03c3c0f2c43..8320b1354dc 100644 --- a/src/apps/molresponse/testing/runners.hpp +++ b/src/apps/molresponse/testing/runners.hpp @@ -100,23 +100,20 @@ struct moldftSchema { std::string mol_name; std::string xc; - moldftSchema(World &world, std::string molecule_name, std::string m_xc, - const runSchema &schema) + moldftSchema(World &world, std::string molecule_name, std::string m_xc, const runSchema &schema) : mol_name(std::move(molecule_name)), xc(std::move(m_xc)) { moldft_path = addPath(schema.xc_path, '/' + mol_name); moldft_restart = addPath(moldft_path, "/moldft.restartdata.00000"); calc_info_json_path = addPath(moldft_path, "/moldft.calc_info.json"); mol_path = addPath(schema.molecules, "/" + mol_name + ".mol"); - if (std::filesystem::exists(moldft_restart) && - std::filesystem::exists(calc_info_json_path)) { + if (std::filesystem::exists(moldft_restart) && std::filesystem::exists(calc_info_json_path)) { // if both exist, read the calc_info json std::ifstream ifs(calc_info_json_path); ifs >> calc_info_json; if (world.rank() == 0) { std::cout << "time: " << calc_info_json["time"] << std::endl; - std::cout << "MOLDFT return energy: " - << calc_info_json["return_energy"] << std::endl; + std::cout << "MOLDFT return energy: " << calc_info_json["return_energy"] << std::endl; } } if (world.rank() == 0) { print(); } @@ -140,8 +137,7 @@ struct frequencySchema { const path moldft_path; vector freq; - frequencySchema(World &world, const runSchema &run_schema, - const moldftSchema &m_schema, std::string r_operator); + frequencySchema(World &world, const runSchema &run_schema, const moldftSchema &m_schema, std::string r_operator); void print_schema() { print("Frequency Calculation"); @@ -152,11 +148,9 @@ struct frequencySchema { print("Frequencies : ", freq); } }; -frequencySchema::frequencySchema(World &world, const runSchema &run_schema, - const moldftSchema &m_schema, +frequencySchema::frequencySchema(World &world, const runSchema &run_schema, const moldftSchema &m_schema, std::string r_operator) - : mol_name(m_schema.mol_name), xc(m_schema.xc), op(std::move(r_operator)), - moldft_path(m_schema.moldft_path) { + : mol_name(m_schema.mol_name), xc(m_schema.xc), op(std::move(r_operator)), moldft_path(m_schema.moldft_path) { if (world.rank() == 0) { print_schema(); freq = run_schema.freq_json_data.get_frequencies(mol_name, xc, op); @@ -174,8 +168,7 @@ frequencySchema::frequencySchema(World &world, const runSchema &run_schema, * @param property * @return */ -size_t set_excited_states(const ResponseDataBase &response_data_base, - const std::string &molecule_name, +size_t set_excited_states(const ResponseDataBase &response_data_base, const std::string &molecule_name, const std::string &xc) { const std::string property = "excited-state"; @@ -184,8 +177,8 @@ size_t set_excited_states(const ResponseDataBase &response_data_base, return response_data_base.get_num_states(molecule_name, xc, property); } catch (json::exception &e) { std::cout << e.what() << std::endl; - std::cout << "did not find the frequency data for [" << molecule_name - << "][" << xc << "][" << property << "]\n"; + std::cout << "did not find the frequency data for [" << molecule_name << "][" << xc << "][" << property + << "]\n"; return 4; } } @@ -202,9 +195,7 @@ size_t set_excited_states(const ResponseDataBase &response_data_base, * @param xc * @return */ -std::filesystem::path -generate_excited_run_path(const std::filesystem::path &moldft_path, - const size_t &num_states) { +std::filesystem::path generate_excited_run_path(const std::filesystem::path &moldft_path, const size_t &num_states) { std::string s_num_states = std::to_string(num_states); std::string run_name = "excited-" + s_num_states; // set r_params to restart true if restart file exist @@ -245,12 +236,9 @@ struct excitedSchema { path rb_json; - excitedSchema(const runSchema &run_schema, const moldftSchema &m_schema) - : xc(m_schema.xc) { - num_states = set_excited_states(run_schema.freq_json_data, - m_schema.mol_name, xc); - excited_state_run_path = - generate_excited_run_path(m_schema.moldft_path, num_states); + excitedSchema(const runSchema &run_schema, const moldftSchema &m_schema) : xc(m_schema.xc) { + num_states = set_excited_states(run_schema.freq_json_data, m_schema.mol_name, xc); + excited_state_run_path = generate_excited_run_path(m_schema.moldft_path, num_states); auto [sp, s] = generate_excited_save_path(excited_state_run_path); save_path = sp; save_string = s; @@ -278,9 +266,7 @@ struct excitedSchema { * @param xc * @return xc_path */ -std::filesystem::path -create_xc_path_and_directory(const std::filesystem::path &root, - const std::string &xc) { +std::filesystem::path create_xc_path_and_directory(const std::filesystem::path &root, const std::string &xc) { // copy construct the root path auto xc_path = std::filesystem::path(root); @@ -306,8 +292,7 @@ create_xc_path_and_directory(const std::filesystem::path &root, * @param frequency_run_path * @return */ -auto generate_frequency_save_path( - const std::filesystem::path &frequency_run_path) +auto generate_frequency_save_path(const std::filesystem::path &frequency_run_path) -> std::pair { auto save_path = std::filesystem::path(frequency_run_path); @@ -332,10 +317,8 @@ auto generate_frequency_save_path( * @param xc * @return */ -auto generate_response_frequency_run_path( - const std::filesystem::path &moldft_path, const std::string &property, - const double &frequency, const std::string &xc) - -> std::filesystem::path { +auto generate_response_frequency_run_path(const std::filesystem::path &moldft_path, const std::string &property, + const double &frequency, const std::string &xc) -> std::filesystem::path { std::string s_frequency = std::to_string(frequency); auto sp = s_frequency.find('.'); s_frequency = s_frequency.replace(sp, sp, "-"); @@ -357,8 +340,8 @@ auto generate_response_frequency_run_path( * @param moldft_filename * @param xc */ -void runMOLDFT(World &world, const moldftSchema &moldftSchema, bool try_run, - bool restart, const std::string &precision) { +void runMOLDFT(World &world, const moldftSchema &moldftSchema, bool try_run, bool restart, + const std::string &precision) { CalculationParameters param1; json calcInfo; @@ -370,12 +353,10 @@ void runMOLDFT(World &world, const moldftSchema &moldftSchema, bool try_run, param1.set_user_defined_value("l", 50); if (precision == "low") { - param1.set_user_defined_value>("protocol", - {1e-4, 1e-6}); + param1.set_user_defined_value>("protocol", {1e-4, 1e-6}); param1.set_user_defined_value("dconv", 1e-4); } else if (precision == "high") { - param1.set_user_defined_value>("protocol", - {1e-4, 1e-6, 1e-7}); + param1.set_user_defined_value>("protocol", {1e-4, 1e-6, 1e-7}); param1.set_user_defined_value("dconv", 1e-6); } else { param1.set_user_defined_value>("protocol", {1e-8}); @@ -401,18 +382,15 @@ void runMOLDFT(World &world, const moldftSchema &moldftSchema, bool try_run, // if parameters are different or if I want to run no matter what run // if I want to restart and if I can. restart if (try_run) { - if (world.rank() == 0) { - print("-------------Running moldft------------"); - } + if (world.rank() == 0) { print("-------------Running moldft------------"); } // if params are different run and if restart exists and if im asking to restar if (std::filesystem::exists(moldftSchema.moldft_restart) && restart) { param1.set_user_defined_value("restart", true); } world.gop.fence(); if (world.rank() == 0) { - molresponse::write_test_input test_input( - param1, "moldft.in", - moldftSchema.mol_path);// molecule HF + molresponse::write_test_input test_input(param1, "moldft.in", + moldftSchema.mol_path);// molecule HF } world.gop.fence(); commandlineparser parser; @@ -451,8 +429,7 @@ void runMOLDFT(World &world, const moldftSchema &moldftSchema, bool try_run, // double energy=ME.value(calc.molecule.get_all_coords().flat()); // ugh! ME.value(calc.molecule.get_all_coords().flat());// ugh! world.gop.fence(); - const real_function_3d rho = - 2.0 * calc.make_density(world, calc.aocc, calc.amo); + const real_function_3d rho = 2.0 * calc.make_density(world, calc.aocc, calc.amo); auto dipole_t = calc.dipole(world, rho); std::map results; results["scf_energy"] = calc.current_energy; @@ -478,19 +455,16 @@ void runMOLDFT(World &world, const moldftSchema &moldftSchema, bool try_run, * @param xc * @param frequency */ -void set_excited_parameters(World &world, ResponseParameters &r_params, - const std::string &xc, const size_t &num_states, +void set_excited_parameters(World &world, ResponseParameters &r_params, const std::string &xc, const size_t &num_states, const std::string &precision) { if (world.rank() == 0) { if (precision == "high") { - r_params.set_user_defined_value>("protocol", - {1e-4, 1e-6, 1e-7}); + r_params.set_user_defined_value>("protocol", {1e-4, 1e-6, 1e-7}); r_params.set_user_defined_value("dconv", 1e-6); } else if (precision == "low") { - r_params.set_user_defined_value>("protocol", - {1e-4, 1e-6}); + r_params.set_user_defined_value>("protocol", {1e-4, 1e-6}); r_params.set_user_defined_value("dconv", 1e-4); } else { r_params.set_user_defined_value>("protocol", {1e-9}); @@ -522,20 +496,43 @@ void set_excited_parameters(World &world, ResponseParameters &r_params, * @param xc * @param frequency */ -void set_frequency_response_parameters(World &world, - ResponseParameters &r_params, - const std::string &property, - const std::string &xc, - const double &frequency, - const std::string &precision) { +void set_hyperpolarizability_parameters(World &world, ResponseParameters &r_params, const std::string &property, + const std::string &xc, const std::vector &frequency, + const std::string precision) { + if (world.rank() == 0) { + + r_params.set_user_defined_value("quadratic", true); + r_params.set_user_defined_value("freq_range", frequency); + r_params.set_user_defined_value("xc", xc); + // r_params.set_user_defined_value("property", property); + + if (precision == "high") { + r_params.set_user_defined_value>("protocol", {1e-7}); + r_params.set_user_defined_value("dconv", 1e-6); + } else if (precision == "low") { + r_params.set_user_defined_value>("protocol", {1e-6}); + r_params.set_user_defined_value("dconv", 1e-4); + } + } + world.gop.broadcast_serializable(r_params, 0); +} +/** + * Sets the response parameters for a frequency response calculation and writes + * to file + * + * @param r_params + * @param property + * @param xc + * @param frequency + */ +void set_frequency_response_parameters(World &world, ResponseParameters &r_params, const std::string &property, + const std::string &xc, const double &frequency, const std::string &precision) { if (world.rank() == 0) { if (precision == "high") { - r_params.set_user_defined_value>("protocol", - {1e-4, 1e-6, 1e-7}); + r_params.set_user_defined_value>("protocol", {1e-4, 1e-6, 1e-7}); r_params.set_user_defined_value("dconv", 1e-6); } else if (precision == "low") { - r_params.set_user_defined_value>("protocol", - {1e-4, 1e-6}); + r_params.set_user_defined_value>("protocol", {1e-4, 1e-6}); r_params.set_user_defined_value("dconv", 1e-4); } else { r_params.set_user_defined_value>("protocol", {1e-9}); @@ -571,37 +568,27 @@ void set_frequency_response_parameters(World &world, * @param frequency * @param moldft_path */ -static auto -set_frequency_path_and_restart(World &world, ResponseParameters ¶meters, - const std::string &property, - const double &frequency, const std::string &xc, - const std::filesystem::path &moldft_path, - std::filesystem::path &restart_path, - bool restart) -> std::filesystem::path { +static auto set_frequency_path_and_restart(World &world, ResponseParameters ¶meters, const std::string &property, + const double &frequency, const std::string &xc, + const std::filesystem::path &moldft_path, + std::filesystem::path &restart_path, bool restart) -> std::filesystem::path { if (world.rank() == 0) { print("restart path", restart_path); } // change the logic create save path - auto frequency_run_path = generate_response_frequency_run_path( - moldft_path, property, frequency, xc); + auto frequency_run_path = generate_response_frequency_run_path(moldft_path, property, frequency, xc); world.gop.fence(); if (world.rank() == 0) { print("frequency run path", frequency_run_path); } - // Crea if (std::filesystem::is_directory(frequency_run_path)) { - if (world.rank() == 0) { - cout << "Response directory found " << std::endl; - } + if (world.rank() == 0) { cout << "Response directory found " << std::endl; } } else {// create the file std::filesystem::create_directory(frequency_run_path); - if (world.rank() == 0) { - cout << "Creating response_path directory" << std::endl; - } + if (world.rank() == 0) { cout << "Creating response_path directory" << std::endl; } } std::filesystem::current_path(frequency_run_path); // frequency save path - auto [save_path, save_string] = - generate_frequency_save_path(frequency_run_path); + auto [save_path, save_string] = generate_frequency_save_path(frequency_run_path); if (world.rank() == 0) { parameters.set_user_defined_value("save", true); parameters.set_user_defined_value("save_file", save_string); @@ -613,14 +600,10 @@ set_frequency_path_and_restart(World &world, ResponseParameters ¶meters, parameters.set_user_defined_value("restart_file", save_string); } else if (std::filesystem::exists(restart_path)) { parameters.set_user_defined_value("restart", true); - auto split_restart_path = - split(restart_path.replace_extension("").string(), '/'); - std::string restart_file_short = - "../" + - split_restart_path[split_restart_path.size() - 2] + - "/" + split_restart_path[split_restart_path.size() - 1]; - parameters.set_user_defined_value("restart_file", - restart_file_short); + auto split_restart_path = split(restart_path.replace_extension("").string(), '/'); + std::string restart_file_short = "../" + split_restart_path[split_restart_path.size() - 2] + "/" + + split_restart_path[split_restart_path.size() - 1]; + parameters.set_user_defined_value("restart_file", restart_file_short); // Then we restart from the previous file instead } else { parameters.set_user_defined_value("restart", false); @@ -644,29 +627,21 @@ set_frequency_path_and_restart(World &world, ResponseParameters ¶meters, * @param restart_path * @return */ -auto RunResponse(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 { +auto RunResponse(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); + 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); - } + 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 if (std::filesystem::exists("response_base.json")) { std::ifstream ifs("response_base.json"); json response_base; ifs >> response_base; - if (response_base["converged"] && - response_base["precision"]["dconv"] == r_params.dconv()) { + if (response_base["converged"] && response_base["precision"]["dconv"] == r_params.dconv()) { return {save_path, true}; } } @@ -711,9 +686,8 @@ auto RunResponse(World &world, const std::string &filename, double frequency, * @param frequency * @param moldft_path */ -static void -set_and_write_restart_excited_parameters(ResponseParameters ¶meters, - excitedSchema &schema, bool restart) { +static void set_and_write_restart_excited_parameters(ResponseParameters ¶meters, excitedSchema &schema, + bool restart) { parameters.set_user_defined_value("save", true); parameters.set_user_defined_value("save_file", schema.save_string); @@ -759,16 +733,13 @@ static void create_excited_paths(excitedSchema &schema) { * @param restart_path * @return */ -auto runExcited(World &world, excitedSchema schema, bool restart, - const std::string& precision) -> bool { +auto runExcited(World &world, excitedSchema schema, bool restart, const std::string &precision) -> bool { // Set the response parameters ResponseParameters r_params{}; - set_excited_parameters(world, r_params, schema.xc, - schema.num_states, - precision); + set_excited_parameters(world, r_params, schema.xc, schema.num_states, precision); create_excited_paths(schema); std::filesystem::current_path(schema.excited_state_run_path); set_and_write_restart_excited_parameters(r_params, schema, restart); @@ -804,13 +775,11 @@ auto runExcited(World &world, excitedSchema schema, bool restart, * @param xc * @param property */ -void runFrequencyTests(World &world, const frequencySchema &schema, - const std::string &high_prec) { +void runFrequencyTests(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"); + 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) { @@ -825,12 +794,167 @@ void runFrequencyTests(World &world, const frequencySchema &schema, } else { throw Response_Convergence_Error{}; } - success = RunResponse(world, "response.in", freq, schema.op, schema.xc, - schema.moldft_path, restart_path, high_prec); + success = RunResponse(world, "response.in", freq, schema.op, schema.xc, schema.moldft_path, restart_path, + high_prec); if (world.rank() == 0) { print("Frequency ", freq, " completed"); } } } +/** + * Takes in a series of frequencies and runs a quadratic response calculations + * for given property at given frequencies. + * + * The frequencies are given in a vector of doubles + * For quadratic response the set of frequencies that a the first order response calculations have been run. + * + * The assumption is that response vectors are store in restart path. If they aren't already then we need to + * run the first order responnse calculation with the mad-freq executable. This would be equivalent to running + * the runFrequencyTests function. + * + * So the first step of this function is to check if the restart path exists for all frequencies. + * If it does, then we can run the quadratic response calculation. which does a double for loop through the frequencies + * + * Also, this function will save the quadratic response output to a json file formatted + * + * A-freq, B-freq, C-freq, A, B, C, , Beta(value) + * + * where A-freq, B-freq, C-freq are the frequencies used in the first order response calculations + * and A, B, C are the perturbation operator used in the quadratic response calculation + * and Beta is the value of the quadratic response calculation + * + * As an example, if A,B,C are dipole operators then this calculations gives the hyperpolarizability + * + * + * + * + * @param world + * @param moldft_path + * @param frequencies + * @param xc + * @param property + */ +void runQuadraticResponse(World &world, const frequencySchema &schema, const std::string precision) { + std::filesystem::current_path(schema.moldft_path); + + bool run_first_order = false; + + // add a restart path + auto restart_path = + addPath(schema.moldft_path, "/" + schema.op + "_0-000000.00000/restart_" + schema.op + "_0-000000.00000"); + try { + for (const auto &freq: schema.freq) { + std::filesystem::current_path(schema.moldft_path); + ResponseParameters r_params{}; + auto save_path = set_frequency_path_and_restart(world, r_params, schema.op, freq, schema.xc, + schema.moldft_path, restart_path, true); + print(save_path); + + if (std::filesystem::exists("response_base.json")) { + std::ifstream ifs("response_base.json"); + json response_base; + ifs >> response_base; + if (response_base["converged"] && response_base["precision"]["dconv"] == r_params.dconv()) { + if (world.rank() == 0) { print("Response calculation already converged"); } + continue; + } else { + if (world.rank() == 0) { print("Response calculation not converged"); } + run_first_order = true; + break; + } + } + if (!std::filesystem::exists(save_path)) { throw Response_Convergence_Error{}; } + } + // if we get here then all the first order response calculations have been run + // so we can run the quadratic response calculations + + + if (world.rank() == 0) { print("Running quadratic response calculations"); } + std::filesystem::current_path(schema.moldft_path); + RHS_Generator rhs_generator; + if (schema.op == "dipole") { + rhs_generator = dipole_generator; + } else { + rhs_generator = nuclear_generator; + } + ResponseParameters quad_parameters{}; + + set_hyperpolarizability_parameters(world, quad_parameters, schema.op, schema.xc, schema.freq, std::string()); + if (world.rank() == 0) { molresponse::write_response_input(quad_parameters, "quad.in"); } + + //auto calc_params = initialize_calc_params(world, std::string("quad.in")); + commandlineparser parser; + std::string moldft_archive = "moldft.restartdata"; + GroundStateCalculation ground_calculation{world, moldft_archive}; + if (world.rank() == 0) { ground_calculation.print_params(); } + Molecule molecule = ground_calculation.molecule(); + quad_parameters.set_ground_state_calculation_data(ground_calculation); + quad_parameters.set_derived_values(world, molecule); + if (world.rank() == 0) { quad_parameters.print(); } + + world.gop.fence(); + FunctionDefaults<3>::set_pmap(pmapT(new LevelPmap>(world))); + for (const auto &omega_b: schema.freq) { + for (const auto &omega_c: schema.freq) { + + auto generate_omega_restart_path = [&](double frequency) { + auto linear_response_calc_path = + generate_response_frequency_run_path(schema.moldft_path, schema.op, frequency, schema.xc); + auto restart_file_and_path = generate_frequency_save_path(linear_response_calc_path); + return restart_file_and_path.first; + }; + + + auto omega_a = omega_b + omega_c; + if (world.rank() == 0) { + print("New combination of frequencies ", omega_a, " ", omega_b, " ", omega_c); + print("-------------------------------------------"); + } + if (omega_a <= schema.freq.back()) { + + auto restartA = generate_omega_restart_path(omega_a); + auto restartB = generate_omega_restart_path(omega_b); + auto restartC = generate_omega_restart_path(omega_c); + if (world.rank() == 0) { + print("Restart file for A", restartA); + print("Restart file for B", restartB); + print("Restart file for C", restartC); + } + + // check if restartA exists + if (!std::filesystem::exists(restartA)) { + if (world.rank() == 0) { print("Restart file for omega_a = ", omega_a, " doesn't exist"); } + throw Response_Convergence_Error{}; + } else { + if (world.rank() == 0) { print("Restart file for omega_a = ", omega_a, " exists"); } + + + } + + + } else { + if (world.rank() == 0) { + print("Skipping omega_a = ", omega_a, " because it is greater than the max frequency"); + } + continue; + } + } + } + + + } catch (Response_Convergence_Error &e) { + if (true) { + // if the first order response calculations haven't been run then run them + if (world.rank() == 0) { print("Running first order response calculations"); } + + runFrequencyTests(world, schema, precision); + } else { + if (world.rank() == 0) { + print("First order response calculations haven't been run and can't be run"); + print("Quadratic response calculations can't be run"); + } + } + } +} /** * @@ -840,25 +964,19 @@ void runFrequencyTests(World &world, const frequencySchema &schema, * @param restart do we force a restart or not * @param precision high precision or no? */ -void moldft(World &world, moldftSchema &m_schema, bool try_moldft, bool restart, - const std::string &precision) { +void moldft(World &world, moldftSchema &m_schema, bool try_moldft, bool restart, const std::string &precision) { if (std::filesystem::is_directory(m_schema.moldft_path)) { - if (world.rank() == 0) { - cout << "MOLDFT directory found " << m_schema.mol_path << "\n"; - } + if (world.rank() == 0) { cout << "MOLDFT directory found " << m_schema.mol_path << "\n"; } } else {// create the file if (world.rank() == 0) { std::filesystem::create_directory(m_schema.moldft_path); - cout << "Creating MOLDFT directory for " << m_schema.mol_name - << ":/" << m_schema.moldft_path << ":\n"; + cout << "Creating MOLDFT directory for " << m_schema.mol_name << ":/" << m_schema.moldft_path << ":\n"; } world.gop.fence(); } std::filesystem::current_path(m_schema.moldft_path); - if (world.rank() == 0) { - cout << "Entering : " << m_schema.moldft_path << " to run MOLDFT \n\n"; - } + if (world.rank() == 0) { cout << "Entering : " << m_schema.moldft_path << " to run MOLDFT \n\n"; } runMOLDFT(world, m_schema, try_moldft, restart, precision); }