Skip to content

Commit 934e0f4

Browse files
alexvakimovalexvakimov
committed
Merge branch 'devel' of https://github.com/Quantum-Dynamics-Hub/libra-code into devel
2 parents ded6e5a + 79ce3fc commit 934e0f4

4 files changed

Lines changed: 248 additions & 44 deletions

File tree

‎src/dyn/heom/heom.cpp‎

Lines changed: 193 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#include "../../math_linalg/liblinalg.h"
2121
#include "../../math_specialfunctions/libspecialfunctions.h"
2222
#include "libheom.h"
23+
#include <unordered_map>
2324

2425
namespace liblibra{
2526
namespace libdyn{
@@ -28,6 +29,7 @@ namespace libheom{
2829
using namespace libutil;
2930
using namespace liblinalg;
3031
using namespace libspecialfunctions;
32+
using std::unordered_map;
3133

3234

3335

@@ -62,6 +64,178 @@ int compute_nn_tot(int d, int max_tier){
6264

6365

6466

67+
//vector< vector<int> > gen_next_level(vector< vector<int> >& parents){
68+
// /**
69+
// This function is similar to gen_next_level, except it takes a list of
70+
// parent vectors of integers and generates the list of the children for
71+
// all of the parent inputs.
72+
// */
73+
//
74+
// int num_parents = parents.size();
75+
// int num_children = parents[0].size(); // per parent
76+
//
77+
// vector< vector<int> > next_level;
78+
//
79+
// for(int i=0; i<num_parents; i++){
80+
// for(int j=0;j<num_children; j++){
81+
//
82+
// vector<int> child = parents[i];
83+
// child[j] += 1;
84+
//
85+
// if(! is_included(child, next_level) ){ next_level.push_back(child); }
86+
//
87+
// }// for j
88+
// }// for i
89+
//
90+
// return next_level;
91+
//}
92+
//
93+
//
94+
//
95+
//void gen_hierarchy(int d, int max_tier, int verbosity,
96+
// vector< vector<int> >& all_vectors,
97+
// vector< vector<int> >& vec_plus,
98+
// vector< vector<int> >& vec_minus){
99+
// /**
100+
//
101+
// d - is the size of the simplex, in HEOM it will be taken as nquant * (KK+1)
102+
//
103+
// max_tier - is the maximal tier of the vectors to be generated, the tier is
104+
// devined by the sum of the vector elements
105+
//
106+
// all_vectors - a list of the d-dimensional vectors of ints
107+
//
108+
// vec_plus : vec_plus[n][k] - index of the vector for which k-th element is +1 of that of the n-th vector
109+
//
110+
// vec_minus : vec_minus[n][k] - index of the vector for which k-th element is -1 of that of the n-th vector
111+
//
112+
// The vec_minus and vec_plus contain -1 for the situation when there is no such a vector included in
113+
// the current hierarchy structure
114+
//
115+
// */
116+
//
117+
// int i, j, k;
118+
//
119+
//
120+
// vector< vector<int> > all_coordinates; // for each vector - the coordinates are (L, i)
121+
// vector<int> tier_nums; // the number of the nodes for the tier up to a given one
122+
//
123+
// vector< vector<int> > parents(1, vector<int>(d, 0));
124+
//
125+
// for(int tier=0; tier<=max_tier; tier++){
126+
//
127+
// int iparent = 0;
128+
// int num_parents = parents.size();
129+
//
130+
// for(i=0; i<num_parents; i++){
131+
// all_vectors.push_back( parents[i] );
132+
// vector<int> coord(2,0);
133+
// coord[0] = tier;
134+
// coord[1] = iparent;
135+
// all_coordinates.push_back( coord );
136+
// iparent += 1;
137+
//
138+
// }// for i
139+
//
140+
// tier_nums.push_back( all_vectors.size() ) ;
141+
//
142+
// vector< vector<int> > new_parents = gen_next_level(parents);
143+
// parents.resize(new_parents.size());
144+
// parents = new_parents;
145+
//
146+
// }// for tier
147+
//
148+
// //#==============================================
149+
//
150+
//
151+
// int num_nodes = all_vectors.size();
152+
// vec_plus = vector< vector<int> >(num_nodes, vector<int>(d, -1));
153+
// vec_minus = vector< vector<int> >(num_nodes, vector<int>(d, -1));
154+
//
155+
//
156+
// for(i=0; i<num_nodes; i++){
157+
// int L = all_coordinates[i][0];
158+
// int ipos = all_coordinates[i][1];
159+
//
160+
// for(k=0; k<d; k++){
161+
// vector<int> np = all_vectors[i];
162+
// np[k] += 1;
163+
//
164+
// int max_range = max_tier;
165+
// if(L<max_tier){
166+
// max_range = tier_nums[L+1];
167+
// }
168+
//
169+
// for(j = tier_nums[L]; j < max_range; j++){
170+
// if(np == all_vectors[j]){
171+
// vec_plus[i][k] = j;
172+
// }
173+
// }
174+
// }// for k
175+
//
176+
// for(k=0; k<d; k++){
177+
// vector<int> nm = all_vectors[i];
178+
// nm[k] -= 1;
179+
//
180+
// int min_range = 0;
181+
// if(L>=2){
182+
// min_range = tier_nums[L-2];
183+
// }
184+
// for(j=min_range; j<tier_nums[L-1]; j++){
185+
// if(nm == all_vectors[j]){
186+
// vec_minus[i][k] = j;
187+
// }
188+
// }
189+
// }// for k
190+
// }// for i
191+
//
192+
//
193+
// if(verbosity>0){
194+
// cout<<"Number of nodes = "<<num_nodes<<endl;
195+
// cout<<"Tier nums: \n";
196+
// for(i=0; i<tier_nums.size(); i++){
197+
// cout<<"Tier "<<i<<" nums = "<<tier_nums[i]<<endl;
198+
// }
199+
// }
200+
//
201+
// if(verbosity>1){
202+
// for(i=0; i<num_nodes; i++){
203+
// cout<<"index= "<<i<<" ";
204+
// cout<<" vector=["; for(j=0;j<d;j++){ cout<<all_vectors[i][j]<<","; } cout<<"] ";
205+
// cout<<" coordinates=["<<all_coordinates[i][0]<<","<<all_coordinates[i][1]<<"] ";
206+
// cout<<" vec_plus=["; for(j=0;j<d;j++){ cout<<vec_plus[i][j]<<","; } cout<<"] ";
207+
// cout<<" vec_minus=["; for(j=0;j<d;j++){ cout<<vec_minus[i][j]<<","; } cout<<"] ";
208+
// cout<<"\n";
209+
// }
210+
// }
211+
//
212+
//}
213+
214+
215+
struct VectorHash {
216+
size_t operator()(const vector<int>& v) const {
217+
size_t seed = v.size();
218+
for(auto& i : v){
219+
seed ^= i + 0x9e3779b9 + (seed<<6) + (seed>>2);
220+
}
221+
return seed;
222+
}
223+
};
224+
225+
int count_trailing_zeros(const vector<int>& v) {
226+
227+
int count = 0;
228+
229+
for(size_t i = v.size(); i-- > 0; ){
230+
if(v[i] == 0)
231+
count++;
232+
else
233+
break;
234+
}
235+
236+
return count;
237+
}
238+
65239
vector< vector<int> > gen_next_level(vector< vector<int> >& parents){
66240
/**
67241
This function is similar to gen_next_level, except it takes a list of
@@ -71,16 +245,20 @@ vector< vector<int> > gen_next_level(vector< vector<int> >& parents){
71245

72246
int num_parents = parents.size();
73247
int num_children = parents[0].size(); // per parent
248+
int num_start = 0; // starting index to begin adding child excitations
74249

75250
vector< vector<int> > next_level;
76251

77252
for(int i=0; i<num_parents; i++){
78-
for(int j=0;j<num_children; j++){
79-
253+
if (i != 0){
254+
num_start = num_children - count_trailing_zeros(parents[i]) - 1;
255+
}
256+
for(int j=num_start;j<num_children; j++){
80257
vector<int> child = parents[i];
258+
81259
child[j] += 1;
82260

83-
if(! is_included(child, next_level) ){ next_level.push_back(child); }
261+
next_level.push_back(child);
84262

85263
}// for j
86264
}// for i
@@ -89,7 +267,6 @@ vector< vector<int> > gen_next_level(vector< vector<int> >& parents){
89267
}
90268

91269

92-
93270
void gen_hierarchy(int d, int max_tier, int verbosity,
94271
vector< vector<int> >& all_vectors,
95272
vector< vector<int> >& vec_plus,
@@ -147,9 +324,12 @@ void gen_hierarchy(int d, int max_tier, int verbosity,
147324

148325

149326
int num_nodes = all_vectors.size();
150-
vec_plus = vector< vector<int> >(num_nodes, vector<int>(d, -1));
151-
vec_minus = vector< vector<int> >(num_nodes, vector<int>(d, -1));
152-
327+
vec_plus = vector< vector<int> >(num_nodes, vector<int>(d, -1));
328+
vec_minus = vector< vector<int> >(num_nodes, vector<int>(d, -1));
329+
unordered_map<vector<int>, int, VectorHash> vector_to_index;
330+
for(int i=0; i<num_nodes; ++i){
331+
vector_to_index[all_vectors[i]] = i;
332+
}
153333

154334
for(i=0; i<num_nodes; i++){
155335
int L = all_coordinates[i][0];
@@ -158,31 +338,16 @@ void gen_hierarchy(int d, int max_tier, int verbosity,
158338
for(k=0; k<d; k++){
159339
vector<int> np = all_vectors[i];
160340
np[k] += 1;
161-
162-
int max_range = max_tier;
163-
if(L<max_tier){
164-
max_range = tier_nums[L+1];
341+
auto itp = vector_to_index.find(np);
342+
if(itp != vector_to_index.end()){
343+
vec_plus[i][k] = itp->second;
165344
}
166345

167-
for(j = tier_nums[L]; j < max_range; j++){
168-
if(np == all_vectors[j]){
169-
vec_plus[i][k] = j;
170-
}
171-
}
172-
}// for k
173-
174-
for(k=0; k<d; k++){
175346
vector<int> nm = all_vectors[i];
176347
nm[k] -= 1;
177-
178-
int min_range = 0;
179-
if(L>=2){
180-
min_range = tier_nums[L-2];
181-
}
182-
for(j=min_range; j<tier_nums[L-1]; j++){
183-
if(nm == all_vectors[j]){
184-
vec_minus[i][k] = j;
185-
}
348+
auto itm = vector_to_index.find(nm);
349+
if(itm != vector_to_index.end()){
350+
vec_minus[i][k] = itm->second;
186351
}
187352
}// for k
188353
}// for i

‎src/dyn/heom/heom.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,8 @@ namespace libheom{
4747
/// General hierarchy setup
4848
int compute_nn_tot(int d, int max_tier);
4949

50+
int count_trailing_zeros(const vector<int>& v);
51+
5052
vector< vector<int> > gen_next_level(vector< vector<int> >& parents);
5153

5254
void gen_hierarchy(int d, int max_tier, int verbosity,

‎src/libra_py/dynamics/heom/compute.py‎

Lines changed: 21 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -144,7 +144,8 @@ def transform_adm(rho, rho_scaled, aux_memory, params, direction):
144144
pack_mtx(aux_memory["rho_unpacked"], rho)
145145

146146

147-
def run_dynamics(dyn_params, Ham, rho_init):
147+
#def run_dynamics(dyn_params, Ham, rho_init):
148+
def run_dynamics(dyn_params, Ham, rho_init, adm_init=None):
148149
"""
149150
This functions integrates the HEOM for a given system's Hamiltonian, initial conditions, and bath parameters
150151
@@ -188,6 +189,9 @@ def run_dynamics(dyn_params, Ham, rho_init):
188189
* **dyn_params["nsteps"]** ( int )
189190
How many steps of the dynamics to perform [ default: 10 ]
190191
192+
* **dyn_params["nprint"]** ( int )
193+
The number of steps in between data saves [ default: 1 ]
194+
191195
* **dyn_params["verbosity"]** ( int )
192196
The level of the run-time printout of any useful information [ default: -1 ]
193197
@@ -317,7 +321,7 @@ def run_dynamics(dyn_params, Ham, rho_init):
317321
"temperature": 300.0,
318322
"el_phon_couplings": initialize_el_phonon_couplings(Ham.num_of_cols),
319323

320-
"dt": 0.1 * units.fs2au, "nsteps": 10,
324+
"dt": 0.1 * units.fs2au, "nsteps": 10, "nprint": 1,
321325
"verbosity": -1, "progress_frequency": 0.1,
322326

323327
"truncation_scheme": 1, "do_scale": 0,
@@ -334,6 +338,7 @@ def run_dynamics(dyn_params, Ham, rho_init):
334338
comn.check_input(params, default_params, critical_params)
335339

336340
nsteps = params["nsteps"]
341+
#nprint = params["nprint"]
337342
print_freq = int(params["progress_frequency"] * nsteps)
338343

339344
# ============= System ======================
@@ -389,9 +394,16 @@ def run_dynamics(dyn_params, Ham, rho_init):
389394
aux_memory["drho_unpacked_scaled"].append(CMATRIX(nquant, nquant))
390395

391396
# Initial conditions
392-
x_ = Py2Cpp_int(list(range(nquant)))
393-
y_ = Py2Cpp_int(list(range(nquant)))
394-
push_submatrix(rho, rho_init, x_, y_)
397+
if adm_init is None:
398+
x_ = Py2Cpp_int(list(range(nquant)))
399+
y_ = Py2Cpp_int(list(range(nquant)))
400+
push_submatrix(rho, rho_init, x_, y_)
401+
else:
402+
# should test to make sure that len(adm_init) == nn_tot-1:
403+
aux_memory["rho_unpacked"][0] = rho_init
404+
for n in range(1,nn_tot):
405+
aux_memory["rho_unpacked"][n] = adm_init[n-1]
406+
pack_mtx(aux_memory["rho_unpacked"], rho)
395407

396408
# unpack_mtx(aux_memory["rho_unpacked"], rho)
397409

@@ -423,7 +435,8 @@ def run_dynamics(dyn_params, Ham, rho_init):
423435
aux_print_matrices(0, aux_memory["rho_unpacked_scaled"])
424436

425437
# Initialize savers
426-
_savers = save.init_heom_savers(params, nquant)
438+
#_savers = save.init_heom_savers(params, nquant)
439+
_savers = save.init_heom_savers(params, nquant, nn_tot)
427440

428441
# ============== Propagation =============
429442

@@ -460,7 +473,9 @@ def run_dynamics(dyn_params, Ham, rho_init):
460473
update_filters(rho_scaled, params, aux_memory)
461474

462475
# ================= Propagation for one timestep ==================================
476+
#print("MADE IT TO RK4 STEP") # These lines were not here, remove
463477
rho_scaled = RK4(rho_scaled, params["dt"], compute_heom_derivatives, params)
478+
#print("PASSED AN RK4 STEP") # These lines were not here, remove
464479

465480
end = time.time()
466481
print(F"Calculations took {end - start} seconds")

0 commit comments

Comments
 (0)