Repository navigation
Expand file tree
/
Copy pathanalysis.py
More file actions
177 lines (132 loc) · 6.49 KB
/
Copy pathanalysis.py
File metadata and controls
177 lines (132 loc) · 6.49 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
import os
import sys
import math
import time
import numpy as np
import h5py
from liblibra_core import *
import util.libutil as comn
from libra_py import units as units
from libra_py import data_conv, data_stat, data_outs, data_read, data_savers
import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import adi_integrators.recipes
def plots(prefix, imdl, X, Y):
"""
This function produces nice plots of the ground state populations for all calculations
"""
plt.rc('axes', titlesize=12) # fontsize of the axes title
plt.rc('axes', labelsize=12) # fontsize of the x and y labels
plt.rc('legend', fontsize=10) # legend fontsize
plt.rc('xtick', labelsize=8) # fontsize of the tick labels
plt.rc('ytick', labelsize=8) # fontsize of the tick labels
plt.rc('figure.subplot', left=0.2)
plt.rc('figure.subplot', right=0.95)
plt.rc('figure.subplot', bottom=0.13)
plt.rc('figure.subplot', top=0.88)
figure = plt.figure(num=imdl, figsize=(3.21, 2.41), dpi=300, edgecolor='black', frameon=True)
plt.plot(X[imdl, :], Y[imdl, :, 0, 0], color="black", label="reflection on state 0", linewidth=1)
plt.plot(X[imdl, :], Y[imdl, :, 0, 1], color="blue", label="transmission on state 0", linewidth=1)
plt.plot(X[imdl, :], Y[imdl, :, 1, 0], color="red", label="reflection on state 1", linewidth=1)
plt.plot(X[imdl, :], Y[imdl, :, 1, 1], color="green", label="transmission on state 1", linewidth=1)
plt.title("",fontsize=9.5)
plt.legend(fontsize=6.75, ncol=1, loc='upper left')
plt.xlabel('Momentum, a.u.',fontsize=10)
plt.ylabel('Population ',fontsize=10)
plt.tight_layout()
plt.savefig(F"{prefix}.png", dpi=300)
plt.show()
def get_stats(prefix, mdls, P0, batches, ntraj, methods):
nmdls = len(mdls)
npts = len(P0)
nstates = 2
nscatt = 2 # scattering channels
X = np.zeros((nmdls, npts), dtype=float)
Y = np.zeros((nmdls, npts, nstates, nscatt), dtype=float)
cols = list(range(ntraj))
for method in methods:
recipe = recipes.recipe_inv_mapping(method)
rec_name, rec_params = recipes.set_recipe(recipe)
for imdl, mdl in enumerate(mdls):
for ip0, p0 in enumerate(P0):
X[imdl, ip0] = p0
cnt = 0.0
for batch in batches:
#filename = F"{prefix}{mdl}_{int(p0)}/start_s0_BC_FSSH_batch{batch}"
filename = F"{prefix}{mdl}_P0_{p0}_RECIPE_{method}/start_s0_{rec_name}_batch{batch}"
states = data_read.get_data_from_file2(F"{filename}/states.txt", cols)
coords = data_read.get_data_from_file2(F"{filename}/q.txt", cols)
for itraj in range(ntraj):
ist = int(states[itraj][-1])
q = coords[itraj][-1]
if q<-5: # reflection on state ist
Y[imdl, ip0, ist, 0 ] += 1
cnt += 1
elif q>5: # transmission on state ist
Y[imdl, ip0, ist, 1 ] += 1
cnt += 1
print(Y, cnt)
if cnt > 0:
Y[imdl, ip0, :, :] = Y[imdl, ip0, :, :] / cnt
#plots(F"model_{mdl}", imdl, X, Y)
with h5py.File(F"model_{mdl}_method_{method}.hdf", "w") as f:
#g = f.create_group("data")
f.create_dataset("X", data = X)
f.create_dataset("Y", data = Y)
def plot(mdls, methods):
for method in methods:
for imdl, mdl in enumerate(mdls):
with h5py.File(F"model_{mdl}_method_{method}.hdf", "r") as f:
print( list(f.keys() ) )
print(f["X"].shape)
print(f["Y"].shape)
plots(F"model_{mdl}_method_{method}", imdl, f["X"], f["Y"])
def compute_error(x1, model2_path, prop='D_adi', method=0, min_nsteps=1, istate_model1=0, istate_model2=0):
"""
This funtion computes the accumulation error between two set of coefficients.
Args:
x1 (numpy array): The density of the reference model1
model2_path (string): The full path to model 2 hdf file data.
method (integer): This argument can take two values of 0 and 1.
0: The minimum number of steps for time-grid is defined based on the
shape of model 2 and 1 length.
1: The minimum number of steps for time-grid is defined based on the
min_steps by the user.
min_nsteps (integer): The minimum number of steps for the time-grid in case method=1 is selected.
istate_model1 (integer): The initial state of model 1.
istate_model2 (integer): The initial state of model 2.
Returns:
x1 (numpy array): The reference density matrix.
x1 (numpy array): The second density matrix.
error (float): The average value of error computed over a set of time-grid points.
"""
F2 = h5py.File(model2_path)
x2 = np.array(F2[F'{prop}/data'])
nsteps1 = x1.shape[0]
nsteps2 = x2.shape[0]
if nsteps1 >= nsteps2 and nsteps1%nsteps2==0:
if method==1:
coeff2 = int(nsteps2/min_nsteps)
nsteps2 = min_nsteps
else:
coeff2 = 1
coeff1 = int(nsteps1/nsteps2)
diff = x1[::coeff1]-x2[::coeff2]
nsteps = nsteps2
elif nsteps1 < nsteps2 and nsteps2%nsteps1==0:
if method==1:
coeff1 = int(nsteps1/min_nsteps)
nsteps1 = min_nsteps
else:
coeff1 = 1
coeff2 = int(nsteps2/nsteps1)
diff = x1[::coeff1]-x2[::coeff2]
nsteps = nsteps1
else:
print('Error! nsteps1 and nsteps2 are not divisible...')
# Needs to sys.exit(0) as well but I do not do this in Jupyter notebook :)
diff_squared = np.power(np.abs(diff),2)
error = np.sum(diff_squared)/nsteps
F2.close()
return x1, x2, error