diff --git a/include/classicaldynamics.h b/include/classicaldynamics.h index 050762d..3b72d1a 100644 --- a/include/classicaldynamics.h +++ b/include/classicaldynamics.h @@ -9,7 +9,7 @@ class ClassicalDynamics { public: struct TrajectoryData { - Eigen::VectorXi state; Eigen::MatrixXd position, velocity; std::vector diabatic_potential, adiabatic_potential; std::vector hopping_geometry; std::vector hopping_time; + Eigen::VectorXi state; Eigen::MatrixXd position, velocity; std::vector diabatic_potential, adiabatic_potential, tdc; std::vector hopping_geometry; std::vector hopping_time; }; ClassicalDynamics(const Input::ClassicalDynamics& input) : input(input) {} diff --git a/include/fewestswitches.h b/include/fewestswitches.h index ced0c62..8b97c2d 100644 --- a/include/fewestswitches.h +++ b/include/fewestswitches.h @@ -4,6 +4,8 @@ #include #include +// #include "export.h" + class FewestSwitches { public: FewestSwitches(const Input::ClassicalDynamics::SurfaceHopping& input, bool adiabatic, int seed); @@ -14,7 +16,14 @@ class FewestSwitches { static Eigen::VectorXcd calculate_population_product(const Eigen::VectorXcd& population, const Eigen::VectorXd& potential, const Eigen::MatrixXd& derivative_coupling); static Eigen::VectorXcd propagate_population(const Eigen::VectorXcd& population, const Eigen::VectorXd& potential, const Eigen::MatrixXd& derivative_coupling, double time_step); - std::tuple jump(Eigen::VectorXcd population, const std::vector& phi_vector, const std::vector& potential_vector, int iteration, int state, double time_step); + std::tuple jump(Eigen::VectorXcd population, const std::vector& phi_vector, const std::vector& potential_vector, int iteration, int state, double time_step); + + // ~FewestSwitches() { + // Eigen::VectorXd data = Eigen::Map(offtdc.data(), offtdc.size()); + // + // Export::EigenMatrixDouble(std::string("OFFTDC-") + (input.kappa ? "K" : "") + "FSSH.mat", data, Eigen::VectorXd::LinSpaced(offtdc.size(), 1, offtdc.size())); + // } + // std::vector offtdc; private: bool adiabatic; std::uniform_real_distribution dist; std::mt19937 mt; Input::ClassicalDynamics::SurfaceHopping input; diff --git a/include/input.h b/include/input.h index 711e1f3..b73dab7 100644 --- a/include/input.h +++ b/include/input.h @@ -45,7 +45,7 @@ struct Input { } quantum_dynamics; struct ClassicalDynamics { struct DataExport { - bool diabatic_population, adiabatic_population, total_energy, total_energy_mean, position, position_mean, momentum, momentum_mean, potential_energy, potential_energy_mean, kinetic_energy, kinetic_energy_mean, hopping_geometry, hopping_time; + bool diabatic_population, adiabatic_population, total_energy, total_energy_mean, position, position_mean, momentum, momentum_mean, potential_energy, potential_energy_mean, kinetic_energy, kinetic_energy_mean, hopping_geometry, hopping_time, tdc, tdc_mean; } data_export; struct InitialConditions { std::vector> positions, momenta, position_distribution, momentum_distribution; int diabatic_state; int adiabatic_state; int mass; @@ -75,7 +75,7 @@ NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::HartreeFock::MollerPlesset, order); NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::HartreeFock, diis_size, max_iter, threshold, configuration_interaction, coupled_cluster, moller_plesset, generalized, gradient); NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::QuantumDynamics::DataExport, diabatic_wavefunction, adiabatic_wavefunction, diabatic_density, adiabatic_density, diabatic_population, adiabatic_population, diabatic_potential, adiabatic_potential, total_energy, position, momentum, acf, potential_energy, kinetic_energy); NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::QuantumDynamics, potential, constants, imaginary, real, iterations, time_step, data_export); -NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::ClassicalDynamics::DataExport, total_energy, position, momentum, diabatic_population, adiabatic_population, total_energy_mean, position_mean, momentum_mean, potential_energy, potential_energy_mean, kinetic_energy, kinetic_energy_mean, hopping_geometry, hopping_time); +NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::ClassicalDynamics::DataExport, total_energy, position, momentum, diabatic_population, adiabatic_population, total_energy_mean, position_mean, momentum_mean, potential_energy, potential_energy_mean, kinetic_energy, kinetic_energy_mean, hopping_geometry, hopping_time, tdc, tdc_mean); NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::ClassicalDynamics::InitialConditions, positions, momenta, diabatic_state, adiabatic_state, mass, position_distribution, momentum_distribution); NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::ClassicalDynamics::SurfaceHopping, kappa, type, quantum_step_factor); NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Input::ClassicalDynamics, potential, constants, iterations, trajectories, time_step, data_export, seed, adiabatic, surface_hopping, log_interval_step, log_interval_traj, initial_conditions); diff --git a/research/fssh_vs_lzsh/fssh_vs_lzsh.sh b/research/fssh_vs_lzsh/fssh_vs_lzsh.sh index e97803e..72f7b99 100755 --- a/research/fssh_vs_lzsh/fssh_vs_lzsh.sh +++ b/research/fssh_vs_lzsh/fssh_vs_lzsh.sh @@ -256,13 +256,13 @@ for MODEL in ${MODELS[@]}; do $PLOT_1D POTENTIAL-ENERGY_EXACT-REAL_1.mat POTENTIAL-ENERGY-MEAN_FS-ADIABATIC.mat POTENTIAL-ENERGY-MEAN_KFS-ADIABATIC.mat POTENTIAL-ENERGY-MEAN_LZ-ADIABATIC.mat --legend "EXACT" "FSSH" "KFSSH" "LZSH" --title "POTENTIAL ENERGY: ${MODEL}" --xlabel "Time (a.u.)" --ylabel "ENERGY (a.u.)" --output "POTENTIAL-ENERGY_${MODEL}_P=${MOMENTUM}" --png # plot the hopping geometries - $PLOT_1D HOPPING-GEOMETRIES_FS-ADIABATIC.mat HOPPING-GEOMETRIES_KFS-ADIABATIC.mat --bins 100 --legend "FSSH" "KFSSH" --title "HOPPING GEOMETRIES: ${MODEL}" --xlabel "Coordinate (a.u.)" --ylabel "Relative Count" --output "HOPPING-GEOMETRIES_${MODEL}" --histogram --png - $PLOT_1D HOPPING-TIMES_FS-ADIABATIC.mat HOPPING-TIMES_KFS-ADIABATIC.mat --bins 100 --legend "FSSH" "KFSSH" --title "HOPPING TIMES: ${MODEL}" --xlabel "Time (a.u.)" --ylabel "Relative Count" --output "HOPPING-TIMES_${MODEL}" --histogram --png + $PLOT_1D HOPPING-GEOMETRIES_FS-ADIABATIC.mat HOPPING-GEOMETRIES_KFS-ADIABATIC.mat --bins 100 --legend "FSSH" "KFSSH" --title "HOPPING GEOMETRIES: ${MODEL}" --xlabel "Coordinate (a.u.)" --ylabel "Relative Count" --output "HOPPING-GEOMETRIES_${MODEL}_P=${MOMENTUM}" --histogram --png + $PLOT_1D HOPPING-TIMES_FS-ADIABATIC.mat HOPPING-TIMES_KFS-ADIABATIC.mat --bins 100 --legend "FSSH" "KFSSH" --title "HOPPING TIMES: ${MODEL}" --xlabel "Time (a.u.)" --ylabel "Relative Count" --output "HOPPING-TIMES_${MODEL}_P=${MOMENTUM}" --histogram --png # make the trajectory analysis image montage "POTENTIAL-ADIABATIC_${MODEL}.png" "POPULATION-ADIABATIC_${MODEL}_P=${MOMENTUM}.png" "POSITION_${MODEL}_P=${MOMENTUM}.png" "MOMENTUM_${MODEL}_P=${MOMENTUM}.png" -mode concatenate -tile x1 "TRAJECTORIES-GENERAL_${MODEL}_P=${MOMENTUM}.png" montage "POTENTIAL-ADIABATIC_${MODEL}.png" "KINETIC-ENERGY_${MODEL}_P=${MOMENTUM}.png" "POTENTIAL-ENERGY_${MODEL}_P=${MOMENTUM}.png" "TOTAL-ENERGY_${MODEL}_P=${MOMENTUM}.png" -mode concatenate -tile x1 "TRAJECTORIES-ENERGETICS_${MODEL}_P=${MOMENTUM}.png" - montage "POTENTIAL-ADIABATIC_${MODEL}.png" "HOPPING-GEOMETRIES_${MODEL}.png" "HOPPING-TIMES_${MODEL}.png" -mode concatenate -tile x1 "TRAJECTORIES-TRANSITIONS_${MODEL}_P=${MOMENTUM}.png" + montage "POTENTIAL-ADIABATIC_${MODEL}.png" "HOPPING-GEOMETRIES_${MODEL}_P=${MOMENTUM}.png" "HOPPING-TIMES_${MODEL}_P=${MOMENTUM}.png" -mode concatenate -tile x1 "TRAJECTORIES-TRANSITIONS_${MODEL}_P=${MOMENTUM}.png" done # plot the population dependence on momentum diff --git a/research/sodium_iodide/check/fssh-check.json b/research/sodium_iodide/check/fssh-check.json new file mode 100644 index 0000000..e7d7ded --- /dev/null +++ b/research/sodium_iodide/check/fssh-check.json @@ -0,0 +1,36 @@ +{ + "classical_dynamics" : { + "initial_conditions" : { + "positions" : [ + [4.966781744519960462] + ], + "momenta" : [ + [-7.511477096479932669] + ], + "adiabatic_state" : 2, + "mass" : 35480.251398 + }, + "potential" : [ + ["((A2+(B2/(x*au2a))^8)*exp(-(x*au2a)/rho)-e2/(x*au2a)-e2*(lp+lm)/2/(x*au2a)^4-C2/(x*au2a)^6-2*e2*lm*lp/(x*au2a)^7+de)/au2ev", "A12*exp(-beta12*((x*au2a)-Rx)^2)/au2ev"], + ["A12*exp(-beta12*((x*au2a)-Rx)^2)/au2ev", "A1*exp(-beta1*((x*au2a)-R0))/au2ev" ] + ], + "constants" : { + "A2" : 2760, "B2" : 2.398, "C2" : 11.3, "rho" : 0.3489, "e2" : 14.3996, "lp" : 0.408, "lm" : 6.431, + "A1" : 0.813, "beta1" : 4.08, "A12" : 0.055, "beta12" : 0.6931, "R0" : 2.67, "Rx" : 6.93, "de" : 2.075, + + "au2ev" : 27.21138602, "au2a" : 0.529177 + }, + "surface_hopping" : { + "type" : "fewest-switches", + "kappa" : false + }, + "adiabatic" : true, + "iterations" : 12000, + "time_step" : 1, + "trajectories" : 1, + "log_interval_step" : 1, + "data_export" : { + "diabatic_population" : true, "adiabatic_population" : true, "tdc" : true + } + } +} diff --git a/research/sodium_iodide/kfssh-check.json b/research/sodium_iodide/check/kfssh-check.json similarity index 96% rename from research/sodium_iodide/kfssh-check.json rename to research/sodium_iodide/check/kfssh-check.json index 0c22723..631cda0 100644 --- a/research/sodium_iodide/kfssh-check.json +++ b/research/sodium_iodide/check/kfssh-check.json @@ -25,12 +25,12 @@ "kappa" : true }, "adiabatic" : true, - "iterations" : 9000, + "iterations" : 12000, "time_step" : 1, "trajectories" : 1, "log_interval_step" : 1, "data_export" : { - "diabatic_population" : true, "adiabatic_population" : true + "diabatic_population" : true, "adiabatic_population" : true, "tdc" : true } } } diff --git a/research/sodium_iodide/fssh.json b/research/sodium_iodide/fssh.json index a6a90ad..cf33e13 100644 --- a/research/sodium_iodide/fssh.json +++ b/research/sodium_iodide/fssh.json @@ -22,9 +22,9 @@ "type" : "fewest-switches" }, "adiabatic" : true, - "iterations" : 9000, + "iterations" : 100000, "time_step" : 1, - "trajectories" : 1000, + "trajectories" : 100, "log_interval_step" : 500, "data_export" : { "diabatic_population" : true, "adiabatic_population" : true diff --git a/research/sodium_iodide/kfssh.json b/research/sodium_iodide/kfssh.json index 4bd8376..16d7742 100644 --- a/research/sodium_iodide/kfssh.json +++ b/research/sodium_iodide/kfssh.json @@ -23,9 +23,9 @@ "kappa" : true }, "adiabatic" : true, - "iterations" : 9000, + "iterations" : 100000, "time_step" : 1, - "trajectories" : 1000, + "trajectories" : 100, "log_interval_step" : 500, "data_export" : { "diabatic_population" : true, "adiabatic_population" : true diff --git a/research/sodium_iodide/lzsh.json b/research/sodium_iodide/lzsh.json index 7c22ba4..ddd3e0f 100644 --- a/research/sodium_iodide/lzsh.json +++ b/research/sodium_iodide/lzsh.json @@ -22,9 +22,9 @@ "type" : "landau-zener" }, "adiabatic" : true, - "iterations" : 9000, + "iterations" : 100000, "time_step" : 1, - "trajectories" : 1000, + "trajectories" : 100, "log_interval_step" : 500, "data_export" : { "diabatic_population" : true, "adiabatic_population" : true diff --git a/research/sodium_iodide/qd.json b/research/sodium_iodide/qd.json index dc0e010..e1c65f0 100644 --- a/research/sodium_iodide/qd.json +++ b/research/sodium_iodide/qd.json @@ -19,7 +19,7 @@ "au2ev" : 27.21138602, "au2a" : 0.529177 }, "real" : 1, - "iterations" : 9000, + "iterations" : 100000, "time_step" : 1, "data_export" : { "diabatic_population" : true, "adiabatic_population" : true, diff --git a/src/classicaldynamics.cpp b/src/classicaldynamics.cpp index 3b08351..424e18c 100644 --- a/src/classicaldynamics.cpp +++ b/src/classicaldynamics.cpp @@ -189,14 +189,14 @@ void ClassicalDynamics::run(const Input::Wavefunction& initial_diabatic_wavefunc // define and initialize the surface hopping algorithms FewestSwitches fewestswitches(input.surface_hopping, input.adiabatic, seed); LandauZener landauzener(input.surface_hopping, input.adiabatic, seed); - // define the current fssh population vector + // define the current fssh population vector and TDC Eigen::VectorXcd fssh_population = Eigen::VectorXcd::Zero(input.potential.size()); fssh_population(state(0)) = 1; // define the vector of hopping geometries and times std::vector hopping_geometry_vector; std::vector hopping_time_vector; // define a vector to contain diabatic potentials and eigenvectors phi - std::vector diabatic_potential_vector(input.iterations + 1), adiabatic_potential_vector(input.iterations + 1), phi_vector(input.iterations + 1); + std::vector diabatic_potential_vector(input.iterations + 1), adiabatic_potential_vector(input.iterations + 1), phi_vector(input.iterations + 1), tdc_vector(input.iterations); // propagate the current trajectory for (int j = 0; j < input.iterations + 1; j++) { @@ -211,7 +211,7 @@ void ClassicalDynamics::run(const Input::Wavefunction& initial_diabatic_wavefunc if (j) {position.row(j) = position.row(j - 1) + input.time_step * (velocity.row(j) + 0.5 * acceleration.row(j) * input.time_step), state(j) = state(j - 1);} // calculate the potential at the current point and assign it to the container and create the new state variable - Eigen::MatrixXd potential = evaluate_potential(potential_expressions, position.row(j)); diabatic_potential_vector.at(j) = potential; int new_state = state(j); + Eigen::MatrixXd potential = evaluate_potential(potential_expressions, position.row(j)), tdc; diabatic_potential_vector.at(j) = potential; int new_state = state(j); // adiabatization block if (input.adiabatic) { @@ -227,7 +227,7 @@ void ClassicalDynamics::run(const Input::Wavefunction& initial_diabatic_wavefunc if (input.surface_hopping.type == "landau-zener" && j > 1) { new_state = landauzener.jump(input.adiabatic ? adiabatic_potential_vector : diabatic_potential_vector, j, state(j), input.time_step); } else if (input.surface_hopping.type == "fewest-switches" && j) { - std::tie(fssh_population, new_state) = fewestswitches.jump(fssh_population, phi_vector, adiabatic_potential_vector, j, state(j), input.time_step); + std::tie(tdc, fssh_population, new_state) = fewestswitches.jump(fssh_population, phi_vector, adiabatic_potential_vector, j, state(j), input.time_step); } // update the velocity and the state @@ -238,6 +238,9 @@ void ClassicalDynamics::run(const Input::Wavefunction& initial_diabatic_wavefunc // append the hopping geometry and time if (j && state(j) != state(j - 1) && input.data_export.hopping_geometry) hopping_geometry_vector.push_back(position.row(j)); if (j && state(j) != state(j - 1) && input.data_export.hopping_time ) hopping_time_vector.push_back(j * input.time_step); + + // append the TDC + if (j && (input.data_export.tdc || input.data_export.tdc_mean)) tdc_vector.at(j - 1) = tdc; // print the iteration info if ((j % input.log_interval_step == 0 || (j && state(j) != state(j - 1))) && (i ? i + 1 : i) % input.log_interval_traj == 0) { @@ -246,7 +249,7 @@ void ClassicalDynamics::run(const Input::Wavefunction& initial_diabatic_wavefunc } // set the trajectory data - trajectory_data_vector.at(i) = {state, position, velocity, diabatic_potential_vector, adiabatic_potential_vector, hopping_geometry_vector, hopping_time_vector}; + trajectory_data_vector.at(i) = {state, position, velocity, diabatic_potential_vector, adiabatic_potential_vector, tdc_vector, hopping_geometry_vector, hopping_time_vector}; } // create the vector of populations diff --git a/src/export.cpp b/src/export.cpp index a912cec..7290df8 100644 --- a/src/export.cpp +++ b/src/export.cpp @@ -241,10 +241,10 @@ void Export::ClassicalTrajectories(const Input::ClassicalDynamics& input, const } // transform the hopping geometries to a matrix - data_matrix = Eigen::MatrixXd::Zero(hopping_geometries.size(), hopping_geometries.front().rows()); for (size_t i = 0; i < hopping_geometries.size(); i++) data_matrix.row(i) = hopping_geometries.at(i).transpose(); + if (!hopping_geometries.empty()) data_matrix = Eigen::MatrixXd::Zero(hopping_geometries.size(), hopping_geometries.front().rows()); for (size_t i = 0; i < hopping_geometries.size(); i++) data_matrix.row(i) = hopping_geometries.at(i).transpose(); // export the energy mean - Export::EigenMatrixDouble(std::string("HOPPING-GEOMETRIES_") + algorithm + (input.adiabatic ? "-ADIABATIC" : "-DIABATIC") + ".mat", data_matrix); + if (!hopping_geometries.empty()) Export::EigenMatrixDouble(std::string("HOPPING-GEOMETRIES_") + algorithm + (input.adiabatic ? "-ADIABATIC" : "-DIABATIC") + ".mat", data_matrix); } // export the hopping times @@ -259,10 +259,46 @@ void Export::ClassicalTrajectories(const Input::ClassicalDynamics& input, const } // transform the hopping geometries to a matrix - data_matrix = Eigen::MatrixXd::Zero(hopping_times.size(), 1); for (size_t i = 0; i < hopping_times.size(); i++) data_matrix(i, 0) = hopping_times.at(i); + if (!hopping_times.empty()) data_matrix = Eigen::MatrixXd::Zero(hopping_times.size(), 1); for (size_t i = 0; i < hopping_times.size(); i++) data_matrix(i, 0) = hopping_times.at(i); // export the energy mean - Export::EigenMatrixDouble(std::string("HOPPING-TIMES_") + algorithm + (input.adiabatic ? "-ADIABATIC" : "-DIABATIC") + ".mat", data_matrix); + if (!hopping_times.empty()) Export::EigenMatrixDouble(std::string("HOPPING-TIMES_") + algorithm + (input.adiabatic ? "-ADIABATIC" : "-DIABATIC") + ".mat", data_matrix); + } + + // export the TDC + if (Eigen::VectorXd time_variable = Eigen::VectorXd::LinSpaced(input.iterations, 1, input.time_step * input.iterations); input.data_export.tdc) { + + // create the matrix containing the TDC + data_matrix = Eigen::MatrixXd::Zero(input.iterations, trajectory_data_vector.size() * trajectory_data_vector.at(0).tdc.at(0).size()); + + // extract the tdc dimension + int dim = trajectory_data_vector.at(0).tdc.at(0).rows(); + + // fill the matrix with the TDC + for (int i = 0; i < input.trajectories; i++) for (int j = 0; j < input.iterations; j++) for (int k = 0; k < dim * dim; k++) { + data_matrix.row(j)(i * dim * dim + k) = trajectory_data_vector.at(i).tdc.at(j)(k / dim, k % dim); + } + + // export the energy + Export::EigenMatrixDouble(std::string("TDC_") + algorithm + (input.adiabatic ? "-ADIABATIC" : "-DIABATIC") + ".mat", data_matrix, time_variable); + } + + // export the mean TDC + if (Eigen::VectorXd time_variable = Eigen::VectorXd::LinSpaced(input.iterations, 1, input.time_step * input.iterations); input.data_export.tdc_mean) { + + // create the matrix containing the TDC mean + data_matrix = Eigen::MatrixXd::Zero(input.iterations, trajectory_data_vector.at(0).tdc.at(0).size()); + + // extract the tdc dimension + int dim = trajectory_data_vector.at(0).tdc.at(0).rows(); + + // fill the matrix with the the TDC + for (int i = 0; i < input.trajectories; i++) for (int j = 0; j < input.iterations; j++) for (int k = 0; k < dim * dim; k++) { + data_matrix.row(j)(k) += trajectory_data_vector.at(i).tdc.at(j)(k / dim, k % dim); + } + + // export the energy + Export::EigenMatrixDouble(std::string("TDC_") + algorithm + (input.adiabatic ? "-ADIABATIC" : "-DIABATIC") + ".mat", data_matrix / input.trajectories, time_variable); } } diff --git a/src/fewestswitches.cpp b/src/fewestswitches.cpp index e75fa61..8f2a8b1 100644 --- a/src/fewestswitches.cpp +++ b/src/fewestswitches.cpp @@ -48,7 +48,7 @@ Eigen::MatrixXd FewestSwitches::calculate_derivative_coupling_kappa(const std::v } // return the derivative coupling matrix - return derivative_coupling - derivative_coupling.transpose(); + return derivative_coupling.transpose() - derivative_coupling; } std::vector> FewestSwitches::calculate_hopping_probabilities(const Eigen::VectorXcd& ci, const Eigen::MatrixXd& derivative_coupling, int state, double time_step) { @@ -79,7 +79,7 @@ Eigen::VectorXcd FewestSwitches::propagate_population(const Eigen::VectorXcd& po return population + time_step / 6 * (k1 + 2 * k2 + 2 * k3 + k4); } -std::tuple FewestSwitches::jump(Eigen::VectorXcd population, const std::vector& phi_vector, const std::vector& potential_vector, int iteration, int state, double time_step) { +std::tuple FewestSwitches::jump(Eigen::VectorXcd population, const std::vector& phi_vector, const std::vector& potential_vector, int iteration, int state, double time_step) { // define the derivative coupling and the new state Eigen::MatrixXd derivative_coupling; int new_state = state; @@ -88,6 +88,8 @@ std::tuple FewestSwitches::jump(Eigen::VectorXcd populati derivative_coupling = calculate_derivative_coupling_kappa(potential_vector, iteration, time_step); } else derivative_coupling = calculate_derivative_coupling(phi_vector.at(iteration), phi_vector.at(iteration - 1), time_step); + // offtdc.push_back(derivative_coupling(0, 1)); + // pripagate the populations for (int k = 0; k < input.quantum_step_factor; k++) { @@ -113,5 +115,5 @@ std::tuple FewestSwitches::jump(Eigen::VectorXcd populati } // return the propagated population and the new state - return {population, new_state}; + return {derivative_coupling, population, new_state}; } diff --git a/src/input.cpp b/src/input.cpp index fce52f7..15e1689 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -111,7 +111,9 @@ nlohmann::json default_input = R"( "kinetic_energy" : false, "kinetic_energy_mean" : false, "hopping_geometry" : false, - "hopping_time" : false + "hopping_time" : false, + "tdc" : false, + "tdc_mean" : false } } }