Skip to content

Commit

Permalink
Version 1 of MADQC using a createCalcManager as a factory
Browse files Browse the repository at this point in the history
  • Loading branch information
ahurta92 committed Aug 28, 2024
1 parent 02ff1cf commit 50cb053
Show file tree
Hide file tree
Showing 5 changed files with 467 additions and 272 deletions.
171 changes: 171 additions & 0 deletions src/apps/madqc/calc_factory.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
#ifndef CALC_STRATS_HPP
#define CALC_STRATS_HPP

#include "calc_manager.hpp"
#include "parameter_manager.hpp"

// The driver is the calculation type (e.g. energy, gradient, hessian, etc.)
// The model is the combination of theory and basis set
// For us, the basis set is adaptive, so we don't need to worry about that
// we simply compute to a certain precision

// Let's start of defining the energy driver. I need to create a calc_manager
// The models are going to reflect which application we are running with
// - Moldft
// - Response (moldft + molresponse)
// - MP2 (moldft + mp2)
// - CIS (moldft + cis)
// - OEP (moldft + oep)
// - MP3 (moldft + mp3)
//
//
//
//
//

// define the properties that can be calculated for each model
// This example is for the response model
struct response_property_map {
bool alpha;
bool beta;
bool shg;
};

void to_json(nlohmann::json& j, const response_property_map& r) {
j = nlohmann::json{{"alpha", r.alpha}, {"beta", r.beta}, {"shg", r.shg}};
}
void from_json(const nlohmann::json& j, response_property_map& r) {
j.at("alpha").get_to(r.alpha);
j.at("beta").get_to(r.beta);
j.at("shg").get_to(r.shg);
}

// A property map is a map of properties to be calculated for each model
using model = std::string;
using model_properties = std::map<std::string, bool>;
//maybe just use json for this?
using property_map = std::map<model, model_properties>;

std::unique_ptr<CalcManager> createCalcManager(const std::string& model_name, ParameterManager& pm,
property_map properties) {
// Create a new CalcManager
auto calc_manager = std::make_unique<CalcManager>();
auto molecule = pm.get_molecule();
// All calculations start with a reference
auto params = pm.get_moldft_params();
auto moldft_properties = properties["moldft"];
auto moldft = std::make_unique<MoldftCalculationStrategy>(params, molecule, "moldft", moldft_properties);
calc_manager->addStrategy(std::move(moldft));

if (model_name == "moldft") {

return calc_manager;

} else if (model_name == "response") {

auto& response_params = pm.get_molresponse_params();
auto response_properties = properties.at("response"); // map of properties to be calculated for response

auto perturbations = response_params.perturbations();
// say you want to do moldft hf and response lda... well it's possible here TODO: vecotrize this so you can do mixed?
auto xc = response_params.xc();
auto freq_range = response_params.freq_range();

vector<std::string> input_names = {"moldft"};
for (auto const& perturbation : perturbations) {
ResponseInput r_input = std::make_tuple(perturbation, xc, freq_range);
auto response_calc =
std::make_unique<LinearResponseStrategy>(response_params, r_input, perturbation, input_names);
calc_manager->addStrategy(std::move(response_calc));
}

if (response_properties["beta"]) {

auto set_freqs = [&]() {
vector<double> freqs_copy = freq_range;
auto num_freqs = freq_range.size();
auto compare_freqs = [](double x, double y) {
return std::abs(x - y) < 1e-3;
};

for (int i = 0; i < num_freqs; i++) { // for i=0:n-1
for (int j = i; j < num_freqs; j++) { // for j = i omega_3=-(omega_1+omega_2)
auto omega_1 = freq_range[i];
auto omega_2 = freq_range[j];
auto omega_3 = omega_1 + omega_2;

// if you can find omega_3 in freq_copy skip
if (omega_2 == 0.0)
continue;
if (std::find_if(freqs_copy.begin(), freqs_copy.end(),
[&](double x) { return compare_freqs(x, omega_3); }) != freqs_copy.end()) {
continue;
} else {
freqs_copy.push_back(omega_3);
}
}
}
return freqs_copy;
};

freq_range = set_freqs();
print("Frequency Range: ", freq_range);
// this is where I we create our calculation
ResponseInput hyp_input = std::make_tuple("dipole", xc, freq_range);
std::vector<std::string> r_input_names = {"moldft"};
std::vector<std::string> h_input_names = {"moldft", "dipole"};
auto response_hyper =
std::make_unique<LinearResponseStrategy>(response_params, hyp_input, "dipole", r_input_names);
auto hyper_calc = std::make_unique<ResponseHyper>(response_params, hyp_input, "hyper", h_input_names);
calc_manager->addStrategy(std::move(response_hyper));
calc_manager->addStrategy(std::move(hyper_calc));
}
if (response_properties["shg"]) {
throw std::invalid_argument("SHG not implemented yet");
}
} else if (model_name == "MP2") {
throw std::invalid_argument("MP2 not implemented yet");
/*auto moldft_strategy = std::make_unique<MoldftCalculationStrategy>(params, molecule);*/
/*auto mp2_strategy = std::make_unique<MP2CalculationStrategy>(params, molecule);*/
/*calc_manager->addStrategy(std::move(moldft_strategy));*/
/*calc_manager->addStrategy(std::move(mp2_strategy));*/
} else if (model_name == "CIS") {
throw std::invalid_argument("CIS not implemented yet");
/*auto moldft_strategy = std::make_unique<MoldftCalculationStrategy>(params, molecule);*/
/*auto cis_strategy = std::make_unique<CISCalculationStrategy>(params, molecule);*/
/*calc_manager->addStrategy(std::move(moldft_strategy));*/
/*calc_manager->addStrategy(std::move(cis_strategy));*/
} else if (model_name == "OEP") {
throw std::invalid_argument("OEP not implemented yet");
/*auto moldft_strategy = std::make_unique<MoldftCalculationStrategy>(params, molecule);*/
/*auto oep_strategy = std::make_unique<OEPCalculationStrategy>(params, molecule);*/
/*calc_manager->addStrategy(std::move(moldft_strategy));*/
/*calc_manager->addStrategy(std::move(oep_strategy));*/
} else if (model_name == "MP3") {
throw std::invalid_argument("MP3 not implemented yet");
/*auto moldft_strategy = std::make_unique<MoldftCalculationStrategy>(params, molecule);*/
/*auto mp3_strategy = std::make_unique<MP3CalculationStrategy>(params, molecule);*/
/*calc_manager->addStrategy(std::move(moldft_strategy));*/
/*calc_manager->addStrategy(std::move(mp3_strategy));*/
} else {
throw std::invalid_argument("Unknown model name: " + model_name);
}
// Return the configured CalcManager
return calc_manager;
}

// Below here we define more complex calculation strategies
DynamicInput example_calculation(Molecule& molecule) {

// Example setup calculation function
auto setupCalculation = [&](CalcManager& calc_manager, const std::string& iter_name) {
// Create strategies and set them up for the current iteration
auto params = CalculationParameters();
auto moldft_strategy = std::make_unique<MoldftCalculationStrategy>(params, molecule);
calc_manager.addStrategy(std::move(moldft_strategy));
};

return DynamicInput{};
}

#endif
171 changes: 92 additions & 79 deletions src/apps/madqc/calc_manager.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,10 @@ class PathManager {
auto calc_path = calc_paths.get<CalculationTemplate>();
for (const auto& calc_dir : calc_path.calc_paths) {
print("Creating directory: ", calc_dir);
if (!std::filesystem::exists(calc_dir.parent_path())) {
std::filesystem::create_directory(calc_dir.parent_path());
}
// if base direcotry does not exist create it
if (!std::filesystem::exists(calc_dir)) {
std::filesystem::create_directory(calc_dir);
}
Expand Down Expand Up @@ -105,13 +109,13 @@ class CalculationStrategy {
virtual void compute(World& world, PathManager& path_manager) = 0;
virtual ~CalculationStrategy() = default;

CalculationStrategy(std::string calc_name, vector<std::string> properties)
CalculationStrategy(std::string calc_name, std::map<std::string, bool> properties)
: name(std::move(calc_name)), requested_properties(std::move(properties)){};
CalculationStrategy() = default;

protected:
std::string name;
std::vector<std::string> requested_properties;
std::map<std::string, bool> requested_properties;
};

class CompositeCalculationStrategy : public CalculationStrategy {
Expand Down Expand Up @@ -147,7 +151,7 @@ class MoldftCalculationStrategy : public CalculationStrategy {
public:
// Notice here that I am passing the parameters of the calculation as well as name
MoldftCalculationStrategy(const CalculationParameters& params, Molecule mol, std::string calc_name = "moldft",
std::vector<std::string> properties = {"energy"})
std::map<std::string, bool> properties = {{"energy", true}})
: parameters(params),
molecule(std::move(mol)),
CalculationStrategy(std::move(calc_name), std::move(properties)){};
Expand Down Expand Up @@ -278,8 +282,10 @@ class MoldftCalculationStrategy : public CalculationStrategy {
paths[name]["output"]["properties"] = {};
auto& output = paths[name]["output"]["properties"];

for (const auto& property : requested_properties) {
output[property] = compute_property(world, property, calc, energy);
for (const auto& [property, compute] : requested_properties) {
if (compute) {
output[property] = compute_property(world, property, calc, energy);
}
}
if (world.rank() == 0) {
print("output: ", output.dump(4));
Expand Down Expand Up @@ -514,7 +520,7 @@ class LinearResponseStrategy : public CalculationStrategy, InputInterface {
auto moldft_paths = path_manager.getPath(input_names[0]);
auto response_paths = path_manager.getPath(calc_name);

auto moldft_restart = moldft_paths.restarts[0].string();
auto moldft_restart = moldft_paths.restarts[0].replace_extension("").string();
auto alpha_outpath = response_paths.outputs["alpha"][0];

// I need to analyze alpha.json to see which frequencies are missing.
Expand Down Expand Up @@ -845,15 +851,6 @@ class ResponseHyper : public CalculationStrategy, InputInterface {
}
};

// A finite difference startegy will manage it's own calc_manager
class FiniteDifferenceStrategy : public CalculationStrategy, InputInterface {

std::string property;

FiniteDifferenceStrategy(const std::string& property, const std::vector<std::string>& inputs)
: property(property), InputInterface(inputs) {}
};

/**/
/*class WriteResponseVTKOutputStrategy : public CalculationStrategy {*/
/**/
Expand Down Expand Up @@ -994,70 +991,86 @@ class CalcManager {
using SetupCalculationFunction = std::function<void(CalcManager&, std::string)>;
using CheckConditionFunction = std::function<bool(World&, PathManager&, const std::string&)>;
using FinalizeCalculationFunction = std::function<void(World&, PathManager&, const path&)>;

using DynamicInput = std::tuple<SetupCalculationFunction, CheckConditionFunction, FinalizeCalculationFunction>;
//
//
class DynamicCalculationStrategy : public CalculationStrategy {
private:
std::string strategy_name;
int iteration_limit = 100; // Example limit to prevent infinite loops
std::vector<std::unique_ptr<CalculationStrategy>> strategies;
SetupCalculationFunction setup_calculation;
CheckConditionFunction check_condition;
FinalizeCalculationFunction finalize_calculation;

public:
// Constructor
explicit DynamicCalculationStrategy(PathManager& path_manager, std::string strategy_name,
std::vector<std::unique_ptr<CalculationStrategy>> strategies,
std::vector<std::string> properties = {})
: strategy_name(std::move(strategy_name)),
strategies(std::move(strategies)),
CalculationStrategy(std::move(strategy_name), std::move(properties)) {}

void setPaths(PathManager& path_manager, const path& root) override {

// Set the calc paths for the base of the calculation... In this case we don't know beforehand what we are working with.
// All we need to communicate is that we are working a directory specified by the strategy_name
CalculationTemplate dynamic_paths;
dynamic_paths.calc_paths.push_back(root / strategy_name);
path_manager.setPath(strategy_name, dynamic_paths);
}

void compute(World& world, PathManager& path_manager) override {
int iteration = 0;
bool condition_met = false;

auto base_path = path_manager.getPath(strategy_name).calc_paths[0]; // only one base path
// make sure the base path is created
path_manager.createDirectories(strategy_name);

vector<std::string> iteration_names;
while (!condition_met && iteration < iteration_limit) {
// Create subdirectory for this iteration
// now we create a new calculation manager for this iteration
CalcManager calc_manager(base_path); // the base path is the root for this iteration
auto iter_name = strategy_name + "_iter_" + std::to_string(iteration);
iteration_names.push_back(iter_name);

// Set up the calculation for the iteration
setup_calculation(calc_manager, iter_name);

// Run the calculations
calc_manager.runCalculations(world);

condition_met = check_condition(world, path_manager, iter_name);

// Check if the condition is met

if (!condition_met) {
// Prepare for the next iteration
iteration++;
} else {

// Finalize the calculation
finalize_calculation(world, path_manager, iter_name);
}
}
};
};
// This strategy a number of calculations based on a condition evaluated at runtime
//
/*class DynamicCalculationStrategy : public CalculationStrategy {*/
/* private:*/
/* std::string strategy_name;*/
/* int iteration_limit = 100; // Example limit to prevent infinite loops*/
/* SetupCalculationFunction setup_calculation;*/
/* CheckConditionFunction check_condition;*/
/* FinalizeCalculationFunction finalize_calculation;*/
/**/
/* public:*/
/* // Constructor*/
/* explicit DynamicCalculationStrategy(std::string strategy_name, SetupCalculationFunction setup_calculation,*/
/* CheckConditionFunction check_condition,*/
/* FinalizeCalculationFunction finalize_calculation,*/
/* std::vector<std::string> properties = {})*/
/* : strategy_name(std::move(strategy_name)),*/
/* setup_calculation(std::move(setup_calculation)),*/
/* check_condition(std::move(check_condition)),*/
/* finalize_calculation(std::move(finalize_calculation)),*/
/* CalculationStrategy(std::move(strategy_name), std::move(properties)) {}*/
/**/
/* explicit DynamicCalculationStrategy(std::string strategy_name, DynamicInput input,*/
/* std::vector<std::string> properties = {})*/
/* : strategy_name(std::move(strategy_name)),*/
/* setup_calculation(std::get<0>(input)),*/
/* check_condition(std::get<1>(input)),*/
/* finalize_calculation(std::get<2>(input)),*/
/* CalculationStrategy(std::move(strategy_name), std::move(properties)) {}*/
/**/
/* void setPaths(PathManager& path_manager, const path& root) override {*/
/* // Set the calc paths for the base of the calculation... In this case we don't know beforehand what we are working with.*/
/* // All we need to communicate is that we are working a directory specified by the strategy_name*/
/* CalculationTemplate dynamic_paths;*/
/* dynamic_paths.calc_paths.push_back(root / strategy_name);*/
/* path_manager.setPath(strategy_name, dynamic_paths);*/
/* }*/
/**/
/* void compute(World& world, PathManager& path_manager) override {*/
/* int iteration = 0;*/
/* bool condition_met = false;*/
/**/
/* auto base_path = path_manager.getPath(strategy_name).calc_paths[0]; // only one base path*/
/* path_manager.createDirectories(strategy_name);*/
/* vector<std::string> iteration_names;*/
/**/
/* // I need to be able to iterate over variable containers*/
/* // where the container can be either a set of numbers*/
/* // or a set of coordinates*/
/**/
/* for (int i = 0; i < iteration_limit; i++) {*/
/* // Create subdirectory for this iteration*/
/* // now we create a new calculation manager for this iteration*/
/* CalcManager calc_manager(base_path); // the base path is the root for this iteration*/
/* auto iter_name = strategy_name + "_iter_" + std::to_string(iteration);*/
/* iteration_names.push_back(iter_name);*/
/**/
/* // Set up the calculation for the iteration*/
/* setup_calculation(calc_manager, iter_name);*/
/**/
/* // Run the calculations*/
/* calc_manager.runCalculations(world);*/
/**/
/* condition_met = check_condition(world, path_manager, iter_name);*/
/**/
/* // Check if the condition is met*/
/**/
/* if (!condition_met) {*/
/* // Prepare for the next iteration*/
/* iteration++;*/
/* } else {*/
/**/
/* // Finalize the calculation*/
/* finalize_calculation(world, path_manager, iter_name);*/
/* }*/
/* }*/
/* }*/
/*};*/
#endif
Loading

0 comments on commit 50cb053

Please sign in to comment.