In [1]:
import quantumpropagator as qp
import matplotlib.pyplot as plt
import matplotlib
from ipywidgets import interact,fixed #, interactive, fixed, interact_manual
import ipywidgets as widgets
from tqdm import tqdm_notebook as tqdm
import os
import re
import sys
import glob
import pandas as pd
import numpy as np
%matplotlib ipympl
plt.rcParams.update({'font.size': 12})
In [2]:
%matplotlib ipympl
def read_dipole_file(dipole_file_path):
nstates = 8
permanents = []
transitions = []
for lab1 in ['x','y','z']:
for lab2 in range(nstates):
permanents.append('perm_{}_{}'.format(lab1,lab2))
for lab3 in range(lab2+1,nstates):
transitions.append('trans_{}_{}_{}'.format(lab1,lab2,lab3))
all_labels = ['fs_2','dipx','dipy','dipz'] + permanents + transitions
df_dipo_pump = pd.read_csv(dipole_file_path, delim_whitespace=True, names=all_labels)
return(df_dipo_pump)
def copy_paste_output(pump_df,probe_df, out_name):
'''
This function takes a pump and a probe run and returns a new dataframe with pump_probe pasted together
'''
fs_cut = probe_df['fs'].iloc[0]
# the index on the pump sheet is where the first fs value of the probe is.
index_cut = pump_df.index[pump_df['fs'] == fs_cut][0]
pump_cut = pump_df.iloc[:index_cut]
final_paper = pd.concat((pump_cut,probe_df))
final_paper.to_csv(out_name)
return(final_paper)
def copy_paste_dipos(pump_df_d,probe_df_d, out_name):
'''
This function takes a pump and a probe run and returns a new dataframe with pump_probe pasted together
'''
fs_cut = probe_df_d['fs_2'].iloc[0]
# the index on the pump sheet is where the first fs value of the probe is.
index_cut = pump_df_d.index[pump_df_d['fs_2'] == fs_cut][0]
pump_cut = pump_df_d.iloc[:index_cut]
final_paper = pd.concat((pump_cut,probe_df_d))
final_paper.to_csv(out_name)
return(final_paper)
def pythonfft(signal_time,signal,pad_length):
'''
fast fourier transform
signal_time :: np.array(double) <- the time array in femtoseconds
signal :: np.array(double) <- the signal in hartree
pad_length :: the length of the padding to increase resolution
returns :: frequencies in eV and the positive side of the signal fft
'''
signal_fft = np.fft.fft(np.pad(signal,(0,pad_length)))
time_au = qp.fromFsToAu(signal_time)
dt = time_au[1]- time_au[0]
frequancies_from_routine = np.fft.fftfreq(time_au.size + pad_length)
time_au_max = np.amax(time_au)
frequencies = qp.fromHartoEv(frequancies_from_routine*2*np.pi/dt)
indexes = np.where(frequancies_from_routine>=0)
return(frequencies[indexes],signal_fft[indexes])
In [3]:
root = '/home/alessio/w-August-Run/'
pump = 'b-UV-0.22_0000'
probe = 'y-probe-UV_0000'
fol_pump = os.path.join(root, pump)
fol_probe = os.path.join(root, probe)
output_norm_pump = os.path.join(fol_pump, 'output')
output_dipo_pump = os.path.join(fol_pump, 'Output_Dipole')
output_norm_probe = os.path.join(fol_probe, 'output')
output_dipo_probe = os.path.join(fol_probe, 'Output_Dipole')
df_norm2_pump = pd.read_csv(output_norm_pump, delim_whitespace=True, index_col=0, names=['counter', 'steps', 'fs','Norm deviation','Kinetic','Potential','Total','Total Deviation','Xpulse','Ypulse','Zpulse','AbZino'])
df_norm2_probe = pd.read_csv(output_norm_probe, delim_whitespace=True, index_col=0, names=['counter', 'steps', 'fs','Norm deviation','Kinetic','Potential','Total','Total Deviation','Xpulse','Ypulse','Zpulse','AbZino'])
df_dipo_pump = read_dipole_file(output_dipo_pump)
df_dipo_probe = read_dipole_file(output_dipo_probe)
ppname = os.path.join('/home/alessio/w-August-Run/','{}_{}.csv'.format(pump,probe))
ppname_dip = os.path.join('/home/alessio/w-August-Run/','{}_{}_dipole.csv'.format(pump,probe))
df1 = copy_paste_output(df_norm2_pump, df_norm2_probe, ppname)
dip1 = copy_paste_dipos(df_dipo_pump, df_dipo_probe, ppname_dip)
In [4]:
# cut the shortest, please
(df1.shape, dip1.shape)
num1, _ = df1.shape
num2, _ = dip1.shape
cut_here = np.amin([num1,num2])
df = df1[:cut_here]
dip = dip1[:cut_here]
In [5]:
x = df['fs']
pulsex = df['Xpulse']
pulsey = df['Ypulse']
pulsez = df['Zpulse']
together = np.stack((pulsex,pulsey,pulsez), axis=1)
y = np.sum(together,axis=1)
fig, ax0 = plt.subplots(1,1,figsize = (10,6))
ax0.plot(x,y)
ax0.set_xlabel('fs')
fig.tight_layout()
In [6]:
pad_length = 10000
freq, signal_fft = pythonfft(np.array(x),np.array(y),pad_length)
fig, ax0 = plt.subplots(1,1,figsize = (10,6))
#ax0.plot(freq, signal_fft)
ax0.plot(freq, np.real(signal_fft))
ax0.plot(freq, np.imag(signal_fft))
ax0.set_xlim(0,10)
ax0.set_xlabel('eV')
fig.tight_layout()
In [7]:
if False:
plt.close('all')
pad_length = 10000
xx = dip['fs_2']
x = df['fs']
#for ylabel in ['dipx', 'dipy', 'dipz','trans_x_0_1']:
for ylabel in ['dipx', 'dipy', 'dipz', 'trans_x_0_1', 'trans_z_0_1']: # threshold component also.
pulsex = df['Xpulse']
pulsey = df['Ypulse']
pulsez = df['Zpulse']
together = np.stack((pulsex,pulsey,pulsez), axis=1)
y = np.sum(together, axis=1)
yy = dip[ylabel]
freq_dip, signal_fft_dip = pythonfft(np.array(xx),np.array(yy),pad_length)
fig, [[ax0, ax1], [ax2, ax3], [ax4, ax5]] = plt.subplots(3,2,figsize = (10,6))
#ax0.plot(freq_dip, signal_fft_dip)
ax0.plot(xx,yy)
ax0.set_xlim(0,100)
ax0.set_xlabel('fs')
ax0.set_title('{} - signal'.format(ylabel))
ax2.plot(freq, np.real(signal_fft_dip))
ax2.plot(freq, np.imag(signal_fft_dip))
ax2.set_xlim(0,10)
ax2.set_xlabel('eV')
ax2.set_title('{} - FT'.format(ylabel))
ax1.plot(x,y)
ax1.set_xlabel('fs')
ax1.set_title('{} - Pulse'.format(ylabel))
freq, signal_fft = pythonfft(np.array(x),np.array(y),pad_length)
ax3.plot(freq, np.real(signal_fft))
ax3.plot(freq, np.imag(signal_fft))
ax3.set_xlim(0,10)
ax3.set_xlabel('eV')
ax3.set_title('{} - Pulse FFT'.format(ylabel))
S_W = -2 * np.imag(np.conj(signal_fft) * signal_fft_dip)
ax4.plot(freq_dip, S_W)
ax4.set_xlim(0,10)
ax4.set_xlabel('eV')
ax4.set_title('{} - Transient Absorption Spectra'.format(ylabel))
#fig.suptitle(ylabel)
fig.canvas.layout.height = '900px'
fig.tight_layout()
In [8]:
#plt.close('all')
if False:
outputFolder = '/home/alessio/w-August-Run/OUTPUTS'
name_output = os.path.join(outputFolder,'{}_{}.csv'.format(pump,probe))
print(name_output)
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'mediumpurple']
pad_length = 10000
nstates = 8
permanents = []
transitions = []
fran_transition_time_sheet = pd.DataFrame()
fran_transition_freq_sheet = pd.DataFrame()
for lab1 in ['x','y','z']:
for lab2 in range(nstates):
permanents.append('perm_{}_{}'.format(lab1,lab2))
for lab3 in range(lab2+1,nstates):
transitions.append('trans_{}_{}_{}'.format(lab1,lab2,lab3))
pulsex = df['Xpulse']
pulsey = df['Ypulse']
pulsez = df['Zpulse']
fran_transition_time_sheet['Xpulse'] = pulsex
fran_transition_time_sheet['Ypulse'] = pulsey
fran_transition_time_sheet['Zpulse'] = pulsez
for cart in ['x','z']:
#for cart in ['x']:
timefs = dip['fs_2']
time = qp.fromFsToAu(timefs)
this_full_dipole = 'dip{}'.format(cart)
tot_signal = dip[this_full_dipole]
full_list_label = permanents + transitions
full_list_this_cart = [x for x in full_list_label if cart in x]
threshold = np.linalg.norm(tot_signal)*0.05 # 5% of the norm
big_components = [ x for x in full_list_this_cart if np.linalg.norm(dip[x]) > threshold or x in ['trans_x_0_1', 'trans_y_0_1', 'trans_z_0_1'] ] # the one which contributes more than 5% to the norm
# make figure
#fig, [ax0, ax1, ax2, ax3] = plt.subplots(4,1,figsize = (10,6))
fig, [[ax0, ax1], [ax2, ax3], [ax4, ax5], [ax6, ax7]] = plt.subplots(4,2,figsize = (10,6))
# PULSE PANEL
pulse_this = df['{}pulse'.format(cart.upper())]
ax0.plot(timefs,pulse_this)
freq_pulse, pulse_fft = pythonfft(np.array(time),pulse_this,pad_length)
pulse_fft_RE = np.real(pulse_fft)
pulse_fft_IM = np.imag(pulse_fft)
freq_pulse_eV = qp.fromHartoEv(freq_pulse)
# send to Francoise Spreadsheet
freq_label = 'Frequencies eV'
if freq_label not in fran_transition_freq_sheet.columns:
fran_transition_freq_sheet[freq_label] = freq_pulse_eV
fran_transition_freq_sheet['{}_pulse_ft_real'.format(cart)] = pulse_fft_RE
fran_transition_freq_sheet['{}_pulse_ft_imag'.format(cart)] = pulse_fft_IM
ax1.plot(freq_pulse_eV, pulse_fft_RE)
ax1.plot(freq_pulse_eV, pulse_fft_IM)
ax1.set_xlim(0,10)
# DIPOLE PANEL
ax2.plot(timefs,tot_signal)
freq_dipole_tot, fft_dipole_tot = pythonfft(np.array(time),tot_signal,pad_length)
freq_dipole_eV = qp.fromHartoEv(freq_dipole_tot)
dipole_total_fft_RE = np.real(fft_dipole_tot)
dipole_total_fft_IM = np.imag(fft_dipole_tot)
# send to Francoise Spreadsheet
fran_transition_freq_sheet['{}_dipole_ft_real'.format(cart)] = pulse_fft_RE
fran_transition_freq_sheet['{}_dipole_ft_imag'.format(cart)] = pulse_fft_IM
ax3.plot(freq_dipole_eV, dipole_total_fft_RE)
ax3.plot(freq_dipole_eV, dipole_total_fft_IM)
ax3.set_xlim(0,10)
# first spectra
S_W_tot = -2 * np.imag(np.conj(pulse_fft) * fft_dipole_tot)
fran_transition_freq_sheet['TA {}'.format(cart)] = S_W_tot
ax6.plot(freq_dipole_eV, S_W_tot)
ax6.set_xlim(0,10)
# COMPONENT PANEL
print(big_components)
for iii, component in enumerate(big_components):
comp_dip = dip[component]
ax4.plot(timefs,comp_dip, label = component, color=colors[iii])
freq_component_tot, fft_component_tot = pythonfft(np.array(time), comp_dip, pad_length)
freq_component_eV = qp.fromHartoEv(freq_component_tot)
component_fft_RE = np.real(fft_component_tot)
component_fft_IM = np.imag(fft_component_tot)
# send to Francoise Spreadsheet
fran_transition_freq_sheet['{}_{}_component_ft_real'.format(cart,component)] = component_fft_RE
fran_transition_freq_sheet['{}_{}_component_ft_imag'.format(cart,component)] = component_fft_IM
ax5.plot(freq_component_eV, component_fft_RE, color=colors[iii])
ax5.plot(freq_component_eV, component_fft_IM, color=colors[iii], ls=':')
ax5.set_xlim(0,10)
# Spectra
S_W_tot_component = -2 * np.imag(np.conj(pulse_fft) * fft_component_tot)
fran_transition_freq_sheet['TA {}_{}'.format(cart,component)] = S_W_tot_component
ax7.plot(freq_dipole_eV, S_W_tot_component, color=colors[iii])
ax7.set_xlim(0,10)
ax4.legend()
fig.canvas.layout.height = '900px'
fig.tight_layout()
fran_transition_freq_sheet.to_csv(name_output)
In [9]:
def new_pulse(t,Ed,omega,sigma,phi,t0):
num = (t-t0)**2
den = 2*(sigma**2)
if (den == 0):
result = 0.0
else:
num2 = np.sin(omega*(t-t0) + phi) * (t-t0)
den2 = omega * sigma**2
result = Ed * np.exp(-num/den) * (np.cos(omega*(t-t0) + phi) - num2/den2 )
return result
In [10]:
plt.close('all')
pad_length = 10000
length_fs = 100
resolution = 2000
times = np.linspace(0, qp.fromFsToAu(length_fs), resolution)
E = 0.06
omega = 0.22#/(2*np.pi)
sigma = qp.fromFsToAu(1.7)
phi = 0
t0 = 2500
y = qp.pulse(times, E, omega, sigma, phi, t0)
y2 = new_pulse(times, E, omega, sigma, phi, t0)
x = qp.fromAuToFs(times)
freq, signal_fft = pythonfft(np.array(x),np.array(y),pad_length)
freq, signal_fft2 = pythonfft(np.array(x),np.array(y2),pad_length)
colors = ['r','g']
fig, [ ax0, ax1 ] = plt.subplots(2,1,figsize = (10,6))
ax0.plot(x,y, color = colors[0], label='Old Pulse')
ax0.plot(x,y2, color = colors[1], label='New Pulse')
ax0.set_xlabel('fs')
## IMAG
ax1.plot(freq, np.real(signal_fft), color = colors[0])
ax1.plot(freq, np.imag(signal_fft), color = colors[0], ls=':')
ax1.plot(freq, np.real(signal_fft2), color = colors[1])
ax1.plot(freq, np.imag(signal_fft2), color = colors[1], ls=':')
## ABSOLUTE
# ax1.plot(freq, np.abs(signal_fft), color = colors[0])
# ax1.plot(freq, np.abs(signal_fft2), color = colors[1])
major_ticks = np.arange(0.5, 10, 0.5)
minor_ticks = np.arange(0, 10, 0.1)
ax1.set_xticks(major_ticks)
ax1.set_xticks(minor_ticks, minor=True)
ax1.set_xlim(0,10)
ax1.set_xlabel('eV')
ax0.legend()
fig.tight_layout()
In [11]:
from scipy import signal
def process_pump_probe_folder(root, pump, probe, cartesian, do_graph = True, select=None, limit_time = None):
'''
this is a function which does everything for pump-probe-cross graphs
'''
# OUTPUTS
select = select or ['trans_z_0_1']
limit_time = limit_time or 9999999
outputFolder = '/home/alessio/w-August-Run/OUTPUTS'
name_output = os.path.join(outputFolder,'{}_{}.csv'.format(pump,probe))
fol_pump = os.path.join(root, pump)
fol_probe = os.path.join(root, probe)
output_norm_pump = os.path.join(fol_pump, 'output')
output_dipo_pump = os.path.join(fol_pump, 'Output_Dipole')
output_norm_probe = os.path.join(fol_probe, 'output')
output_dipo_probe = os.path.join(fol_probe, 'Output_Dipole')
df_norm2_pump = pd.read_csv(output_norm_pump, delim_whitespace=True, index_col=0, names=['counter', 'steps', 'fs','Norm deviation','Kinetic','Potential','Total','Total Deviation','Xpulse','Ypulse','Zpulse','AbZino'])
df_norm2_probe = pd.read_csv(output_norm_probe, delim_whitespace=True, index_col=0, names=['counter', 'steps', 'fs','Norm deviation','Kinetic','Potential','Total','Total Deviation','Xpulse','Ypulse','Zpulse','AbZino'])
df_dipo_pump = read_dipole_file(output_dipo_pump)
df_dipo_probe = read_dipole_file(output_dipo_probe)
ppname = os.path.join(outputFolder,'{}_{}.csv'.format(pump,probe))
ppname_dip = os.path.join(outputFolder,'{}_{}_dipole.csv'.format(pump,probe))
df1 = copy_paste_output(df_norm2_pump, df_norm2_probe, ppname)
dip1 = copy_paste_dipos(df_dipo_pump, df_dipo_probe, ppname_dip)
# cut the shortest, please
(df1.shape, dip1.shape)
num1, _ = df1.shape
num2, _ = dip1.shape
print('out {}: {} --- dip {}: {}'.format(num1,df1['fs'].iloc[-1], num2, dip1['fs_2'].iloc[-1]))
cut_here = np.amin([num1,num2,limit_time]) # which is the shortest?
if cut_here == limit_time: # it means that we are cutting
print('limit {}: {}'.format(limit_time,df1['fs'].iloc[cut_here]))
df = df1[:cut_here]
dip = dip1[:cut_here]
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'mediumpurple']
pad_length_until = 10000
pad_length = pad_length_until - cut_here
nstates = 8
permanents = []
transitions = []
fran_transition_time_sheet = pd.DataFrame()
fran_transition_freq_sheet = pd.DataFrame()
for lab1 in ['x','y','z']:
for lab2 in range(nstates):
permanents.append('perm_{}_{}'.format(lab1,lab2))
for lab3 in range(lab2+1,nstates):
transitions.append('trans_{}_{}_{}'.format(lab1,lab2,lab3))
pulsex = df['Xpulse']
pulsey = df['Ypulse']
pulsez = df['Zpulse']
fran_transition_time_sheet['Xpulse'] = pulsex
fran_transition_time_sheet['Ypulse'] = pulsey
fran_transition_time_sheet['Zpulse'] = pulsez
#for cart in ['x','y','z']:
for cart in [cartesian]:
timefs = dip['fs_2']
time = qp.fromFsToAu(timefs)
this_full_dipole = 'dip{}'.format(cart)
tot_signal = dip[this_full_dipole]
full_list_label = permanents + transitions
full_list_this_cart = [x for x in full_list_label if cart in x]
threshold = np.linalg.norm(tot_signal)*0.05 # 5% of the norm
big_components = [ x for x in full_list_this_cart if np.linalg.norm(dip[x]) > threshold or x in ['trans_x_0_1', 'trans_y_0_1', 'trans_z_0_1'] ]
# Create Arrays
pulse_this = df['{}pulse'.format(cart.upper())]
freq_pulse, pulse_fft = pythonfft(np.array(time),pulse_this,pad_length)
pulse_fft_RE = np.real(pulse_fft)
pulse_fft_IM = np.imag(pulse_fft)
freq_pulse_eV = qp.fromHartoEv(freq_pulse)
time_normal = np.array(time)
#window = gau_window(timefs, 10, 75)
window = signal.gaussian(timefs.shape[0], std=450)
tot_signal_windowed = tot_signal * window
freq_dipole_tot, fft_dipole_tot = pythonfft(np.array(time), tot_signal_windowed, pad_length)
freq_dipole_eV = qp.fromHartoEv(freq_dipole_tot)
dipole_total_fft_RE = np.real(fft_dipole_tot)
dipole_total_fft_IM = np.imag(fft_dipole_tot)
# first spectra
S_W_tot = -2 * np.imag(np.conj(pulse_fft) * fft_dipole_tot)
if do_graph:
# CREATE GRAPHS
fig, [[ax0, ax1], [ax2, ax3], [ax4, ax5], [ax6, ax7]] = plt.subplots(4,2,figsize = (10,6))
ax0.plot(timefs,pulse_this)
ax0.set_title('Pulse')
ax1.plot(freq_pulse_eV, pulse_fft_RE)
ax1.plot(freq_pulse_eV, pulse_fft_IM)
ax1.set_xlim(0,10)
ax1.set_title('Pulse FT')
ax2.plot(timefs,tot_signal_windowed)
ax2.plot(timefs,window)
ax2.set_title('Signal and window')
ax3.plot(freq_dipole_eV, dipole_total_fft_RE)
ax3.plot(freq_dipole_eV, dipole_total_fft_IM)
ax3.set_xlim(0,10)
ax3.set_title('Signal windower FT')
ax6.plot(freq_dipole_eV, S_W_tot)
ax6.set_xlim(0,10)
ax6.set_title('Transient Absorbtion spectra total')
# send to Francoise Spreadsheet
freq_label = 'Frequencies eV'
if freq_label not in fran_transition_freq_sheet.columns:
fran_transition_freq_sheet[freq_label] = freq_pulse_eV
fran_transition_freq_sheet['{}_pulse_ft_real'.format(cart)] = pulse_fft_RE
fran_transition_freq_sheet['{}_pulse_ft_imag'.format(cart)] = pulse_fft_IM
fran_transition_freq_sheet['{}_dipole_ft_real'.format(cart)] = dipole_total_fft_RE
fran_transition_freq_sheet['{}_dipole_ft_imag'.format(cart)] = dipole_total_fft_IM
fran_transition_freq_sheet['TA {}'.format(cart)] = S_W_tot
# COMPONENT PANEL
# for iii, component in enumerate(big_components):
for iii, component in enumerate(select):
comp_dip = dip[component]
freq_component_tot, fft_component_tot = pythonfft(np.array(time), comp_dip, pad_length)
freq_component_eV = qp.fromHartoEv(freq_component_tot)
component_fft_RE = np.real(fft_component_tot)
component_fft_IM = np.imag(fft_component_tot)
# Spectra
S_W_tot_component = -2 * np.imag(np.conj(pulse_fft) * fft_component_tot)
if do_graph:
ax4.plot(timefs,comp_dip, color=colors[iii])
ax4.set_title('{} component on the dipole'.format(component))
ax5.plot(freq_component_eV, component_fft_RE, color=colors[iii])
ax5.plot(freq_component_eV, component_fft_IM, color=colors[iii], ls=':')
ax5.set_xlim(0,10)
ax5.set_title('{} component FT'.format(component))
ax7.plot(freq_dipole_eV, S_W_tot_component, color=colors[iii])
ax7.set_xlim(0,10)
ax7.set_title('{} Transient Absoption'.format(component))
plt.suptitle('{}'.format(probe))
# send to Francoise Spreadsheet
fran_transition_freq_sheet['{}_{}_component_ft_real'.format(cart,component)] = component_fft_RE
fran_transition_freq_sheet['{}_{}_component_ft_imag'.format(cart,component)] = component_fft_IM
fran_transition_freq_sheet['TA {}_{}'.format(cart,component)] = S_W_tot_component
if do_graph:
fig.canvas.layout.height = '900px'
fig.tight_layout()
fran_transition_freq_sheet.to_csv(name_output)
return(freq_dipole_eV, S_W_tot_component, S_W_tot)
In [12]:
def gau_window(t,sigma,t0):
'''
It returns the value of the gaussian envelope for the pulse at time t
now it works with arrays, too
'''
num = (t-t0)**2
den = 2*(sigma**2)
if (den == 0):
if type(t) == float:
result = 0.0
else:
result = np.zeros_like(t)
else:
result = np.exp(-num/den)
norm = np.linalg.norm(result)
return result/norm
In [13]:
import re
root = '/home/alessio/w-August-Run/'
pump = 'b-UV-0.22_0000'
folders = sorted(glob.glob('/home/alessio/w-August-Run/a-UV-Pump-Probe*'))
print(folders)
delays_au = [ -float(re.findall("[-+]?[.]?[\d]+(?:,\d\d\d)*[\.]?\d*(?:[eE][-+]?\d+)?", s)[0]) for s in folders]
In [17]:
graph = True
all_of_the_1_0 = []
all_of_the_TOT = []
all_the_freq = []
for this_one in tqdm(folders):
freq_dipole_eV, S_W_tot_component, S_W_tot = process_pump_probe_folder(root, pump, os.path.basename(this_one), 'z', do_graph = graph, select=['trans_z_0_1'])#, limit_time = 2400)
all_of_the_1_0.append(S_W_tot_component)
all_of_the_TOT.append(S_W_tot)
all_the_freq.append(freq_dipole_eV)
In [18]:
# save in files
save_this = False
folder_output = '/STORE/alessio/things/dox/Alessio-Francoise/Results_MD_October2019/Transient_absorption/UVx'
appended_vectors = [all_of_the_TOT, all_of_the_1_0]
vector_labels = ['ALL','0_1_component']
delays_fs = [ '{:.1f}'.format(qp.fromAuToFs(x)) for x in delays_au ]
energies_vector = np.array(all_the_freq[0])
if save_this:
for i, this_label in enumerate(vector_labels):
all_freq = np.array(appended_vectors[i])
df = pd.DataFrame(all_freq.T, index=energies_vector, columns=delays_fs)
df.index.name = 'Freq eV'
output_name = os.path.join(folder_output, 'Transient_abs_spectra_UVx_{}.txt'.format(this_label))
df.to_csv(output_name)
print('File {} written...'.format(output_name))
In [16]:
from matplotlib import cm
cmap = 'hot'#'hsv'
viridis = cm.get_cmap(cmap,12)
color_map = viridis(np.linspace(0,1, len(delays_au)+2))
fig, ax0 = plt.subplots(1,1,figsize = (10,6))
ax1 = ax0.twinx()
how_many = len(delays_au)
for i,(x,y,z) in enumerate(zip(all_the_freq[:how_many], all_of_the_1_0[:how_many], all_of_the_TOT[:how_many])):
ax0.plot(x,y, label=delays_fs[i], ls=':', color=color_map[i])
ax1.plot(x,z, color=color_map[i])
ax0.set_xlim(0.0,0.8)
ax1.set_xlim(0.0,0.8)
ax0.legend()
ax0.set_xlabel('eV')
fig.tight_layout()
In [17]:
together = np.stack(all_of_the_1_0, axis=0)
together2 = np.stack(all_of_the_TOT, axis=0)
pp_runs = zip(folders,delays_fs)
together.shape, freq_dipole_eV.shape, len(list(pp_runs)), together2.shape, delays_fs
Out[17]:
In [18]:
vector_not_cut = together
cut_from = 0
cut_at = 130
vector = vector_not_cut[:,cut_from:cut_at]
y = freq_dipole_eV[cut_from:cut_at]
_, vector_y_dim = vector.shape
fig, [ ax1, ax2 ] = plt.subplots(1,2,figsize = (10,6))
z_min, z_max = -np.abs(vector).max(), np.abs(vector).max()
inter = 'none' # 'lanczos','none','gaussian','bilinear' 'spline16'
inter2 = 'lanczos'
# I put the c variable to set colorbar
c = ax1.imshow(vector.T, interpolation=inter, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
d = ax2.imshow(vector.T, interpolation=inter2, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
# this is easy because I want to put a label for each field I have
where_xtiks = range(len(delays_fs))
label_x = delays_fs
# this is more cumbersome, as I want to put a limited amount of tiks into the y axis.
how_many_y_ticks = 7 # how many I want
where_ytiks = range(0,vector_y_dim,int((vector_y_dim/(how_many_y_ticks-1)))) # this are the indexes at which I want to put yticks
label_y = [ '{:.2f}'.format(x) for x in y[where_ytiks]] # value of y axis at index where_yticks
fig.canvas.layout.height = '900px'
fig.colorbar(c, ax=ax1)
fig.colorbar(d, ax=ax2)
ax1.set_xlabel(r'$\tau$')
ax1.set_ylabel(r'eV')
ax1.set_xticks(where_xtiks)
ax1.set_xticklabels(delays_fs)
ax1.set_yticks(where_ytiks)
ax1.set_yticklabels(label_y)
ax1.set_title('Row data')
ax2.set_xlabel(r'$\tau$')
ax2.set_ylabel(r'eV')
ax2.set_xticks(where_xtiks)
ax2.set_xticklabels(delays_fs)
ax2.set_yticks(where_ytiks)
ax2.set_yticklabels(label_y)
ax2.set_title('interpolated')
fig.tight_layout()
In [19]:
vector_not_cut = together2
cut_from = 10
cut_at = 120
vector = vector_not_cut[:,cut_from:cut_at]
y = freq_dipole_eV[cut_from:cut_at]
_, vector_y_dim = vector.shape
fig, [ ax1, ax2 ] = plt.subplots(1,2,figsize = (2,6))
z_min, z_max = -np.abs(vector).max(), np.abs(vector).max()
# https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/interpolation_methods.html
inter = 'none' # 'lanczos' 'none' 'gaussian' 'bilinear'
inter2 = 'spline16'
# I put the c variable to set colorbar
c = ax1.imshow(vector.T, interpolation=inter, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
d = ax2.imshow(vector.T, interpolation=inter2, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
# this is easy because I want to put a label for each field I have
where_xtiks = range(len(delays_fs))
label_x = delays_fs
# this is more cumbersome, as I want to put a limited amount of tiks into the y axis.
how_many_y_ticks = 7 # how many I want
where_ytiks = range(0,vector_y_dim,int((vector_y_dim/(how_many_y_ticks-1)))) # this are the indexes at which I want to put yticks
label_y = [ '{:.2f}'.format(x) for x in y[where_ytiks]] # value of y axis at index where_yticks
fig.canvas.layout.height = '900px'
fig.colorbar(c, ax=ax1)
fig.colorbar(d, ax=ax2)
ax1.set_xlabel(r'$\tau$')
ax1.set_ylabel(r'eV')
ax1.set_xticks(where_xtiks)
ax1.set_xticklabels(delays_fs)
ax1.set_yticks(where_ytiks)
ax1.set_yticklabels(label_y)
ax1.set_title('Row data')
ax2.set_xlabel(r'$\tau$')
ax2.set_ylabel(r'eV')
ax2.set_xticks(where_xtiks)
ax2.set_xticklabels(delays_fs)
ax2.set_yticks(where_ytiks)
ax2.set_yticklabels(label_y)
ax2.set_title('interpolated')
fig.tight_layout()
In [48]:
# This passes the vector in gnuplot format
if False:
xL, yL = vector.shape
print(vector.shape)
xs = np.arange(xL)
ys = np.arange(yL)
output_file = '/home/alessio/Desktop/vectorGnup'
with open(output_file,'w') as out:
for x in xs:
for y in ys:
out.write('{} {} {}\n'.format(x,y,vector[x,y]))
out.write('\n')
In [49]:
root = '/home/alessio/w-August-Run/'
#root = '/home/alessio/Desktop/Dropbox/Cast_helper/Pump_Probe'
pump = 'u-target-2-3-only_0000'
folders = sorted(glob.glob('/home/alessio/w-August-Run/d-UVz-Pump-Probe*'))
#folders = sorted(glob.glob('/home/alessio/Desktop/Dropbox/Cast_helper/Pump_Probe/d-UVz-Pump-Probe*'))
print(folders)
delays_au = [ -float(re.findall("[-+]?[.]?[\d]+(?:,\d\d\d)*[\.]?\d*(?:[eE][-+]?\d+)?", s)[0]) for s in folders]
In [50]:
graph = False
all_of_the_2_3 = []
all_of_the_TOT = []
all_the_freq = []
for this_one in tqdm(folders):
freq_dipole_eV, S_W_tot_component, S_W_tot = process_pump_probe_folder(root, pump, os.path.basename(this_one), 'y', do_graph = graph, select=['trans_y_2_3'])
all_of_the_2_3.append(S_W_tot_component)
all_of_the_TOT.append(S_W_tot)
all_the_freq.append(freq_dipole_eV)
In [51]:
# save in files
folder_output = '/STORE/alessio/things/dox/Alessio-Francoise/Results_MD_October2019/Transient_absorption/UVz'
delays_fs = [ '{:.1f}'.format(qp.fromAuToFs(x)) for x in delays_au ]
if False:
appended_vectors = [all_of_the_TOT, all_of_the_2_3]
vector_labels = ['ALL','2_3_component']
energies_vector = np.array(all_the_freq[0])
for i, this_label in enumerate(vector_labels):
all_freq = np.array(appended_vectors[i])
df = pd.DataFrame(all_freq.T, index=energies_vector, columns=delays_fs)
df.index.name = 'Freq eV'
output_name = os.path.join(folder_output, 'Transient_abs_spectra_UVz_{}.txt'.format(this_label))
df.to_csv(output_name)
print('File {} written...'.format(output_name))
In [52]:
from matplotlib import cm
cmap = 'hot'#'hsv'
viridis = cm.get_cmap(cmap,12)
color_map = viridis(np.linspace(0,1, len(delays_au)+2))
fig, ax0 = plt.subplots(1,1,figsize = (10,6))
ax1 = ax0.twinx()
how_many = len(delays_au)
for i,(x,y,z) in enumerate(zip(all_the_freq[:how_many], all_of_the_2_3[:how_many], all_of_the_TOT[:how_many])):
ax0.plot(x,y, label=delays_fs[i], ls=':', color=color_map[i])
ax1.plot(x,z, color=color_map[i])
ax0.set_xlim(0,10)
ax1.set_xlim(0,10)
ax0.legend()
ax0.set_xlabel('eV')
fig.tight_layout()
In [53]:
together = np.stack(all_of_the_2_3, axis=0)
together2 = np.stack(all_of_the_TOT, axis=0)
pp_runs = zip(folders,delays_fs)
together.shape, freq_dipole_eV.shape, len(list(pp_runs)), together2.shape, delays_fs
Out[53]:
In [54]:
vector_not_cut = together
cut_from = 0
cut_at = 120
vector = vector_not_cut[:,cut_from:cut_at]
y = freq_dipole_eV[cut_from:cut_at]
_, vector_y_dim = vector.shape
fig, [ ax1, ax2 ] = plt.subplots(1,2,figsize = (10,6))
z_min, z_max = -np.abs(vector).max(), np.abs(vector).max()
inter = 'none' # 'lanczos','none','gaussian','bilinear' 'spline16'
inter2 = 'lanczos'
# I put the c variable to set colorbar
c = ax1.imshow(vector.T, interpolation=inter, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
d = ax2.imshow(vector.T, interpolation=inter2, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
# this is easy because I want to put a label for each field I have
where_xtiks = range(len(delays_fs))
label_x = delays_fs
# this is more cumbersome, as I want to put a limited amount of tiks into the y axis.
how_many_y_ticks = 7 # how many I want
where_ytiks = range(0,vector_y_dim,int((vector_y_dim/(how_many_y_ticks-1)))) # this are the indexes at which I want to put yticks
label_y = [ '{:.2f}'.format(x) for x in y[where_ytiks]] # value of y axis at index where_yticks
fig.canvas.layout.height = '900px'
fig.colorbar(c, ax=ax1)
fig.colorbar(d, ax=ax2)
ax1.set_xlabel(r'$\tau$')
ax1.set_ylabel(r'eV')
ax1.set_xticks(where_xtiks)
ax1.set_xticklabels(delays_fs)
ax1.set_yticks(where_ytiks)
ax1.set_yticklabels(label_y)
ax1.set_title('Row data')
ax2.set_xlabel(r'$\tau$')
ax2.set_ylabel(r'eV')
ax2.set_xticks(where_xtiks)
ax2.set_xticklabels(delays_fs)
ax2.set_yticks(where_ytiks)
ax2.set_yticklabels(label_y)
ax2.set_title('interpolated')
fig.tight_layout()
In [55]:
vector_not_cut = together2
cut_from = 0
cut_at = 120
vector = vector_not_cut[:,cut_from:cut_at]
y = freq_dipole_eV[cut_from:cut_at]
_, vector_y_dim = vector.shape
fig, [ ax1, ax2 ] = plt.subplots(1,2,figsize = (2,6))
z_min, z_max = -np.abs(vector).max(), np.abs(vector).max()
# https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/interpolation_methods.html
inter = 'none' # 'lanczos' 'none' 'gaussian' 'bilinear'
inter2 = 'spline16'
# I put the c variable to set colorbar
c = ax1.imshow(vector.T, interpolation=inter, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
d = ax2.imshow(vector.T, interpolation=inter2, cmap='RdBu', aspect='auto', origin='lower', vmin=z_min, vmax=z_max)
# this is easy because I want to put a label for each field I have
where_xtiks = range(len(delays_fs))
label_x = delays_fs
# this is more cumbersome, as I want to put a limited amount of tiks into the y axis.
how_many_y_ticks = 7 # how many I want
where_ytiks = range(0,vector_y_dim,int((vector_y_dim/(how_many_y_ticks-1)))) # this are the indexes at which I want to put yticks
label_y = [ '{:.2f}'.format(x) for x in y[where_ytiks]] # value of y axis at index where_yticks
fig.canvas.layout.height = '900px'
fig.colorbar(c, ax=ax1)
fig.colorbar(d, ax=ax2)
ax1.set_xlabel(r'$\tau$')
ax1.set_ylabel(r'eV')
ax1.set_xticks(where_xtiks)
ax1.set_xticklabels(delays_fs)
ax1.set_yticks(where_ytiks)
ax1.set_yticklabels(label_y)
ax1.set_title('Row data')
ax2.set_xlabel(r'$\tau$')
ax2.set_ylabel(r'eV')
ax2.set_xticks(where_xtiks)
ax2.set_xticklabels(delays_fs)
ax2.set_yticks(where_ytiks)
ax2.set_yticklabels(label_y)
ax2.set_title('interpolated')
fig.tight_layout()
In [14]:
all_probes_folders = sorted(glob.glob(root + 'o-UV*'))
delays = [68.0,68.2,68.4,68.6,68.8,69.0,69.2,69.4,69.6,69.8,70.0,70.2,70.4,64.0,65.0,66.0,67.0,71.0,72.0,73.0,74.0,75.0]
graph = False
all_of_the_1_0 = []
all_of_the_TOT = []
all_the_freq = []
all_probes_folders_selected = all_probes_folders[13:17] + [ all_probes_folders[0], all_probes_folders[5], all_probes_folders[10] ] + all_probes_folders[17:]
for this_one in tqdm(all_probes_folders_selected):
freq_dipole_eV, S_W_tot_component, S_W_tot = process_pump_probe_folder(root, pump, os.path.basename(this_one), do_graph = graph)
all_of_the_1_0.append(S_W_tot_component)
all_of_the_TOT.append(S_W_tot)
all_the_freq.append(freq_dipole_eV)
In [15]:
from matplotlib import cm
delays = [64.0,65.0,66.0,67.0,68.0,69.0,70.0,71.0,72.0,73.0,74.0,75.0]
cmap = 'hot'#'hsv'
viridis = cm.get_cmap(cmap,12)
color_map = viridis(np.linspace(0,1, len(delays)+2))
fig, ax0 = plt.subplots(1,1,figsize = (10,6))
ax1 = ax0.twinx()
how_many = len(delays)
for i,(x,y,z) in enumerate(zip(all_the_freq[:how_many], all_of_the_1_0[:how_many], all_of_the_TOT[:how_many])):
ax0.plot(x,y, label=delays[i], ls=':', color=color_map[i])
ax1.plot(x,z, color=color_map[i])
ax0.set_xlim(0.3,0.8)
ax1.set_xlim(0.3,0.8)
ax0.legend()
ax0.set_xlabel('eV')
fig.tight_layout()
In [16]:
together = np.stack(all_of_the_1_0, axis=0)
together2 = np.stack(all_of_the_TOT, axis=0)
pp_runs = zip(all_probes_folders_selected,delays)
together.shape, freq_dipole_eV.shape, len(list(pp_runs)), together2.shape
Out[16]:
In [17]:
vector = together
cut_ev = 250
fig, ax = plt.subplots(1,1,figsize = (2,6))
z_min, z_max = -np.abs(vector).max(), np.abs(vector).max()
y, x = np.meshgrid(freq_dipole_eV[:cut_ev],delays)
c = ax.pcolormesh(x,y, vector[:,:cut_ev], cmap='RdBu', vmin=z_min, vmax=z_max)
fig.canvas.layout.height = '900px'
fig.colorbar(c, ax=ax)
ax.set_xlabel(r'$\tau$')
ax.set_ylabel(r'eV')
fig.tight_layout()
In [18]:
vector = together2
cut_ev = 250
fig, ax = plt.subplots(1,1,figsize = (2,6))
z_min, z_max = -np.abs(vector).max(), np.abs(vector).max()
#z_min, z_max = -10,10
y_axis = freq_dipole_eV[:cut_ev]
x_axis = delays
z_axis = vector[:,:cut_ev]
y, x = np.meshgrid(freq_dipole_eV[:cut_ev],delays)
c = ax.pcolormesh(x,y, vector[:,:cut_ev], cmap='RdBu', vmin=z_min, vmax=z_max)
fig.canvas.layout.height = '900px'
fig.colorbar(c, ax=ax)
ax.set_xlabel(r'$\tau$')
ax.set_ylabel(r'eV')
# ax.set_ylim(0.3,0.8)
fig.tight_layout()