In [3]:
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
#dataDict = np.load('datainput_FinerGrid.npy')[()]
dataDict = np.load('/home/alessio/n-Propagation/datanewoneWithNACnow.npy',allow_pickle=True)[()]
qp.printDictKeys(dataDict)
print(dataDict['kinCube'].shape, dataDict['potCube'].shape, dataDict['dipCUBE'].shape, dataDict['geoCUBE'].shape, dataDict['nacCUBE'].shape)
%matplotlib notebook
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from ipywidgets import interact,fixed #, interactive, fixed, interact_manual
import ipywidgets as widgets
In [4]:
gsm_phi_ind = dataDict['phis'].index('P000-000')
gsm_gam_ind = dataDict['gams'].index('P016-923')
gsm_the_ind = dataDict['thes'].index('P114-804')
gsm_phi_ind, gsm_gam_ind, gsm_the_ind
Out[4]:
In [5]:
import quantumpropagator as qp
from quantumpropagator import labTranformA,fromLabelsToFloats
phiV, gamV, theV = fromLabelsToFloats(dataDict)
phis = labTranformA(dataDict['phis'])
gams = labTranformA(dataDict['gams'])
thes = labTranformA(dataDict['thes'])
phiL = len(phis)
gamL = len(gams)
theL = len(thes)
In [6]:
phis,phiV,phiV[0]-phiV[-1],phiV[0]-phiV[1],phiL, np.pi/((phiV[1]-phiV[0])/4)
Out[6]:
In [7]:
gams,gamV,gamV[0]-gamV[-1],gamV[0]-gamV[1],gamL
Out[7]:
In [8]:
thes,theV,theV[0]-theV[-1],theV[0]-theV[1], theL
Out[8]:
In [7]:
def do3dplot(xs,ys,zss):
'with mesh function'
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X,Y = np.meshgrid(ys,xs)
ax.view_init(elev=1., azim=0.)
ax.plot_wireframe(X, Y, zss)
#ax.scatter(X, Y, zss)
def do3dplot2(X,Y,Z):
'without mesh function'
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
#ax.plot_wireframe(X, Y, Z)
ax.scatter(X, Y, Z)
In [8]:
def kinElem(mkelem,deriv,phi):
expl = '{}\'th order derivative coefficient in K matrix'.format(deriv)
tis = ['Tpp','Tpg','Tpt','Tgp','Tgg','Tgt','Ttp','Ttg','Ttt']
tisexpl = '{} -> '.format(tis[mkelem])
print(tisexpl + expl)
x = gams
y = thes
z = dataDict['kinCube'][phi,:,:,mkelem,deriv]
do3dplot(x,y,z)
interact(kinElem, mkelem = widgets.IntSlider(min=0,max=8,step=1,value=0), deriv = widgets.IntSlider(min=0,max=2,step=1,value=0),phi = widgets.IntSlider(min=0,max=phiL-1,step=1,value=0));
In [9]:
dataDict['kinCube'][11,10,14],dataDict['kinCube'][8,10,14]
Out[9]:
In [9]:
import bottleneck as bn
a = nac[:,:,:,4,5,7:,:].flatten()
print(np.amax(a))
biggest = np.sort(bn.partition(a, a.size-30)[-30:])
string = [ '{}'.format(x) for x in biggest ]
print(string)
In [10]:
nac = dataDict['nacCUBE']
print(nac.shape)
cutter = nac[15,:,:,4,5,1,1]
#print(cutter)
def do3dplot3(xs,ys,zss):
'with mesh function'
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X,Y = np.meshgrid(ys,xs)
#ax.plot_wireframe(X, Y, zss)
ax.set_zlim(-1, 1)
ax.scatter(X, Y, zss)
do3dplot3(gams,thes,cutter)
In [11]:
potential = dataDict['potCube']- np.amin(dataDict['potCube'])
In [12]:
def expARR(array,howmany,direction=None):
direction = direction or 'both'
dx = array[1]-array[0]
after = (array[-1] + (dx * np.arange(howmany+1)))[1:]
before = np.flip((array[0] - (dx * np.arange(howmany+1)))[1:],0)
if direction == 'both':
return np.concatenate((before,array,after))
elif direction == 'dx':
return np.concatenate((array,after))
elif direction == 'sx':
return np.concatenate((before,array))
def doubleAxespoins(Y):
N = len(Y)
X = np.arange(0, 2*N, 2)
X_new = np.arange(2*N-1) # Where you want to interpolate
Y_new = np.interp(X_new, X, Y)
return(Y_new)
In [13]:
def pot2D(potential,state,kind,dim):
print('Doing cut in {} with value ({:8.4f},{:8.4f})'.format(kind,dimV[kind][dim],dims[kind][dim]))
if kind == 'Phi':
pot = potential[dim,:,:,state]
do3dplot(gams,thes,pot)
elif kind == 'Gam':
pot = potential[:,dim,:,state]
do3dplot(phis,thes,pot)
elif kind == 'The':
pot = potential[:,:,dim,state]
do3dplot(phis,gams,pot)
dimV = { 'Phi': phiV, 'Gam': gamV, 'The': theV } # real values
dims = { 'Phi': phis, 'Gam': gams, 'The': thes } # for labels
state = 0
kinds = ['Phi','Gam','The']
def fun_pot2D(kind,state,dim):
pot2D(potential, state,kind, dim)
def nested(kinds):
dimensionV = dimV[kinds]
interact(fun_pot2D, kind=fixed(kinds), state = widgets.IntSlider(min=0,max=7,step=1,value=0), dim = widgets.IntSlider(min=0,max=(len(dimensionV)-1),step=1,value=0))
interact(nested, kinds = ['Phi','Gam','The']);
In [21]:
line = potential[14,14,21]
line[6]-line[1]
Out[21]:
In [15]:
def polinomial_fitting(true_lines,fit_lines,both_cut, phi_cut, gam_cut, expand_N, degree_fit):
nstates = 8
cut_sx, cut_dx = both_cut
fig = plt.figure(figsize=(15,10))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'mediumpurple']
y = potential[phi_cut,gam_cut,cut_sx:cut_dx,:]
x = thes[cut_sx:cut_dx]
tit = 'Phi = {:8.4f} | Gamma = {:8.4f}'.format(phis[phi_cut], gams[gam_cut])
plt.title(tit)
if true_lines:
for iii in range(nstates):
plt.plot(x, y[:,iii], ls='-', lw=1.5, color=colors[iii]);
#plt.plot(x,y, ls='-', lw=1)# , marker = 'o',markersize = '1');
coefficients = np.apply_along_axis(np.polyfit,0,x,y,degree_fit)
ncoe, states = coefficients.shape
new_x = expARR(x,expand_N)
functs = np.array([ np.poly1d(coefficients[:,i])(new_x) for i in range(states)])
new_y = functs.T
if fit_lines:
for iii in range(nstates):
plt.plot(new_x, new_y[:,iii], ls='--', lw=1, color=colors[iii]);
interact(polinomial_fitting,
true_lines = widgets.Checkbox(value=True, description='True values'),
fit_lines = widgets.Checkbox(value=True, description='Fit values'),
both_cut = widgets.IntRangeSlider(min=0,max=theL-1,step=1,value=[0, theL-1]),
#cut_sx = widgets.IntSlider(min=0,max=theL-1,step=1,value=0),
#cut_dx = widgets.IntSlider(min=0,max=theL-1,step=1,value=theL),
phi_cut = widgets.IntSlider(min=0,max=phiL-1,step=1,value=16),
gam_cut = widgets.IntSlider(min=0,max=gamL-1,step=1,value=20),
expand_N = widgets.IntSlider(min=0,max=40,step=1,value=0),
degree_fit = widgets.IntSlider(min=0,max=10,step=1,value=2));
In [15]:
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
ground = potential[:,:,:,1]
tupl = ground.shape
ground.shape,phiV.shape,gamV.shape,theV.shape
phi_mesh, gam_mesh, the_mesh = np.meshgrid(phiV,gamV,theV, indexing='ij')
phi_prime = phi_mesh.reshape(-1,1)
gam_prime = gam_mesh.reshape(-1,1)
the_prime = the_mesh.reshape(-1,1)
X = np.concatenate((phi_prime,gam_prime,the_prime), axis=1)
y = ground.reshape(-1,1)
poly = PolynomialFeatures(degree=10)
X_ = poly.fit_transform(X)
clf = LinearRegression()
clf.fit(X_, y)
prediction = clf.predict(X_)
prediction_reshaped = prediction.reshape(tupl)
corto = 0
def two_d_see_it(phi_cut):
ground_corto = ground[phi_cut,:,:]
pred_corto = prediction_reshaped[phi_cut,:,:]
gam_mesh_corto = gam_mesh[phi_cut,:,:]
the_mesh_corto = the_mesh[phi_cut,:,:]
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(gam_mesh_corto, the_mesh_corto, pred_corto, s=3, alpha=.5, color='k', marker='o')
ax.scatter(gam_mesh_corto, the_mesh_corto, ground_corto, s=3, alpha=.5, color='r', marker='o')
#ax.scatter(X_mesh_ext, Y_mesh_ext, prediction_mesh_ext, s=3, alpha=.5, color='r', marker='o')
fig.tight_layout()
interact(two_d_see_it,
phi_cut = widgets.IntSlider(min=0,max=phiL-1,step=1,value=0));
In [15]:
import scipy.io as sio
saveNew_matlab_Vector = False
phi_L, gam_L, the_L, nstates = potential.shape
if saveNew_matlab_Vector:
allstates = {}
for s in range(nstates):
ground = potential[:,:,:,s]
tupl = ground.shape
ground.shape,phiV.shape,gamV.shape,theV.shape
phi_mesh, gam_mesh, the_mesh = np.meshgrid(phiV,gamV,theV, indexing='ij')
phi_prime = phi_mesh.reshape(-1,1)
gam_prime = gam_mesh.reshape(-1,1)
the_prime = the_mesh.reshape(-1,1)
X = np.concatenate((phi_prime,gam_prime,the_prime), axis=1)
y = ground.reshape(-1,1)
poly = PolynomialFeatures(degree=10)
X_ = poly.fit_transform(X)
clf = LinearRegression()
clf.fit(X_, y)
prediction = clf.predict(X_)
prediction_reshaped = prediction.reshape(tupl)
stringName = 'TDM{}'.format(s)
allstates[stringName] = prediction_reshaped
sio.savemat('np_vector.mat', allstates)
In [16]:
aaa = sio.loadmat('/home/alessio/n-Propagation/tdm.mat')
new = aaa['v0']
new_phi_L, new_gam_L, new_the_L = new.shape
ext_potential = np.empty((new_phi_L, new_gam_L, new_the_L, nstates))
phi_cayo = expARR(phiV,15,'both')
gam_cayo = expARR(gamV,15,'both')
the_cayo = expARR(theV,30,'both')
phi_mesh_cayo, gam_mesh_cayo, the_mesh_cayo = np.meshgrid(phi_cayo,gam_cayo,the_cayo, indexing='ij')
print(nstates)
for s in range(nstates):
label = 'v{}'.format(s)
ext_potential[:,:,:,s] = aaa[label]
def two_d_see_it(phi_cut,state):
corto = phi_cut
ground_corto = potential[corto,:,:,state]
pred_corto_cayo = ext_potential[corto+15,:,:,state]
gam_mesh_corto = gam_mesh[corto,:,:]
the_mesh_corto = the_mesh[corto,:,:]
gam_mesh_corto_cayo = gam_mesh_cayo[corto+15,:,:]
the_mesh_corto_cayo = the_mesh_cayo[corto+15,:,:]
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(gam_mesh_corto_cayo, the_mesh_corto_cayo, pred_corto_cayo, s=3, alpha=.2, color='k', marker='o')
ax.scatter(gam_mesh_corto, the_mesh_corto, ground_corto, s=3, alpha=.5, color='r', marker='o')
#ax.scatter(X_mesh_ext, Y_mesh_ext, prediction_mesh_ext, s=3, alpha=.5, color='r', marker='o')
fig.tight_layout()
interact(two_d_see_it,
phi_cut = widgets.IntSlider(min=0,max=phiL-1,step=1,value=0),
state = widgets.IntSlider(min=0,max=7,step=1,value=0));
In [17]:
ppL, ggL, ttL = new.shape
phis_n = phi_cayo
gams_n = gam_cayo
thes_n = the_cayo
pp_cayo,gg_cayo,tt_cayo = qp.fromFloatsToLabels(phi_cayo, gam_cayo, the_cayo)
In [18]:
a = np.arange(6) + 0.5
b = a
In [ ]:
import pickle
doublALL = np.zeros((ppL,ggL,ttL))
doublALL2 = np.zeros((ppL,ggL,ttL,8))
kineti = np.zeros((ppL,ggL,ttL,9,3))
smo = np.zeros((ppL,ggL,ttL,8,8,3))
geom = np.zeros((ppL,ggL,ttL,15,3))
dipo = np.zeros((ppL,ggL,ttL,3,8,8))
#nac = np.zeros((ppL,ggL,ttL,8,8,15,3))
nacINSIDE = np.pad(dataDict['nacCUBE'], [(15, 15), (15, 15), (30, 30), (0, 0), (0, 0), (0, 0), (0, 0)], mode='constant')
umass = 1836
C_mass = 12 * umass
H_mass = 1 * umass
# Toggle this
saveNew = False
saveFile = False
if saveNew:
for pi, p in qp.log_progress(enumerate(phis_n),every=1,size=(len(phis_n))):
for gi, g in enumerate(gams_n):
for ti,t in enumerate(thes_n):
kineti[pi,gi,ti] = qp.calc_g_G(p,g,t)
smo_element = qp.calc_s_mat(p,g,t)
geom[pi,gi,ti] = qp.generateNorbGeometry(p,g,t,True) # True is the "vector" mode that does not save any file
for i in range(8):
for j in range(8):
thisnacelement = nacINSIDE[pi,gi,ti,i,j] # I want nac, this matrix is already antisymmetrized.
# I need nac to be divided by mass
thisnacelement_linear = np.concatenate(thisnacelement[7:],axis=0) # nac put in the same order as s matrix
massVector = np.array([C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass])
final_scaled_vector_tau_prime = thisnacelement_linear/massVector
smo[pi,gi,ti,i,j] = smo_element@final_scaled_vector_tau_prime
newdict = {}
newdict['smoCube'] = smo
newdict['kinCube'] = kineti
newdict['potCube'] = ext_potential
newdict['dipCUBE'] = np.pad(dataDict['dipCUBE'], [(15, 15), (15, 15), (30, 30), (0, 0), (0, 0), (0, 0)], mode='constant')
newdict['geoCUBE'] = geom
newdict['phis'] = pp_cayo
newdict['gams'] = gg_cayo
newdict['thes'] = tt_cayo
if saveFile:
with open('newExtrapolated_allCorrectionZZZ.pickle', 'wb') as f:
pickle.dump(newdict, f, pickle.HIGHEST_PROTOCOL)
#np.save('newExtrapolated_allCorrection.npy', newdict)
qp.printDict(newdict)
In [19]:
umass = 1836
C_mass = 12 * umass
H_mass = 1 * umass
a = (5,5,10)
pi = 20
gi = 20
ti = 40
i = 0
j = 1
p = phis_n[pi]
g = gams_n[gi]
t = thes_n[ti]
s = qp.calc_s_mat(p,g,t,True) #verbose
thisnacelement = nacINSIDE[pi,gi,ti,i,j] # I want nac
thisnacelement_linear = np.concatenate(thisnacelement[7:],axis=0) # nac put in the smae order as s matrix
massVector = np.array([C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,C_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass,H_mass])
thisnacelement_linear/massVector
In [23]:
a = newdict['smoCube'][15:-15,15:-15,30:-30,0,1,0]
print(a.shape)
a.tofile('porchetta.txt')
do3dplot(gams,thes,np.abs(a[3]))
#do3dplot(gams,thes,a)
In [ ]:
print(smo.shape)
ah = smo[15:-15,15:-15,30:-30,:,:,:]
def two_NAC_see_it(phi_cut, save=None):
save = save or False
cut_phi = ah[phi_cut,:,:,:,:,:]
X,Y = np.meshgrid(thes,gams)
maxi=3
mini=-3
fig = plt.figure(figsize=(40,15))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'mediumpurple']
all_states_indexes = np.arange(nstates)
for i in all_states_indexes:
others = np.delete(all_states_indexes, [i], axis=0)
ax1 = fig.add_subplot(3,8,i+1, projection='3d')
ax1.set_title('Phi: {} | State {} along X'.format(phi_cut,i))
ax1.set_xlabel('Theta')
ax1.set_ylabel('Gamma')
ax1.set_zlabel('DipoleValue')
#ax1.set_zlim(mini,maxi)
ax2 = fig.add_subplot(3,8,i+9, projection='3d')
ax2.set_title('State {} along Y'.format(i))
ax2.set_xlabel('Theta')
ax2.set_ylabel('Gamma')
ax2.set_zlabel('DipoleValue')
#ax2.set_zlim(mini,maxi)
ax3 = fig.add_subplot(3,8,i+17, projection='3d')
ax3.set_title('State {} along Z'.format(i))
ax3.set_xlabel('Theta')
ax3.set_ylabel('Gamma')
ax3.set_zlabel('DipoleValue')
#ax3.set_zlim(mini,maxi)
for j in others:
cutfirst = cut_phi[:,:,i,j,:]
cutx = cutfirst[:,:,0]
cuty = cutfirst[:,:,1]
cutz = cutfirst[:,:,2]
col = colors[j]
ax1.plot_wireframe(X,Y,cutx, alpha=.5, color=col)
ax2.plot_wireframe(X,Y,cuty, alpha=.5, color=col)
ax3.plot_wireframe(X,Y,cutz, alpha=.5, color=col)
fig.tight_layout()
if save:
name_fig = 'NACMATRIX_phi{:02d}.png'.format(phi_cut)
fig.savefig(name_fig)
savefigures = False
if savefigures:
for phi_cut in range(phiL):
two_NAC_see_it(phi_cut,True)
In [ ]:
dipo.shape
In [ ]:
dipo = dataDict['dipCUBE']
def two_d_see_it(phi_cut, save=None):
save = save or False
cut_phi = dipo[phi_cut,:,:,:,:,:]
X,Y = np.meshgrid(thes,gams)
maxi=3
mini=-3
fig = plt.figure(figsize=(40,15))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'mediumpurple']
all_states_indexes = np.arange(nstates)
for i in all_states_indexes:
others = np.delete(all_states_indexes, [i], axis=0)
ax1 = fig.add_subplot(3,8,i+1, projection='3d')
ax1.set_title('Phi: {} | State {} along X'.format(phi_cut,i))
ax1.set_xlabel('Theta')
ax1.set_ylabel('Gamma')
ax1.set_zlabel('DipoleValue')
ax1.set_zlim(mini,maxi)
ax2 = fig.add_subplot(3,8,i+9, projection='3d')
ax2.set_title('State {} along Y'.format(i))
ax2.set_xlabel('Theta')
ax2.set_ylabel('Gamma')
ax2.set_zlabel('DipoleValue')
ax2.set_zlim(mini,maxi)
ax3 = fig.add_subplot(3,8,i+17, projection='3d')
ax3.set_title('State {} along Z'.format(i))
ax3.set_xlabel('Theta')
ax3.set_ylabel('Gamma')
ax3.set_zlabel('DipoleValue')
ax3.set_zlim(mini,maxi)
for j in others:
cutfirst = cut_phi[:,:,:,i,j]
cutx = cutfirst[:,:,0]
cuty = cutfirst[:,:,1]
cutz = cutfirst[:,:,2]
col = colors[j]
ax1.plot_wireframe(X,Y,cutx, alpha=.5, color=col)
ax2.plot_wireframe(X,Y,cuty, alpha=.5, color=col)
ax3.plot_wireframe(X,Y,cutz, alpha=.5, color=col)
fig.tight_layout()
if save:
name_fig = 'dipolematrix_phi{:02d}.png'.format(phi_cut)
fig.savefig(name_fig)
savefigures = False
if savefigures:
for phi_cut in range(phiL):
two_d_see_it(phi_cut,True)
In [ ]:
def two_d_see_it_Only_One(phi_cut,state1,xyz):
cut_phi = dipo[phi_cut,:,:,:,:,:]
X,Y = np.meshgrid(thes,gams)
maxi=3
mini=-3
fig = plt.figure(figsize=(15,10))
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'mediumpurple']
all_states_indexes = np.arange(nstates)
i=state1
others = np.delete(all_states_indexes, [i], axis=0)
ax1 = fig.add_subplot(111, projection='3d')
stringXYZ = ['X','Y','Z']
ax1.set_title('State {} along {}'.format(i,stringXYZ[xyz]))
ax1.set_xlabel('Theta')
ax1.set_ylabel('Gamma')
ax1.set_zlabel('DipoleValue')
ax1.set_zlim(mini,maxi)
for j in others:
cutfirst = cut_phi[:,:,:,i,j]
cutx = cutfirst[:,:,xyz]
col = colors[j]
laby = r'S_{}'.format(j)
#ax1.plot_wireframe(X,Y,cutx, alpha=.5, color=col, label=laby)
ax1.scatter(X,Y,cutx, alpha=.5, color=col, label=laby)
fig.tight_layout()
ax1.legend()
phi_cut = 6
state1 = 6
xyz = 1
two_d_see_it_Only_One(phi_cut,state1,xyz)
In [ ]:
from scipy.interpolate import RegularGridInterpolator
from numpy import linspace, zeros, array
x = linspace(1,4,11)
y = linspace(4,7,22)
z = linspace(7,9,33)
V = zeros((11,22,33))
for i in range(11):
for j in range(22):
for k in range(33):
V[i,j,k] = 100*x[i] + 10*y[j] + z[k]
fn = RegularGridInterpolator((x,y,z), V)
pts = array([[2,6,8],[3,5,7]])
print(fn(pts))
In [ ]:
nstates = 8
# h is for "here" i.e. in this point of the code (not extrapolated)
phi_h = phiV
gam_h = gamV
the_h = np.flip(theV,0)
pp_h,gg_h,tt_h = qp.fromFloatsToLabels(phi_h, gam_h, the_h)
ppL_h = len(pp_h)
ggL_h = len(gg_h)
ttL_h = len(tt_h)
int_pot_linear = np.empty((ppL_h,ggL_h,ttL_h,nstates))
int_pot_nearest = np.empty((ppL_h,ggL_h,ttL_h,nstates))
dointerpolation = False
if dointerpolation:
for s in range(nstates):
print('Doing state {}.'.format(s))
aa = potential[:,:,:,s]
fnn1 = RegularGridInterpolator((phi_h,gam_h,the_h), aa, method='linear')
fnn2 = RegularGridInterpolator((phi_h,gam_h,the_h), aa, method='nearest')
for pi, p in qp.log_progress(enumerate(phi_h),every=1,size=(len(phi_h))):
for gi, g in enumerate(gam_h):
for ti,t in enumerate(the_h):
int_pot_linear[pi,gi,ti,s] = fnn1(array([p,g,t]))
int_pot_nearest[pi,gi,ti,s] = fnn2(array([p,g,t]))
In [ ]:
def fun_pot2D_inter(kind,state,dim):
pot2D(int_pot_linear, state,kind, dim)
#pot2D(int_pot_nearest, state,kind, dim)
fun_pot2D_inter('Phi',1,0)
In [ ]:
# phiV,gamV,theV
In [ ]:
# # gestionate Phi
# refi_phi = 1
# expandi_phi = 0
# phis_n = phiV
# phis_n = expARR(phis_n,expandi_phi)
# for _ in range(refi_phi):
# phis_n = doubleAxespoins(phis_n)
# # gestionate Gam
# refi_gam = 1
# expandi_gam = 0
# gams_n = gamV
# gams_n = expARR(gams_n,expandi_gam)
# for _ in range(refi_gam):
# gams_n = doubleAxespoins(gams_n)
# # gestionate The
# refi_the = 1
# expandi_the = 0
# thes_n = theV
# thes_n = expARR(thes_n,expandi_the)
# for _ in range(refi_the):
# thes_n = doubleAxespoins(thes_n)
# # equilibrium indices and equilibrium points
# p0 = phis_n[(gsm_phi_ind+expandi_phi)*2**refi_phi]
# g0 = gams_n[(gsm_gam_ind+expandi_gam)*2**refi_gam]
# t0 = thes_n[(gsm_the_ind+expandi_the)*2**refi_the]
# print(p0,g0,t0)
# pp,gg,tt = qp.fromFloatsToLabels(phis_n, gams_n, thes_n)
# ppL = len(pp)
# ggL = len(gg)
# ttL = len(tt)
In [ ]:
# phis_n,gams_n,thes_n
In [ ]:
# print(' | eq. #p. dx dL')
# strngOut = 'Phi: {:8.4f} {:4} {:10.3e} {:10.3e} \nGamma: {:8.4f} {:4} {:10.3e} {:10.3e} \nTheta: {:8.4f} {:4} {:10.3e} {:10.3e}'
# print(strngOut.format(p0,len(phis_n),phis_n[1]-phis_n[0],phis_n[-1]-phis_n[0],g0,len(gams_n), gams_n[1]-gams_n[0],gams_n[-1]-gams_n[0],t0,len(thes_n),thes_n[1]-thes_n[0],thes_n[-1]-thes_n[0]))
In [ ]:
# doublALL = np.zeros((ppL,ggL,ttL))
# doublALL2 = np.zeros((ppL,ggL,ttL,8))
# kineti = np.zeros((ppL,ggL,ttL,9,3))
# smo = np.zeros((ppL,ggL,ttL,3,24))
# geom = np.zeros((ppL,ggL,ttL,15,3))
# p_s = 100
# g_s = 100
# t_s = 100
# # Toggle this to calculate harmonic potentials
# saveNew = False
# if saveNew:
# for pi, p in qp.log_progress(enumerate(phis_n),every=1,size=(len(phis_n))):
# for gi, g in enumerate(gams_n):
# for ti,t in enumerate(thes_n):
# doublALL[pi,gi,ti] = p_s*(p-p0)**2 + g_s*(g-g0)**2 + t_s*(t-t0)**2
# kineti[pi,gi,ti] = qp.calc_g_G(p,g,t)
# smo[pi,gi,ti] = qp.calc_s_mat(p,g,t)
# geom = qp.generateNorbGeometry(p,g,t,True) # True is the "vector" mode that does not save any file
# doublALL2[:,:,:,0] = doublALL
# newdict = {}
# newdict['smoCube'] = smo
# newdict['kinCube'] = kineti
# newdict['potCube'] = doublALL2
# newdict['dipCUBE'] = np.zeros((ppL,ggL,ttL,3,8,8))
# newdict['geoCUBE'] = geom
# newdict['phis'] = pp
# newdict['gams'] = gg
# newdict['thes'] = tt
# np.save('doubledoublefinerArmonic.npy', newdict)
# qp.printDict(newdict)
In [ ]:
def unpack_things(vector):
return(vector[0],vector[-1],vector.size,vector[1]-vector[0])
the_0, the_last, theL, the_dt = unpack_things(theV)
gam_0, gam_last, gamL, gam_dt = unpack_things(gamV)
phi_0, phi_last, phiL, phi_dt = unpack_things(phiV)
In [ ]:
#%matplotlib inline
from ipywidgets import interactive,interact, HBox, Layout,VBox
def drawing_abs(a,b,c,reP,length,function):
re = np.pi+ reP
print(re)
# if kinds == 'Phi':
# dim_0, dim_last, dimL, dim_dt = unpack_things(phiV)
# elif kinds == 'Gam':
# dim_0, dim_last, dimL, dim_dt = unpack_things(gamV)
# elif kinds == 'The':
# dim_0, dim_last, dimL, dim_dt = unpack_things(theV)
# elif kinds == 'mon':
# dim_0, dim_last, dimL, dim_dt = -10,10,160,0.5
# print(dim_0, dim_last, dimL, dim_dt)
# x = np.linspace(dim_0,dim_last,dimL)
x = np.linspace(5,10,int(length))
#y = a * (np.arctan(b*(x + c)) + np.arctan(b*(x - c)))
if function == 'symmetric':
y = a * ((-np.arctan(b*(x + c)) + np.arctan(b*(x - c)))+re)
elif function == 'oneside_posi':
y = a * ((np.arctan(b*(x - c)))+re)
elif function == 'oneside_nega':
y = a * (-np.arctan(b*(x + c))+re)
y = np.where(y<0,0,y)
fig, ax = plt.subplots(1,1,figsize=(15,6))
ax.plot(x,y)
ax.set_title('Absorbing Potential')
ax.set_xlabel('q')
ax.set_ylabel('Hartree')
print(y)
widget = interactive(drawing_abs, a = widgets.FloatText(2), b = widgets.FloatText(90), c = widgets.FloatText(10), reP = widgets.FloatText(0.31), length = widgets.FloatText(160), function = ['symmetric','oneside_posi','oneside_nega']);
display(VBox([HBox(widget.children[:-1], layout = Layout(flex_flow='row wrap')), widget.children[-1]]))
In [ ]:
minimo = np.amin(dataDict['potCube'])
potZero = dataDict['potCube'] - minimo
np.amax(potZero[:,:,:,0])
In [14]:
def do3dplot_abs(xs,ys,zss):
'with mesh function'
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X,Y = np.meshgrid(ys,xs)
ax.view_init(elev=1., azim=0.)
ax.set_xlabel('theta')
ax.set_ylabel('gamma')
ax.plot_wireframe(X, Y, zss)
phiLen,gamLen,theLen,nstates = 55, 56, 160, 8
tan1_a, tan1_b, tan1_c, tan1_reP, tan2_len = 1, 0.05, 9.9, 0, 56
tan2_a, tan2_b, tan2_c, tan2_reP, tan2_len = 1, 80, 9.9, -0.01, 160
re_1 = np.pi + tan1_reP
re_2 = np.pi + tan2_reP
#tan1_a,tan1_b,tan1_c,tan1_len = 0.1, 0.01, 9.4, 56
#tan2_a,tan2_b,tan2_c,tan2_len = 2, 10, 8, 160
tan1_X = np.linspace(5,10,tan1_len)
tan2_X = np.linspace(5,10,tan2_len)
y1_neg = tan1_a * ((-np.arctan(tan1_b*(tan1_X + tan1_c)) + np.arctan(tan1_b*(tan1_X - tan1_c)))+re_1)
y1 = np.where(y1_neg<0,0,y1_neg)
y2_neg = tan2_a * ((-np.arctan(tan2_b*(tan2_X + tan2_c)) + np.arctan(tan2_b*(tan2_X - tan2_c)))+re_2)
y2 = np.where(y2_neg<0,0,y2_neg)
yy1,yy2 = np.meshgrid(y2,np.flip(y1,0))
finalGridEV = yy1*yy2
finalGrid = qp.fromEvtoHart(finalGridEV)
print(tan1_X.shape,tan2_X.shape,finalGrid.shape)
print(finalGrid[30,129])
do3dplot_abs(tan1_X,tan2_X,finalGrid)
In [10]:
import os
folder_abs = '/home/alessio/Desktop/Noise_Or_Not/IMPORTANTS/VISUALIZE_ABSORBING'
file_namez = os.path.join(folder_abs,'Abs_THETA_EXAGERATED.h5')
absorbing = np.zeros((phiLen,gamLen,theLen,nstates))
for s in range(nstates):
for p in range(phiLen):
absorbing[p,:,:,s] = finalGrid
diction = {'absorb': absorbing}
save_this_thing = True
if save_this_thing:
qp.writeH5fileDict(file_namez, diction)
In [28]:
dataDict['geoCUBE'].shape
Out[28]:
In [29]:
Mc = qp.massOf('C')
Mh = qp.massOf('H')
def massCenter(geom):
masses = np.array([Mc,Mc,Mc,Mh,Mh,Mh,Mh,Mc,Mc,Mc,Mc,Mh,Mh,Mh,Mh])
total_mass = sum(masses)
zcoor = geom[:,2]
z_cm = sum(masses * zcoor)
return(z_cm/total_mass)
equi = dataDict['geoCUBE'][gsm_phi_ind, gsm_gam_ind, gsm_the_ind]
In [30]:
centerofMass = np.empty((phiL,gamL,theL))
for pi, p in qp.log_progress(enumerate(phiV),every=1,size=(phiL)):
for gi, g in enumerate(gamV):
for ti,t in enumerate(theV):
centerofMass[pi,gi,ti] = massCenter(dataDict['geoCUBE'][pi,gi,ti])
centerofMass.shape
Out[30]:
In [31]:
%matplotlib notebook
do3dplot(gamV,theV,centerofMass[gsm_phi_ind,:,:])
do3dplot(phiV,theV,centerofMass[:,gsm_gam_ind,:])
do3dplot(phiV,gamV,centerofMass[:,:,gsm_the_ind])
In [ ]: