Skip to content

Commit

Permalink
add first attempt at LZSH
Browse files Browse the repository at this point in the history
  • Loading branch information
@tjira
tjira committed Apr 26, 2024
1 parent 1e9db1f commit e708027
Show file tree
Hide file tree
Showing 10 changed files with 150 additions and 103 deletions.
9 changes: 3 additions & 6 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,19 +22,16 @@
<a href="https://github.com/tjira/acorn/stargazers">
<img src="https://img.shields.io/github/stars/tjira/acorn?style=for-the-badge"/>
</a>
<a href="https://github.com/tjira/acorn/releases/latest">
<img src="https://img.shields.io/github/downloads/tjira/acorn/total?style=for-the-badge"/>
</a>
<br>
<a href="https://github.com/tjira/acorn">
<img src="https://tokei.rs/b1/github/tjira/acorn?category=code&style=for-the-badge"/>
</a>
<a href="https://github.com/tjira/acorn">
<img src="https://img.shields.io/github/languages/code-size/tjira/acorn?style=for-the-badge"/>
</a>
<a href="https://github.com/tjira/acorn/releases/latest">
<img src="https://img.shields.io/github/v/release/tjira/acorn?display_name=tag&style=for-the-badge"/>
</a>
<a href="https://github.com/tjira/acorn/releases/latest">
<img src="https://img.shields.io/github/downloads/tjira/acorn/total?style=for-the-badge"/>
</a>
</p>

<p align="center">
Expand Down
2 changes: 1 addition & 1 deletion docs/pages/configurationinteraction.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ where $\otimes_K$ stands for the Kronecker product. Similarly, the transformatio
J_{pqrs}^{MS}=C_{\mu p}C_{\nu q}(\mathbf{I}\_{2}\otimes_K(\mathbf{I}\_{2}\otimes_K\mathbf{J})^{(4,3,2,1)})\_{\mu\nu\kappa\lambda}C_{\kappa r}C_{\lambda s}
\end{equation}

This notation accounts for the spin modifications and ensures that the transformations adhere to quantum mechanical principles.
where the superscript $(4,3,2,1)$ denotes the axes transposition. This notation accounts for the spin modifications and ensures that the transformations adhere to quantum mechanical principles.

### The Determinants Generation

Expand Down
16 changes: 8 additions & 8 deletions example/input/ad1d.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,13 @@
"ngrid" : 1024,
"mass" : 1
},
"solve" : {
"dynamics" : {
"gradient" : ["x"],
"iters" : 1000,
"step" : 0.05,
"position" : [-5],
"state" : 1
}
"dynamics" : {
"gradient" : ["x"],
"iters" : 1000,
"step" : 0.05,
"position" : [-5],
"velocity" : [0],
"state" : 1,
"seed" : 1
}
}
18 changes: 9 additions & 9 deletions example/input/nad1d.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,14 @@
"ngrid" : 4096,
"mass" : 2000
},
"solve" : {
"dynamics" : {
"gradient" : ["0.006/cosh(0.6*x)^2", "-0.006/cosh(0.6*x)^2"],
"iters" : 1000,
"step" : 10,
"position" : [-10],
"velocity" : [0.005475],
"state" : 1
}
"dynamics" : {
"gradient" : ["0.006/cosh(0.6*x)^2", "-0.006/cosh(0.6*x)^2"],
"iters" : 400,
"step" : 10,
"position" : [-10],
"velocity" : [0.005475],
"state" : 1,
"trajs" : 100,
"seed" : 1
}
}
1 change: 0 additions & 1 deletion example/input/qnd1d_real.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
"cap" : "0.1*exp(-0.1*(x+24)^2)+0.1*exp(-0.1*(x-24)^2)",
"momentum" : 10.95,
"iters" : 350,
"real" : true,
"step" : 10,
"savewfn" : true,
"adiabatic" : true
Expand Down
14 changes: 7 additions & 7 deletions include/default.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,15 +21,15 @@ inline std::string moloptstr = R"({
})";

inline std::string msaoptstr = R"({
"real" : false, "step" : 0.1, "guess" : "-x^2", "nstate" : 1, "iters" : 1000, "savewfn" : false,
"optimize" : {"step" : 1, "iters" : 1000},
"spectrum" : {"normalize" : false, "zpesub" : false, "potential" : "", "window" : "1", "zeropad" : 0},
"dynamics" : {"iters" : 100, "step" : 20, "state" : 1, "position" : [0], "velocity" : [0], "gradient" : []}
"spectrum" : {"normalize" : false, "zpesub" : false, "potential" : "", "window" : "1", "zeropad" : 0}
})";

inline std::string msdoptstr = R"({
"trajs" : 1
})";

inline std::string msnoptstr = R"({
"step" : 0.1, "guess" : ["-x^2"], "iters" : 1000, "savewfn" : false, "cap" : "0", "momentum" : 0, "adiabatic" : true,
"dynamics" : {"iters" : 100, "step" : 20, "state" : 1, "position" : [0], "velocity" : [0], "gradient" : []}
})";

inline std::string rmpoptstr = R"({
Expand All @@ -53,9 +53,9 @@ inline std::string uhfoptstr = R"({
})";

inline std::string bgloptstr = R"({
"dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
"dynamics" : {"berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
})";

inline std::string orcoptstr = R"({
"dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
"dynamics" : {"berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
})";
23 changes: 9 additions & 14 deletions include/modelsolver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,6 @@ class ModelSolver {
public:
struct OptionsAdiabatic {OptionsAdiabatic(){};
// options structures
struct Dynamics {
int iters=100, state=0; double step=1.0;
std::vector<double> position, velocity;
std::vector<std::string> gradient;
} dynamics={};
struct Optimize {
double step=1; int iters=1000;
} optimize={};
Expand All @@ -26,30 +21,30 @@ class ModelSolver {
std::string guess;
};
struct OptionsNonadiabatic {OptionsNonadiabatic(){};
// dynamics structure
struct Dynamics {
int iters=100, state=0; double step=1.0;
std::vector<double> position, velocity;
std::vector<std::string> gradient;
} dynamics={};

// variables
int iters=1000; std::vector<std::string> guess; std::string cap="0";
bool savewfn=false, adiabatic=true; double step=0.1, momentum=0;
};
struct OptionsDynamics {
int iters=100, seed=1, state=0, trajs=1; double step=1.0;
std::vector<double> position, velocity;
std::vector<std::string> gradient;
} dynamics={};

public:
// constructors and destructors
ModelSolver(const OptionsNonadiabatic& optn) : optn(optn) {};
ModelSolver(const OptionsAdiabatic& opta) : opta(opta) {};
ModelSolver(const OptionsDynamics& optd) : optd(optd) {};

// methods and solvers
Result run(const ModelSystem& system, Result res = {}, bool print = true);

private:
template <typename T> Result runcd(const ModelSystem& system, const T& optdyn, Result res = {}, bool print = true);
Result runnad(const ModelSystem& system, Result res = {}, bool print = true);
Result runad(const ModelSystem& system, Result res = {}, bool print = true);
Result runcd(const ModelSystem& system, Result res = {}, bool print = true);

private:
OptionsAdiabatic opta; OptionsNonadiabatic optn;
OptionsAdiabatic opta; OptionsNonadiabatic optn; OptionsDynamics optd;
};
16 changes: 8 additions & 8 deletions src/determinant.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,12 @@ std::vector<Determinant> Determinant::excitations(const std::vector<int>& excs)

// single excitations
if (std::find(excs.begin(), excs.end(), 1) != excs.end()) {
for (int i = 0; i < a.size(); i++) {
for (size_t i = 0; i < a.size(); i++) {
for (int j = a.size(); j < norb; j++) {
std::vector<int> det = this->a; det.at(i) = j; excitations.push_back(Determinant(norb, det, b));
}
}
for (int i = 0; i < b.size(); i++) {
for (size_t i = 0; i < b.size(); i++) {
for (int j = b.size(); j < norb; j++) {
std::vector<int> det = this->b; det.at(i) = j; excitations.push_back(Determinant(norb, a, det));
}
Expand All @@ -61,27 +61,27 @@ std::vector<Determinant> Determinant::excitations(const std::vector<int>& excs)

// double excitations
if (std::find(excs.begin(), excs.end(), 2) != excs.end()) {
for (int i = 0; i < a.size(); i++) {
for (size_t i = 0; i < a.size(); i++) {
for (int j = a.size(); j < norb; j++) {
for (int k = 0; k < b.size(); k++) {
for (size_t k = 0; k < b.size(); k++) {
for (int l = b.size(); l < norb; l++) {
std::vector<int> deta = this->a, detb = this->b; deta.at(i) = j, detb.at(k) = l; excitations.push_back(Determinant(norb, deta, detb));
}
}
}
}
for (int i = 0; i < a.size(); i++) {
for (size_t i = 0; i < a.size(); i++) {
for (int j = a.size(); j < norb; j++) {
for (int k = i + 1; k < a.size(); k++) {
for (size_t k = i + 1; k < a.size(); k++) {
for (int l = j + 1; l < norb; l++) {
std::vector<int> det = this->a; det.at(i) = j, det.at(k) = l; excitations.push_back(Determinant(norb, det, b));
}
}
}
}
for (int i = 0; i < b.size(); i++) {
for (size_t i = 0; i < b.size(); i++) {
for (int j = b.size(); j < norb; j++) {
for (int k = i + 1; k < b.size(); k++) {
for (size_t k = i + 1; k < b.size(); k++) {
for (int l = j + 1; l < norb; l++) {
std::vector<int> det = this->b; det.at(i) = j, det.at(k) = l; excitations.push_back(Determinant(norb, a, det));
}
Expand Down
37 changes: 27 additions & 10 deletions src/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,13 +56,12 @@ NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options::Dynamics, iters, step, berends
// option structures loaders for model methods
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Optimize, step, iters);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Spectrum, potential, window, normalize, zpesub, zeropad);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic::Dynamics, iters, step, state, position, gradient, velocity);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Dynamics, iters, step, state, position, gradient, velocity);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsDynamics, iters, step, state, position, gradient, velocity, seed, trajs);

// option loaders
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic, dynamics, real, step, iters, nstate, optimize, guess, savewfn, spectrum);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options, dynamics, gradient, hessian, excitations, nstate, state);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic, dynamics, step, iters, guess, savewfn, cap, momentum, adiabatic);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic, real, step, iters, nstate, optimize, guess, savewfn, spectrum);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic, step, iters, guess, savewfn, cap, momentum, adiabatic);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(UnrestrictedHartreeFock::Options, dynamics, gradient, hessian, maxiter, thresh);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options, dynamics, gradient, hessian, maxiter, thresh);
NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options, dynamics, gradient, hessian, order);
Expand Down Expand Up @@ -114,6 +113,7 @@ int main(int argc, char** argv) {
auto mdlopt = nlohmann::json::parse(mdloptstr);
auto molopt = nlohmann::json::parse(moloptstr);
auto msaopt = nlohmann::json::parse(msaoptstr);
auto msdopt = nlohmann::json::parse(msdoptstr);
auto msnopt = nlohmann::json::parse(msnoptstr);
auto orcopt = nlohmann::json::parse(orcoptstr);
auto bglopt = nlohmann::json::parse(bgloptstr);
Expand All @@ -125,6 +125,7 @@ int main(int argc, char** argv) {
// patch the input json file and apply defaults
if (input.contains("integral")) intopt.merge_patch(input.at("integral")), fullinput["integral"] = intopt;
if (input.contains("molecule")) molopt.merge_patch(input.at("molecule")), fullinput["molecule"] = molopt;
if (input.contains("dynamics")) msdopt.merge_patch(input.at("dynamics")), fullinput["dynamics"] = msdopt;
if (input.contains("bagel")) bglopt.merge_patch(input.at("bagel")), fullinput["bagel"] = bglopt;
if (input.contains("model")) mdlopt.merge_patch(input.at("model")), fullinput["model"] = mdlopt;
if (input.contains("orca")) orcopt.merge_patch(input.at("orca")), fullinput["orca"] = orcopt;
Expand Down Expand Up @@ -424,11 +425,8 @@ int main(int argc, char** argv) {
ModelSystem model(mdlopt.at("mass"), mdlopt.at("potential"), mdlopt.at("variables"), mdlopt.at("limits"), mdlopt.at("ngrid"));

if (input.contains("solve")) {Printer::Title("EXACT QUANTUM DYNAMICS");
// create the adiabatic solver object
ModelSolver msv(msaopt.get<ModelSolver::OptionsAdiabatic>());

// if nonadiabatic model was specified assign to the solver the nonadiabatic version
if (mdlopt.at("potential").size() > 1) msv = ModelSolver(msnopt.get<ModelSolver::OptionsNonadiabatic>());
ModelSolver msv = mdlopt.at("potential").size() == 1 ? ModelSolver(msaopt.get<ModelSolver::OptionsAdiabatic>()) : ModelSolver(msnopt.get<ModelSolver::OptionsNonadiabatic>());

// run the optimization if requested
if (mdlopt.at("potential").size() == 1 && input.at("solve").contains("optimize")) {
Expand Down Expand Up @@ -462,10 +460,29 @@ int main(int argc, char** argv) {
// save the spectral analysis results if available
if (res.msv.spectra.size()) {
for (size_t i = 0; i < res.msv.spectra.size(); i++) {
Matrix<> A(res.msv.t.size(), 3); A << res.msv.t, res.msv.acfs.at(i).real(), res.msv.acfs.at(i).imag(); EigenWrite(std::filesystem::path(ip) / ("acf" + std::to_string(i) + ".mat"), A);
Matrix<> F(res.msv.f.size(), 2); F << res.msv.f, res.msv.spectra.at(i).cwiseAbs(); EigenWrite(std::filesystem::path(ip) / ("spectrum" + std::to_string(i) + ".mat"), F);
Matrix<> A(res.msv.t.size(), 3); A << res.msv.t, res.msv.acfs.at(i).real(), res.msv.acfs.at(i).imag(); EigenWrite(std::filesystem::path(ip) / ("acf" + std::to_string(i + 1) + ".mat"), A);
Matrix<> F(res.msv.f.size(), 2); F << res.msv.f, res.msv.spectra.at(i).cwiseAbs(); EigenWrite(std::filesystem::path(ip) / ("spectrum" + std::to_string(i + 1) + ".mat"), F);
}
}
} else if (input.contains("dynamics")) {Printer::Title("MODEL CLASSICAL DYNAMICS");
// create the classical solver object
ModelSolver msv(msdopt.get<ModelSolver::OptionsDynamics>());

// run the calculation
res = msv.run(model, res);

// extract and save the potential points
if (model.vars().size() == 2) {
Matrix<> U(res.msv.r.size() * res.msv.r.size(), res.msv.U.cols() + 2);
for (int i = 0; i < res.msv.r.size(); i++) {
for (int j = 0; j < res.msv.r.size(); j++) {
U(i * res.msv.r.size() + j, 0) = res.msv.r(i), U(i * res.msv.r.size() + j, 1) = res.msv.r(j);
for (int k = 0; k < res.msv.U.cols(); k++) {
U(i * res.msv.r.size() + j, k + 2) = res.msv.U(i * res.msv.r.size() + j, k);
}
}
} EigenWrite(std::filesystem::path(ip) / "U.mat", U);
} else {Matrix<> U(res.msv.r.size(), res.msv.U.cols() + 1); U << res.msv.r, res.msv.U; EigenWrite(std::filesystem::path(ip) / "U.mat", U);}
}
}

Expand Down
Loading

0 comments on commit e708027

Please sign in to comment.