Skip to content

Commit

Permalink
remove classical adiabatic dynamics
Browse files Browse the repository at this point in the history 
The adiabatic populations are now transformed from the diabatic ones.
Also add an example of three state crossing potential and a simple LZ
fix for that (LZ probabilities are sorted).
  • Loading branch information
@tjira
tjira committed Jun 19, 2024
1 parent d38833a commit 762626d
Show file tree
Hide file tree
Showing 9 changed files with 215 additions and 125 deletions.
30 changes: 12 additions & 18 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,23 +118,17 @@ set_property(TEST qdyn_3d_HO_imaginary PROPERTY PASS_REGULAR_EXPRESSION "FINAL W
set_property(TEST qdyn_3d_HO_real PROPERTY PASS_REGULAR_EXPRESSION "FINAL WFN 0 ENERGY: 3.374795" )

# add random tests
add_test(NAME random_qdyn_1d_ds_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_qdyn_1d_ds_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_qdyn_1d_ts_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_qdyn_1d_ts_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_qdyn_1d_tully_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_qdyn_1d_tully_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_diabatic_ds_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_diabatic_ds_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_diabatic_ts_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_diabatic_ts_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_diabatic_tully_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_diabatic_tully_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_adiabatic_ds_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_adiabatic_ds_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_adiabatic_ts_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_adiabatic_ts_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_adiabatic_tully_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_adiabatic_tully_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_qdyn_1d_tully_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_qdyn_1d_tully_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_tully_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_tully_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_qdyn_1d_ds_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_qdyn_1d_ds_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_ds_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_ds_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_cdyn_1d_ts_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_cdyn_1d_ts_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)
add_test(NAME random_qdyn_1d_ts_1 COMMAND ${CMAKE_MAKE_PROGRAM} random_qdyn_1d_ts_1 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/example)

# add expected results of random tests
set_property(TEST random_qdyn_1d_ds_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.052824")
set_property(TEST random_qdyn_1d_ts_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.720077")
set_property(TEST random_qdyn_1d_tully_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.676839")
set_property(TEST random_cdyn_1d_adiabatic_ds_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.072000")
set_property(TEST random_cdyn_1d_adiabatic_ts_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.713000")
set_property(TEST random_cdyn_1d_adiabatic_tully_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.504000")
set_property(TEST random_cdyn_1d_diabatic_ds_1 PROPERTY PASS_REGULAR_EXPRESSION "DIABATIC STATE 0 POP: 0.917000" )
set_property(TEST random_cdyn_1d_diabatic_ts_1 PROPERTY PASS_REGULAR_EXPRESSION "DIABATIC STATE 0 POP: 0.177000" )
set_property(TEST random_cdyn_1d_diabatic_tully_1 PROPERTY PASS_REGULAR_EXPRESSION "DIABATIC STATE 0 POP: 0.445000" )
set_property(TEST random_qdyn_1d_tully_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.676839")
set_property(TEST random_cdyn_1d_tully_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.496000")
set_property(TEST random_qdyn_1d_ds_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.052824")
set_property(TEST random_cdyn_1d_ds_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.066000")
set_property(TEST random_qdyn_1d_ts_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.720077")
set_property(TEST random_cdyn_1d_ts_1 PROPERTY PASS_REGULAR_EXPRESSION "ADIABATIC STATE 0 POP: 0.708000")
14 changes: 14 additions & 0 deletions education/mathematica/adiabatize.wls
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
#!/usr/bin/env wolframscript
(* ::Package:: *)

ClearAll["Global`*"];


Ud[x_]:={
{0.01*Tanh[0.6*x],0.001*Exp[-x^2],0.001*Exp[-x^2]},
{0.001*Exp[-x^2],0,0.001*Exp[-x^2]},
{0.001*Exp[-x^2],0.001*Exp[-x^2],-0.01*Tanh[0.6*x]}
};Ud[x]//MatrixForm


GraphicsRow[{Plot[Evaluate@Diagonal@Ud[x],{x,-8,8}],Plot[Evaluate@Sort@Eigenvalues[Ud[x]],{x,-8,8}]}]
14 changes: 14 additions & 0 deletions education/mathematica/findiff.wls
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
#!/usr/bin/env wolframscript
(* ::Package:: *)

ClearAll["Global`*"]


points={0,-1,-2,-3,-4};order=2;


A=Join[{ConstantArray[1,Length@points]},Table[points^k,{k,1,Length@points-1}]];A//MatrixForm
b=D[Table[x^k,{k,0,Length@points-1}],{x,order}]/.x->points[[1]];b//MatrixForm


LinearSolve[A,b]
10 changes: 10 additions & 0 deletions education/mathematica/wigner.wls
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
#!/usr/bin/env wolframscript
(* ::Package:: *)

ClearAll["Global`*"]


\[Psi][x_]=Exp[-(x-x0)^2]*Exp[+I*p0*x]


W[x_,p_]=Factor@Simplify@Integrate[Simplify[Refine[\[Psi][x+y]\[Conjugate],{p0\[Element]Reals,x\[Element]Reals,x0\[Element]Reals,y\[Element]Reals}]*\[Psi][x-y]]*Exp[2*I*p*y],{y,-Infinity,Infinity}]
56 changes: 35 additions & 21 deletions example/makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,50 +74,64 @@ mp2: molecule
mp3: molecule
"../bin/acorn_integral" -b sto-3g && ../bin/acorn_hf -d 5 -i 100 -t 1e-8 && ../bin/acorn_transform -s && ../bin/acorn_mp -o 3

# target to perform a MP3 calculation
# target to perform a FCI calculation
fci: molecule
"../bin/acorn_integral" -b sto-3g && ../bin/acorn_hf -d 5 -i 100 -t 1e-8 && ../bin/acorn_transform -s && ../bin/acorn_ci

# RANDOM TARGETS ===================================================================================

random_cdyn_1d_adiabatic_ds_1:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u "0.01*tanh(0.6*x)" "0.001*exp(-x^2)" "0.001*exp(-x^2)" "(-0.01)*tanh(0.6*x)" --adiabatic
TS_1 = "0.03*(tanh(1.6*x)+tanh(1.6*(x+7)))" "0.005*exp(-x^2)" "0.005*exp(-(x+7)^2)" \
"0.005*exp(-x^2)" "(-0.03)*(tanh(1.6*x)+tanh(1.6*(x-7)))" "0.005*exp(-(x-7)^2)" \
"0.005*exp(-(x+7)^2)" "0.005*exp(-(x-7)^2)" "(-0.03)*(tanh(1.6*(x+7))-tanh(1.6*(x-7)))"

TS_2 = "0.01*tanh(0.6*x)" "0.001*exp(-x^2)" "0.001*exp(-x^2)" \
"0.001*exp(-x^2)" "0" "0.001*exp(-x^2)" \
"0.001*exp(-x^2)" "0.001*exp(-x^2)" "(-0.01)*tanh(0.6*x)"

TULLY_1 = "if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))" "0.005*exp(-x^2)" \
"0.005*exp(-x^2)" "(-1)*if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))"

random_cdyn_1d_diabatic_ds_1:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u "0.01*tanh(0.6*x)" "0.001*exp(-x^2)" "0.001*exp(-x^2)" "(-0.01)*tanh(0.6*x)"
DS_1 = "0.01*tanh(0.6*x)" "0.001*exp(-x^2)" \
"0.001*exp(-x^2)" "(-0.01)*tanh(0.6*x)"

random_cdyn_1d_adiabatic_ts_1:
"../bin/acorn_cdyn" -c -10 -e 2 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u "0.03*(tanh(1.6*x)+tanh(1.6*(x+7)))" "0.005*exp(-x^2)" "0.005*exp(-(x+7)^2)" "0.005*exp(-x^2)" "(-0.03)*(tanh(1.6*x)+tanh(1.6*(x-7)))" "0.005*exp(-(x-7)^2)" "0.005*exp(-(x+7)^2)" "0.005*exp(-(x-7)^2)" "(-0.03)*(tanh(1.6*(x+7))-tanh(1.6*(x-7)))" --adiabatic
random_cdyn_1d_tully_1:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u $(TULLY_1) --adiabatic

random_cdyn_1d_diabatic_ts_1:
"../bin/acorn_cdyn" -c -10 -e 1 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u "0.03*(tanh(1.6*x)+tanh(1.6*(x+7)))" "0.005*exp(-x^2)" "0.005*exp(-(x+7)^2)" "0.005*exp(-x^2)" "(-0.03)*(tanh(1.6*x)+tanh(1.6*(x-7)))" "0.005*exp(-(x-7)^2)" "0.005*exp(-(x+7)^2)" "0.005*exp(-(x-7)^2)" "(-0.03)*(tanh(1.6*(x+7))-tanh(1.6*(x-7)))"
random_cdyn_1d_ds_1:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u $(DS_1) --adiabatic

random_cdyn_1d_adiabatic_tully_1:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u "if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))" "0.005*exp(-x^2)" "0.005*exp(-x^2)" "(-1)*if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))" --adiabatic
random_cdyn_1d_ts_1:
"../bin/acorn_cdyn" -c -10 -e 1 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u $(TS_1) --adiabatic

random_cdyn_1d_diabatic_tully_1:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u "if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))" "0.005*exp(-x^2)" "0.005*exp(-x^2)" "(-1)*if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))"
random_cdyn_1d_ts_2:
"../bin/acorn_cdyn" -c -10 -e 0 -i 3500 -l 500 -m 2000 -p 15 -r 1 -s 1 -t 1000 -u $(TS_2) --adiabatic

random_qdyn_1d_tully_1:
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e $(TULLY_1)
"../bin/acorn_expression" -g -24 40 -o PSI_DIA_GUESS.mat -p 1024 -e "exp(-(x+10)^2)" "0" "0" "0"
"../bin/acorn_qdyn" -i 350 -m 2000 -p 15 -s 10 --adiabatic --savewfn

random_qdyn_1d_ds_1:
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e "0.01*tanh(0.6*x)" "0.001*exp(-x^2)" "0.001*exp(-x^2)" "(-0.01)*tanh(0.6*x)"
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e $(DS_1)
"../bin/acorn_expression" -g -24 40 -o PSI_DIA_GUESS.mat -p 1024 -e "exp(-(x+10)^2)" "0" "0" "0"
"../bin/acorn_qdyn" -i 350 -m 2000 -p 15 -s 10 --adiabatic --savewfn

random_qdyn_1d_ts_1:
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e "0.03*(tanh(1.6*x)+tanh(1.6*(x+7)))" "0.005*exp(-x^2)" "0.005*exp(-(x+7)^2)" "0.005*exp(-x^2)" "(-0.03)*(tanh(1.6*x)+tanh(1.6*(x-7)))" "0.005*exp(-(x-7)^2)" "0.005*exp(-(x+7)^2)" "0.005*exp(-(x-7)^2)" "(-0.03)*(tanh(1.6*(x+7))-tanh(1.6*(x-7)))"
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e $(TS_1)
"../bin/acorn_expression" -g -24 40 -o PSI_DIA_GUESS.mat -p 1024 -e "0" "0" "exp(-(x+10)^2)" "0" "0" "0"
"../bin/acorn_qdyn" -i 350 -m 2000 -p 15 -s 10 --adiabatic --savewfn

random_qdyn_1d_tully_1:
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e "if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))" "0.005*exp(-x^2)" "0.005*exp(-x^2)" "(-1)*if(x>0,0.01*(1-exp(-1.6*x)),-0.01*(1-exp(1.6*x)))"
"../bin/acorn_expression" -g -24 40 -o PSI_DIA_GUESS.mat -p 1024 -e "exp(-(x+10)^2)" "0" "0" "0"
random_qdyn_1d_ts_2:
"../bin/acorn_expression" -g -24 40 -o U_DIA.mat -p 1024 -e $(TS_2)
"../bin/acorn_expression" -g -24 40 -o PSI_DIA_GUESS.mat -p 1024 -e "exp(-(x+10)^2)" "0" "0" "0" "0" "0"
"../bin/acorn_qdyn" -i 350 -m 2000 -p 15 -s 10 --adiabatic --savewfn

random_qdyn_hydrogen_imaginary:
"../bin/acorn_expression" -d 3 -g -16 16 -o U_DIA.mat -p 128 -e "(-1)/sqrt(x^2+y^2+z^2)"
"../bin/acorn_expression" -d 3 -g -16 16 -o PSI_DIA_GUESS.mat -p 128 -e "exp(-(x^2+y^2+z^2))" "0"
"../bin/acorn_qdyn" -d 3 -i 1000 -m 0.999455 -o 1 -p 0 -s 0.1

random_dyn_1d_ds_1: random_qdyn_1d_ds_1 random_cdyn_1d_diabatic_ds_1 random_cdyn_1d_adiabatic_ds_1
random_dyn_1d_ts_1: random_qdyn_1d_ts_1 random_cdyn_1d_diabatic_ts_1 random_cdyn_1d_adiabatic_ts_1
random_dyn_1d_tully_1: random_qdyn_1d_tully_1 random_cdyn_1d_diabatic_tully_1 random_cdyn_1d_adiabatic_tully_1
random_dyn_1d_tully_1: random_qdyn_1d_tully_1 random_cdyn_1d_tully_1
random_dyn_1d_ds_1: random_qdyn_1d_ds_1 random_cdyn_1d_ds_1
random_dyn_1d_ts_1: random_qdyn_1d_ts_1 random_cdyn_1d_ts_1
random_dyn_1d_ts_2: random_qdyn_1d_ts_2 random_cdyn_1d_ts_2
13 changes: 6 additions & 7 deletions program/cdyn/include/landauzener.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,14 @@

class LandauZener {
public:
LandauZener(int nstate, bool adiabatic, int points, std::mt19937& mt); LandauZener() = default;
LandauZener(int nstate, int points); LandauZener();

// function to perform the Landau-Zener jump
int jump(const Matrix& U, int state, int i, double tstep);
std::vector<std::tuple<int, double, bool>> jump(const Matrix& U, int state, int i, double tstep);

private:
// general variables necessary for the Landau-Zener algorithm
Matrix de, dz, ddz; std::vector<std::vector<int>> combs; bool adiabatic;
// matrix getters
const Matrix& getEd() const {return ed;} const Matrix& getDed() const {return ded;}

// random number generator with uniform distribution
std::mt19937 mt; std::uniform_real_distribution<double> dist;
private:
Matrix ed, ded; std::vector<std::vector<int>> combs;
};
45 changes: 23 additions & 22 deletions program/cdyn/src/landauzener.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#include "landauzener.h"

LandauZener::LandauZener() = default;

inline std::vector<std::vector<int>> Combinations(int n, int k) {
// create the bitmask that will get permuted and the resulting vector
std::string bitmask(k, 1); bitmask.resize(n, 0); std::vector<std::vector<int>> combs;
Expand All @@ -15,44 +17,43 @@ inline std::vector<std::vector<int>> Combinations(int n, int k) {
return combs;
}

LandauZener::LandauZener(int nstate, bool adiabatic, int points, std::mt19937& mt) : adiabatic(adiabatic), mt(mt), dist(0, 1) {
combs = Combinations(nstate, 2); de = Matrix::Zero(points, combs.size()), dz = Matrix::Zero(points, combs.size()), ddz = Matrix::Zero(points, combs.size());
LandauZener::LandauZener(int nstate, int points) {
combs = Combinations(nstate, 2); ed = Matrix::Zero(points, combs.size()), ded = Matrix::Zero(points, combs.size());
}

int LandauZener::jump(const Matrix& U, int state, int i, double tstep) {
std::vector<std::tuple<int, double, bool>> LandauZener::jump(const Matrix& U, int state, int i, double tstep) {
// loop over all the state combinations
for (size_t j = 0; j < combs.size(); j++) {

// calculate the state energy differences
de(i, j) = U(combs.at(j).at(1), combs.at(j).at(1)) - U(combs.at(j).at(0), combs.at(j).at(0));

// calculate the second derivative
if (i > 1) ddz(i, j) = (de(i, j) - 2 * de(i - 1, j) + de(i - 2, j)) / (tstep * tstep);
// if (i > 2) ddz(i, j) = (2 * de(i, j) - 5 * de(i - 1, j) + 4 * de(i - 2, j) - de(i - 3, j)) / (tstep * tstep);
// if (i > 3) ddz(i, j) = ((35.0/12) * de(i, j) - (26.0/3) * de(i - 1, j) + 9.5 * de(i - 2, j) - (14.0/3) * de(i - 3, j) + (11.0/12) * de(i - 4, j)) / (tstep * tstep);
ed(i, j) = U(combs.at(j).at(1), combs.at(j).at(1)) - U(combs.at(j).at(0), combs.at(j).at(0));

// calculate the first derivative
if (i > 0) dz(i, j) = (de(i, j) - de(i - 1, j)) / tstep;
// if (i > 1) dz(i, j) = (1.5 * de(i, j) - 2 * de(i - 1, j) + 0.5 * de(i - 2, j)) / tstep;
// if (i > 2) dz(i, j) = ((11.0/6) * de(i, j) - 3 * de(i - 1, j) + 1.5 * de(i - 2, j) - (1.0/3) * de(i - 3, j)) / tstep;
// calculate the first derivative of the energy difference
if (i > 0) ded(i, j) = (ed(i, j) - ed(i - 1, j)) / tstep;
}

// define the container for the accepted transitions
std::vector<std::tuple<int, double, bool>> transitions;

// loop over all the state combinations
for (size_t j = 0; j < combs.size() && i > 0; j++) {

// skip if the current state is not in the combination and define the probability
double P; if (state != combs.at(j).at(0) && state != combs.at(j).at(1)) continue;
if (state != combs.at(j).at(0) && state != combs.at(j).at(1)) continue;

// calculate the probability of state change
if (!adiabatic) P = 1 - std::exp(-2 * M_PI * std::pow(U(combs.at(j).at(0), combs.at(j).at(1)), 2) / std::abs(dz(i, j)));
else P = std::exp(-0.5 * M_PI * std::sqrt(std::pow(de(i, j), 3) / ddz(i, j)));
double P = 1 - std::exp(-2 * M_PI * std::pow(U(combs.at(j).at(0), combs.at(j).at(1)), 2) / std::abs(ded(i, j)));

// return the new state if the jump is made
if ((adiabatic && dz(i, j) * dz(i - 1, j) < 0 && dist(mt) < P) || (!adiabatic && de(i, j) * de(i - 1, j) < 0 && dist(mt) < P)) {
return state == combs.at(j).at(0) ? combs.at(j).at(1) : combs.at(j).at(0);
}
// add the transition to the container with the probability of the jump
transitions.push_back(std::make_tuple(state == combs.at(j).at(0) ? combs.at(j).at(1) : combs.at(j).at(0), std::isnan(P) ? 0 : P, false));

// if the trajectory is at a crossing point, enable the roll for the jump
if (ed(i, j) * ed(i - 1, j) < 0) std::get<2>(transitions.back()) = true;
}

// return the current state if no jump is made
return state;
// sort the transitions by probabilities from smallest to largest
std::sort(transitions.begin(), transitions.end(), [](const auto& a, const auto& b) {return std::get<1>(a) < std::get<1>(b);});

// return transitions
return transitions;
}
Loading

0 comments on commit 762626d

Please sign in to comment.