From e708027495075a2abdddfa893f32a2ecc04b69f3 Mon Sep 17 00:00:00 2001 From: Tomas Jira Date: Fri, 26 Apr 2024 15:30:53 +0200 Subject: [PATCH] add first attempt at LZSH --- README.md | 9 +- docs/pages/configurationinteraction.md | 2 +- example/input/ad1d.json | 16 ++-- example/input/nad1d.json | 18 ++-- example/input/qnd1d_real.json | 1 - include/default.h | 14 +-- include/modelsolver.h | 23 ++--- src/determinant.cpp | 16 ++-- src/main.cpp | 37 +++++--- src/modelsolver.cpp | 117 ++++++++++++++++--------- 10 files changed, 150 insertions(+), 103 deletions(-) diff --git a/README.md b/README.md index 03e18c1..8a1fa52 100644 --- a/README.md +++ b/README.md @@ -22,19 +22,16 @@ - - -
- - - + + +

diff --git a/docs/pages/configurationinteraction.md b/docs/pages/configurationinteraction.md index 1e769e1..616d6aa 100644 --- a/docs/pages/configurationinteraction.md +++ b/docs/pages/configurationinteraction.md @@ -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 diff --git a/example/input/ad1d.json b/example/input/ad1d.json index d0e500b..f666b35 100644 --- a/example/input/ad1d.json +++ b/example/input/ad1d.json @@ -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 } } diff --git a/example/input/nad1d.json b/example/input/nad1d.json index 74214d7..2162551 100644 --- a/example/input/nad1d.json +++ b/example/input/nad1d.json @@ -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 } } diff --git a/example/input/qnd1d_real.json b/example/input/qnd1d_real.json index 87a8664..c4c1d3c 100644 --- a/example/input/qnd1d_real.json +++ b/example/input/qnd1d_real.json @@ -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 diff --git a/include/default.h b/include/default.h index ee34175..a870ee2 100644 --- a/include/default.h +++ b/include/default.h @@ -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"({ @@ -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}} })"; diff --git a/include/modelsolver.h b/include/modelsolver.h index 29481ef..ac7bdd1 100644 --- a/include/modelsolver.h +++ b/include/modelsolver.h @@ -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 position, velocity; - std::vector gradient; - } dynamics={}; struct Optimize { double step=1; int iters=1000; } optimize={}; @@ -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 position, velocity; - std::vector gradient; - } dynamics={}; - // variables int iters=1000; std::vector 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 position, velocity; + std::vector 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 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; }; diff --git a/src/determinant.cpp b/src/determinant.cpp index 30710c4..4b5c0f3 100644 --- a/src/determinant.cpp +++ b/src/determinant.cpp @@ -47,12 +47,12 @@ std::vector Determinant::excitations(const std::vector& 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 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 det = this->b; det.at(i) = j; excitations.push_back(Determinant(norb, a, det)); } @@ -61,27 +61,27 @@ std::vector Determinant::excitations(const std::vector& 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 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 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 det = this->b; det.at(i) = j, det.at(k) = l; excitations.push_back(Determinant(norb, a, det)); } diff --git a/src/main.cpp b/src/main.cpp index bedf68b..aa8aacc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -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); @@ -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); @@ -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; @@ -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()); - // if nonadiabatic model was specified assign to the solver the nonadiabatic version - if (mdlopt.at("potential").size() > 1) msv = ModelSolver(msnopt.get()); + ModelSolver msv = mdlopt.at("potential").size() == 1 ? ModelSolver(msaopt.get()) : ModelSolver(msnopt.get()); // run the optimization if requested if (mdlopt.at("potential").size() == 1 && input.at("solve").contains("optimize")) { @@ -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()); + + // 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);} } } diff --git a/src/modelsolver.cpp b/src/modelsolver.cpp index e333e9e..10cce8b 100644 --- a/src/modelsolver.cpp +++ b/src/modelsolver.cpp @@ -10,7 +10,7 @@ Result ModelSolver::run(const ModelSystem& system, Result res, bool print) { Result ModelSolver::runad(const ModelSystem& system, Result res, bool print) { // run the classical dynamics if requested - if (!opta.dynamics.gradient.empty()) return runcd(system, opta.dynamics, res, print); + if (!optd.gradient.empty()) return runcd(system, res, print); // define the real and Fourier space along with time and frequency domains res.msv.r = Vector<>(system.ngrid), res.msv.k = Vector<>(system.ngrid), res.msv.t = Vector<>(opta.iters + 1), res.msv.f = Vector<>(opta.iters + opta.spectrum.zeropad + 1); @@ -165,8 +165,8 @@ Result ModelSolver::runad(const ModelSystem& system, Result res, bool print) { } // save the state wavefunction - if (opta.savewfn && system.vars().size() == 2) ModelSystem::SaveWavefunction(ip / ("state" + std::to_string(i) + ".dat"), x.col(0), y.row(0), wfns, energies); - if (opta.savewfn && system.vars().size() == 1) ModelSystem::SaveWavefunction(ip / ("state" + std::to_string(i) + ".dat"), x, wfns, energies); + if (opta.savewfn && system.vars().size() == 2) ModelSystem::SaveWavefunction(ip / ("state" + std::to_string(i + 1) + ".dat"), x.col(0), y.row(0), wfns, energies); + if (opta.savewfn && system.vars().size() == 1) ModelSystem::SaveWavefunction(ip / ("state" + std::to_string(i + 1) + ".dat"), x, wfns, energies); } // print the newline @@ -178,7 +178,7 @@ Result ModelSolver::runad(const ModelSystem& system, Result res, bool print) { Result ModelSolver::runnad(const ModelSystem& system, Result res, bool print) { // run the classical dynamics if requested - if (!optn.dynamics.gradient.empty()) return runcd(system, optn.dynamics, res, print); + if (!optd.gradient.empty()) return runcd(system, res, print); // define the containers for energies, wfns, transforms, density matrices and coordinates and initialize wfn std::vector energies; Matrix> psi(system.ngrid, 2); std::vector>> psis; @@ -196,7 +196,7 @@ Result ModelSolver::runnad(const ModelSystem& system, Result res, bool print) { // fill in the initial wavefunction and normalize it for (size_t i = 0; i < optn.guess.size(); i++) { - psi.block(0, i, system.ngrid, 1) = (Expression(optn.guess.at(i), system.vars()).eval(res.msv.r).array() * (I * optn.momentum * res.msv.r.array()).exp()); + psi.block(0, i, system.ngrid, 1) = (Expression(optn.guess.at(i), system.vars()).eval(res.msv.r).array() * (I * optn.momentum * (res.msv.r.array() + 10)).exp()); } psi.block(0, 0, system.ngrid, 1).array() /= std::sqrt(psi.block(0, 0, system.ngrid, 1).array().abs2().sum() * dr); // create the potential @@ -330,7 +330,7 @@ Result ModelSolver::runnad(const ModelSystem& system, Result res, bool print) { std::vector>> state; for (auto psi : psis) state.push_back(psi.col(i)); // save the wavefunction - ModelSystem::SaveWavefunction(ip / ("state" + std::to_string(i) + ".dat"), res.msv.r, state, energies); + ModelSystem::SaveWavefunction(ip / ("state" + std::to_string(i + 1) + ".dat"), res.msv.r, state, energies); } // print the newline @@ -340,13 +340,14 @@ Result ModelSolver::runnad(const ModelSystem& system, Result res, bool print) { return res; } -template -Result ModelSolver::runcd(const ModelSystem& system, const T& optdyn, Result res, bool print) { +Result ModelSolver::runcd(const ModelSystem& system, Result res, bool print) { // print the header - if (print) std::printf(" ITER TIME [fs] STATE POT [Eh] KIN [Eh] E [Eh] |GRAD| TIME\n"); + if (print) std::printf(" TRAJ ITER TIME [fs] STATE PROB POT [Eh] KIN [Eh] E [Eh] |GRAD| TIME\n"); - // define energy and gradient functions - std::vector energy, grad; + // define energy, gradient and coupling functions and define the array of state populations + std::vector energy, coupling, grad; std::vector> states; + + // fill the real space coordinates res.msv.r = Vector<>(system.ngrid), res.msv.U = Matrix<>(system.ngrid, system.potential.size()); // fill the real space @@ -357,50 +358,88 @@ Result ModelSolver::runcd(const ModelSystem& system, const T& optdyn, Result res // fill the gradient and energy expressions for (size_t i = 0; i < system.potential.size(); i++) { energy.push_back(Expression(system.potential.at(i).at(i), system.vars())); - grad.push_back(Expression(optdyn.gradient.at(i), system.vars())); + grad.push_back(Expression(optd.gradient.at(i), system.vars())); res.msv.U.col(i) = energy.at(i).eval(res.msv.r); } + + // fill the coupling expressions + for (size_t i = 0; i < system.potential.size(); i++) { + for (size_t j = 0; j < system.potential.size(); j++) { + if (i != j) coupling.push_back(Expression(system.potential.at(i).at(j), system.vars())); + } + } + + for (int i = 0; i < optd.trajs; i++) { + // random number generator + std::mt19937 mt(optd.seed + i); std::uniform_real_distribution dist(0, 1); - // define the initial conditions - Matrix<> r(optdyn.iters + 1, system.vars().size()); Vector<> v(system.vars().size()), a(system.vars().size()), m(system.vars().size()); Vector state(optdyn.iters + 1); m.fill(system.mass()); + // define and fill the initial conditions + Matrix<> r(optd.iters + 1, system.vars().size()); Vector<> v(system.vars().size()), a(system.vars().size()), m(system.vars().size()); m.fill(system.mass()); + r.row(0) = Eigen::Map>(optd.position.data(), r.cols()), v = Eigen::Map>(optd.velocity.data(), r.cols()), a.fill(0); - // start the timer, fill the initial conditions - auto start = Timer::Now(); r.row(0) = Eigen::Map>(optdyn.position.data(), optdyn.position.size()), a.fill(0); state(0) = optdyn.state - 1; + r.row(0)(0) = r.row(0)(0) + dist(mt) * 2 - 1; - // fill the initial velocity - v = Eigen::Map>(optdyn.velocity.data(), optdyn.velocity.size()); + // define the initial state, state container and start the timer + Vector state(optd.iters + 1); state(0) = optd.state - 1; auto start = Timer::Now(); - // calculate the force - Vector<> F = -grad.at(state(0)).eval(r.row(0)); double Epot = energy.at(state(0)).get(r.row(0)), Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); + // calculate the force and potential with kinetic energy + Vector<> F = -grad.at(state(0)).eval(r.row(0)); double Epot = energy.at(state(0)).get(r.row(0)), Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); - // print the zeroth iteration - if (print) std::printf("%6d %9.4f %5d %14.8f %14.8f %14.8f %.2e %s\n", 0, 0.0, state(0) + 1, Epot, Ekin, Epot + Ekin, F.norm(), Timer::Format(Timer::Elapsed(start)).c_str()); + // define the energy difference container and fill the zeroth iteration + std::vector Ediff; Ediff.push_back(energy.at(1).get(r.row(0)) - energy.at(0).get(r.row(0))); - for (int i = 0; i < optdyn.iters; i++) { - // start the timer and store the previous v and a - start = Timer::Now(); Vector<> vp = v, ap = a; + // print the zeroth iteration + if (print) std::printf("%6d %6d %9.4f %5d %1.2f %14.8f %14.8f %14.8f %.2e %s\n", i + 1, 0, 0.0, state(0) + 1, 0.0, Epot, Ekin, Epot + Ekin, F.norm(), Timer::Format(Timer::Elapsed(start)).c_str()); - // calculate the velocity and accceleration - a = F.array() / m.array(); v = vp + 0.5 * (ap + a) * optdyn.step; + for (int j = 0; j < optd.iters; j++) { + // start the timer and store the previous v and a + start = Timer::Now(); Vector<> vp = v, ap = a; - // move the system - r.row(i + 1) = r.row(i) + optdyn.step * (v + 0.5 * a * optdyn.step); + // calculate the velocity and accceleration + a = F.array() / m.array(); v = vp + 0.5 * (ap + a) * optd.step; - // fill the state - state(i + 1) = state(i); + // move the system + r.row(j + 1) = r.row(j) + optd.step * (v + 0.5 * a * optd.step); - // calculate the potential and kinetic energy - double Epot = energy.at(state(i + 1)).get(r.row(i + 1)), Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); + // fill the state + state(j + 1) = state(j); - // calculate the force - F = -grad.at(state(i + 1)).eval(r.row(i + 1)); + // calculate the potential and kinetic energy + double Epot = energy.at(state(j + 1)).get(r.row(j + 1)), Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); - // print the iteration - if (print) std::printf("%6d %9.4f %5d %14.8f %14.8f %14.8f %.2e %s\n", i + 1, AU2FS * optdyn.step * (i + 1), state(i + 1) + 1, Epot, Ekin, Epot + Ekin, F.norm(), Timer::Format(Timer::Elapsed(start)).c_str()); + // append the energy and calculate the derivative of the energy difference + Ediff.push_back(energy.at(1).get(r.row(j + 1)) - energy.at(0).get(r.row(j + 1))); double dE = (Ediff.at(j + 1) - Ediff.at(j)) / optd.step; + + // calculate the probability of state change + double gamma = std::pow(coupling.at(0).get(r.row(j + 1)), 2) / std::abs(dE); double P = 1 - std::exp(-2 * M_PI * (std::isnan(gamma) ? 0 : gamma)); + + // change the state if the jump is accepted + if (dist(mt) < P) state(j + 1) = state(j + 1) == 1 ? 0 : 1; + + // calculate the force + F = -grad.at(state(j + 1)).eval(r.row(j + 1)); + + // print the iteration + if (print) std::printf("%6d %6d %9.4f %5d %1.2f %14.8f %14.8f %14.8f %.2e %s\n", i + 1, j + 1, AU2FS * optd.step * (j + 1), state(j + 1) + 1, P, Epot, Ekin, Epot + Ekin, F.norm(), Timer::Format(Timer::Elapsed(start)).c_str()); + } + + // save the states and define the density diagonal + states.push_back(state); + + // save the matrix of states and coordinates + Matrix<> SR(r.rows(), r.cols() + 1); SR << (state.array() + 1).cast(), r; EigenWrite(ip / ("trajectory" + std::to_string(i + 1) + ".mat"), SR); } - // save the matrix of states and coordinates - Matrix<> SR(r.rows(), r.cols() + 1); SR << (state.array() + 1).cast(), r; EigenWrite(ip / "trajectory.mat", SR); + // define the density matrix diagonal array + Matrix<> rho = Matrix<>(optd.iters + 1, system.potential.size() + 1); + + // fill the density matrix diagonal + for (int j = 0; j < optd.iters + 1; j++) { + for (const auto& state : states) { + rho(j, 1) += state(j) == 0 ? 1.0 / states.size() : 0; + } rho(j, 0) = j * optd.step, rho(j, 2) = 1 - rho(j, 1); + } EigenWrite(ip / ("rho.mat"), rho); + // print the newline if (print) std::printf("\n");