From dfa8f6d428fe49222c2ee44208d850b59a80dc40 Mon Sep 17 00:00:00 2001 From: myrta75 Date: Mon, 8 Dec 2025 17:42:49 +0000 Subject: [PATCH 01/23] new file: yambopy/nl/nl_analysis.py => abstract class for analysis of multiple runs to extract response new file: yambopy/nl/sin_analysis.py => class for analysis of multiple runs with monochromatic signal to extract response --- yambopy/nl/nl_analysis.py | 113 +++++++++++++++++++++++++++++++++++++ yambopy/nl/sin_analysis.py | 112 ++++++++++++++++++++++++++++++++++++ 2 files changed, 225 insertions(+) create mode 100644 yambopy/nl/nl_analysis.py create mode 100644 yambopy/nl/sin_analysis.py diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py new file mode 100644 index 00000000..c658cabd --- /dev/null +++ b/yambopy/nl/nl_analysis.py @@ -0,0 +1,113 @@ +# Copyright (c) 2023-2025, Claudio Attaccalite, +# Myrta Gruening, Mao Yunchen +# All rights reserved. +# +# This file is part of the yambopy project +# Calculate linear response from real-time calculations (yambo_nl) +# Modified to allow generalisation +# +import numpy as np +from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 +from yambopy.nl.external_efield import Divide_by_the_Field +from tqdm import tqdm +import sys +import os +from abc import ABC,abstractmethod +# +# +# Template class +# +class Xn_from_signal(ABC): + def __init__(self,nldb,X_order=4,T_urange=[-1, -1],l_out_current=False): # note I need to add user opportunity to set time-range + self.time = nldb.IO_TIME_points # Time series + self.T_step =self.time[1] - self.time[0] # Time step of the simulation + self.efield = nldb.Efield[0] # External field of the first run + self.efields = nldb.Efield + self.n_runs = len(nldb.Polarization) # Number of external laser frequencies + self.polarization = nldb.Polarization # Array of polarizations for each laser frequency + self.current =nldb.Current # Array of currents for each laser frequency + self.l_eval_current=nldb.l_eval_CURRENT + self.l_out_current = l_out_current and self.l_eval_current + self.X_order = X_order + self.T_urange = T_urange + self.freqs = np.array([efield["freq_range"][0] for efield in nldb.Efield], dtype=np.double) # Harmonic frequencies + self.prefix = f'-{nldb.calc}' if nldb.calc != 'SAVE' else '' + super().__init__() + + def __str__(self): + """ + Print info of the class + """ + s="\n * * * Xn from signal class * * * \n\n" + s+="Max time: "+str(self.time[-1])+"\n" + s+="Time step : "+str(self.T_step)+"\n" + s+="Type Efield : "+str(self.efield["name"])+"\n" + s+="Number of runs : "+str(self.n_runs)+"\n" + s+="Current evaluated : "+str(self.l_eval_current)+"\n" + s+="Max harmonic order : "+str(self.X_order)+"\n" + if self.T_urange!=[-1, -1]: + s+="User time range : "+str(self.T_urange)+"\n" + s+="Frequency range: ["+str(self.freqs[0])+","+str(self.freqs[-1])+"] [au] \n" + if self.l_out_current: + s+="Output is set to current." + else: + s+="Output is set to polarization." + return s + + @abstractmethod + def get_sampling(self,idir,ifrq): + pass + + @abstractmethod + def define_matrix(self,samp,ifrq): + pass + + @abstractmethod + def update_time_range(self,T_period): + pass + + @abstractmethod + def get_Unit_of_Measure(self,i_order): + pass + + def solve_lin_system(self,mat,samp,solver): + mat_dim = int(mat.shape[1]) + out=np.zeros(mat_dim,dtype=np.cdouble) + SOLVE_MODES = ['stnd', 'lstsq', 'svd'] + if solver not in SOLVE_MODES: + raise ValueError("Invalid inversion mode. Expected one of: %s" % SOLVE_MODES) + if solver=="stdn" and ((not IsSquare(mat)) or (not IsWellDefined(mat))): + solver = "lstsq" + if solver=="stnd": + out = np.linalg.solve(mat,samp) + if solver=="lstsq": + out = np.linalg.lstsq(mat,samp,rcond=tol)[0] + if solver=="svd": + inv = np.linalg.pinv(M,rcond=tol) + for i_n in range(mat_dim): + out[i_n]=out[i_n]+np.sum(inv[i_n,:]*samp[:]) + return out + + def perform_analysis(self,solver="stnd"): + out = np.zeros((self.X_order + 1, self.n_runs, 3), dtype=np.cdouble) + for i_f in tqdm(range(self.n_runs)): + for i_d in range(3): + samp_time, samp_sig= self.get_sampling(i_d,i_f) + matrix = self.define_matrix(samp_time,i_f) + raw = self.solve_lin_system(matrix,samp_sig,solver) + out[:, i_f, i_d] = raw[:self.X_order + 1] + return out + + @abstractmethod + def output_analysis(self,out,to_file): + pass + + @abstractmethod + def reconstruct_signal(self,out,to_file): + pass + + def IsSquare(self,m): + return m.shape[0] == m.shape[1] + + def IsWellConditioned(self,m): # with this I am trying to avoid inverting a matrix + return np.linalg.cond(m) < 1/sys.float_info.epsilon diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py new file mode 100644 index 00000000..77be31a8 --- /dev/null +++ b/yambopy/nl/sin_analysis.py @@ -0,0 +1,112 @@ +# Copyright (c) 2023-2025, Claudio Attaccalite, +# Myrta Gruening, Mao Yunchen +# All rights reserved. +# +# This file is part of the yambopy project +# Calculate linear response from real-time calculations (yambo_nl) +# Modified to allow generalisation +# +import numpy as np +from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 +from yambopy.nl.external_efield import Divide_by_the_Field +from tqdm import tqdm +import sys +import os +from abc import ABC,abstractmethod +from nl_analysis import * +# +# +# Derived class for monochromatic signal +# +class Xn_from_sine(Xn_from_signal): + def get_sampling(self,idir,ifrq): + samp_dim = 2*self.X_order + 1 + T_period = 2.0 * np.pi / self.freqs[ifrq] + T_range, out_of_bounds = self.update_time_range(T_period) + if (out_of_bounds): + print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') + print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") + i_t_start = int(np.round(T_range[0]/self.T_step)) + i_deltaT = int(np.round(T_period/self.T_step)/samp_dim) + T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(samp_dim)]) + if self.l_out_current: + S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(samp_dim)]) # **CURRENT + else: + S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(samp_dim)]) + return T_i,S_i + + def update_time_range(self,T_period): # not sure if this is a general or specific method - let it here for the moment + T_range = self.T_urange + out_of_bounds = False + if T_range[0] <= 0.0: + T_range[0] = self.time[-1] - T_period + if T_range[1] > 0.0: + T_range[1] = T_range[0] + T_period + else: + T_range[1] = self.time[-1] + if T_range[1] > self.time[-1]: + T_range[1] = slef.time[-1] + T_range[0] = T_range[1] - T_period + out_of_bound = True + return T_range, out_of_bounds + + def define_matrix(self,T_i,ifrq): + M_size = len(T_i) + M = np.zeros((M_size, M_size), dtype=np.cdouble) + M[:, 0] = 1.0 + W = self.freqs[ifrq] + for i_n in range(1, self.X_order+1): + exp_neg = np.exp(-1j * i_n* W * T_i, dtype=np.cdouble) + exp_pos = np.exp(1j * i_n * W * T_i, dtype=np.cdouble) + M[:, i_n] = exp_neg + M[:, i_n +self.X_order] = exp_pos + return M + + def get_Unit_of_Measure(self,i_order): # not sure if this is a general or specific method - let it here for the moment + linear = 1.0 + ratio = SVCMm12VMm1 / AU2VMm1 # is there a better way than this? + if self.l_out_current: + linear = Junit/EFunit + ratio = 1.0/EFunit + if i_order == 0: + return np.power(ratio, 1, dtype=np.float64)*linear + return np.power(ratio, i_order - 1, dtype=np.float64)*linear + + def output_analysis(self,out,to_file=True): + for i_order in range(self.X_order + 1): + for i_f in range(self.n_runs): + out[i_order, i_f, :] *= Divide_by_the_Field(self.efields[i_f], i_order) + out[i_order,:,:]*=self.get_Unit_of_Measure(i_order) + if (to_file): + output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_order}' + header = "[eV] " + " ".join([f"X/Im(z){i_order} X/Re(z){i_order}" for _ in range(3)]) + if self.l_out_current: + output_file = f'o{self.prefix}.YamboPy-Sigma_probe_order_{i_order}' + header = "[eV] " + " ".join([f"S/Im(z){i_order} S/Re(z){i_order}" for _ in range(3)]) + values = np.column_stack((self.freqs * ha2ev, out[i_order, :, 0].imag, out[i_order, :, 0].real, + out[i_order, :, 1].imag, out[i_order, :, 1].real, + out[i_order, :, 2].imag, out[i_order, :, 2].real)) + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Harmonic analysis results") + else: + return (self.freqs, out) + + def reconstruct_signal(self,out,to_file=True): + Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) + for i_f in tqdm(range(self.n_runs)): + for i_d in range(3): + for i_order in range(self.X_order + 1): + freq_term = np.exp(-1j * i_order * self.freqs[i_f] * self.time) + Seff[i_f, i_d, :] += out[i_order, i_f, i_d] * freq_term + Seff[i_f, i_d, :] += np.conj(out[i_order, i_f, i_d]) * np.conj(freq_term) + if (to_file): + for i_f in tqdm(range(self.n_runs)): + values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) + output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' + header="[fs] Px Py Pz" + if self.l_out_current: + output_file = f'o{self.prefix}.YamboPy-curr_reconstructed_F{i_f + 1}' + header="[fs] Jx Jy Jz" + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") + else: + return values + From 6cfb7c1d65fd542ee40d3340e333ae23b97481a1 Mon Sep 17 00:00:00 2001 From: myrta75 Date: Thu, 11 Dec 2025 15:54:25 +0000 Subject: [PATCH 02/23] Changes to be committed: modified: yambopy/__init__.py=> freqmix_analysis added new file: yambopy/nl/freqmix_analysis.py=> analysis of two monochromatic signals (sum_frequencies) modified: yambopy/nl/nl_analysis.py=> changes to further generalise modified: yambopy/nl/sin_analysis.py=> changes as consequence of freqmix_analysis addition --- yambopy/__init__.py | 4 +- yambopy/nl/freqmix_analysis.py | 131 +++++++++++++++++++++++++++++++++ yambopy/nl/nl_analysis.py | 70 +++++++++++++----- yambopy/nl/sin_analysis.py | 39 +++++----- 4 files changed, 205 insertions(+), 39 deletions(-) create mode 100644 yambopy/nl/freqmix_analysis.py diff --git a/yambopy/__init__.py b/yambopy/__init__.py index 5cf68556..be033b6b 100644 --- a/yambopy/__init__.py +++ b/yambopy/__init__.py @@ -119,7 +119,9 @@ from yambopy.nl.harmonic_analysis import * from yambopy.nl.sum_frequencies import * from yambopy.nl.hhg_tools import * - +from yambopy.nl.nl_analysis import * +from yambopy.nl.sin_analysis import * +from yambopy.nl.freqmix_analysis import * #doublegrid files from yambopy.double_grid.dg_convergence import * diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py new file mode 100644 index 00000000..638a5662 --- /dev/null +++ b/yambopy/nl/freqmix_analysis.py @@ -0,0 +1,131 @@ +# Copyright (c) 2023-2025, Claudio Attaccalite,Mike Nico Pionteck +# Myrta Gruening +# All rights reserved. +# +# This file is part of the yambopy project +# Calculate linear response from real-time calculations (yambo_nl) +# Modified to allow generalisation +# +import numpy as np +from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 +from yambopy.nl.external_efield import Divide_by_the_Field +from tqdm import tqdm +import sys +import os +from abc import ABC,abstractmethod +from yambopy.nl.nl_analysis import Xn_from_signal +# +# +# Derived class for frequency-mixed signal +# +class Xn_from_freqmix(Xn_from_signal): + # + def set_defaults(self): + self.l_out_current = False # there is no implementation for current yet though it can be easily introduced + if self.solver == '': + self.solver = 'svd' + if self.samp_mode == '': + self.samp_mode = 'log' + EFIELDS = ["SIN","SOFTSIN"] + if self.efield["name"] not in EFIELDS: + raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") + EFIELDS = ["","SIN","SOFTSIN"] + if self.pumps[0]["name"] not in EFIELDS: + raise ValueError(f"Invalid pump field for frequency mixing analysis. Expected one of: {EFIELDS}") + self.pump_freq = 0.0 + if self.pumps[0]["name"] != "": + self.pump_freq = self.pumps[0]["freq_range"][0] + if self.pumps[1]["name"] != "": + raise ValueError("Only mixing of two fields implemented.") + if isinstance (self.X_order,int): # for frequency mixing it expects 2 orders, if one is given the same is given + self.X_order = [self.X_order,self.X_order] + if self.pump_freq == 0.0: + self.X_order[1] = 0 + self.out_dim = (2*self.X_order[0] + 1)*(2*self.X_order[1] + 1) + # + def update_time_range(self): + T_range = self.T_urange + if T_range[0] <= 0.0: + T_range[0] = self.T_deph + if T_range[1] <= 0.0: + T_range[1]=time[-1] + return T_range + # + def get_sampling(self,idir,ifrq): + samp_order = self.out_dim + if self.nsamp == -1 or self.nsamp < samp_order: + self.nsamp = int(samp_order**2) + # + T_range = update_time_range() + i_t_start = int(np.round( T_range[0] / T_step)) + i_deltaT = int(np.round((T_range[1]-T_range[0])/T_step)/self.nsamp) + if self.samp_mode=='linear': + T_i = (i_t_start + i_deltaT * np.arange(self.nsamp))*T_step + P_i = P[i_t_start + i_deltaT * np.arange(self.nsamp)] + elif self.samp_mode=='log': + T_i = np.geomspace(i_t_start * T_step, T_range[1], self.nsamp, endpoint=False) + for i in range(len(T_i)): + i_t=int(np.round(T_i[i]/T_step)) + T_i[i]=i_t*T_step + P_i[i]=P[i_t] + elif self.samp_mode=='random': + T_i = np.random.uniform(i_t_start * T_step, T_range[1], self.nsamp) + P_i = [P[int(np.round(t / T_step))] for t in T_i] + return T_i,P_i + + def define_matrix(self,T_i,ifrq): + NX,MX = self.X_order[:] + W1 = self.freqs(ifrq) + W2 = self.pump_freq + M_size = self.output_dim + M = np.zeros((self.nsamp, M_size), dtype=np.cdouble) + for i_t in range(self.nsamp): + for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): + M[i_t, i_c] = np.exp(-1j * (i_n*W1+i_n2*W2) * T_i[i_t],dtype=np.cdouble) + return M + + def output_analysis(self,out,to_file=True): + # + NX,MX = self.X_order[:] + # + for i_f in range(self.n_runs): + for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): + field_fac = 1.0 + if i_n !=0: field_fac *= Divide_by_the_Field(self.efields[i_f],abs(i_n)) #check this + if self.pump_freq !=0: + if i_m !=0: field_fac *= Divide_by_the_Field(self.pumps,abs(i_m)) #check this! what is nldb.Efield2? + out[i_c,i_f,:] *= field_fac + out[i_c,i_f,:] *= self.get_Unit_of_Measure(abs(i_n)+abs(i_m)) + if (to_file): + for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): + output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_n}_{i_m}' + iorder = (i_n-i_m,i_n,i_m) + header = "E[eV] " + " ".join([f"X{i_order}/Im_{d} X{i_order}/Re_{d}" for d in ('x','y','z')]) + values = np.column_stack((self.freqs * ha2ev, out[i_c, :, 0].imag, out[i_c, :, 0].real, + out[i_c, :, 1].imag, out[i_c, :, 1].real, + out[i_c, :, 2].imag, out[i_c, :, 2].real)) + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Frequency mixing analysis with pump frequency {self.pump_freq*ha2ev}") + else: + return (self.freqs, out) + + def reconstruct_signal(self,out,to_file=True): + Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) + for i_f in tqdm(range(self.n_runs)): + for i_d in range(3): + for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): + freq_term = np.exp(-1j * (i_n * self.freqs[i_f] + i_m*self.pump_freq) * self.time) + Seff[i_f, i_d, :] += out[i_c, i_f, i_d] * freq_term + if (to_file): + for i_f in tqdm(range(self.n_runs)): + values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) + output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' + header="[fs] Px Py Pz" + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") + else: + return values + + def spike_correction(X_eff):# Response function to check for spike + return 'Not implemented yet' + + def residuals_func(x): + return 'Not implemented yet' diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index c658cabd..0bedc998 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -1,5 +1,5 @@ # Copyright (c) 2023-2025, Claudio Attaccalite, -# Myrta Gruening, Mao Yunchen +# Myrta Gruening # All rights reserved. # # This file is part of the yambopy project @@ -18,10 +18,12 @@ # Template class # class Xn_from_signal(ABC): - def __init__(self,nldb,X_order=4,T_urange=[-1, -1],l_out_current=False): # note I need to add user opportunity to set time-range + def __init__(self,nldb,X_order=4,T_urange=[-1, -1],l_out_current=False,nsamp=-1,samp_mod='',solver=''): # note I need to add user opportunity to set time-range self.time = nldb.IO_TIME_points # Time series self.T_step =self.time[1] - self.time[0] # Time step of the simulation + self.T_deph =12.0/nldb.NL_damping self.efield = nldb.Efield[0] # External field of the first run + self.pumps = nldb.Efield_general[1:] self.efields = nldb.Efield self.n_runs = len(nldb.Polarization) # Number of external laser frequencies self.polarization = nldb.Polarization # Array of polarizations for each laser frequency @@ -29,12 +31,21 @@ def __init__(self,nldb,X_order=4,T_urange=[-1, -1],l_out_current=False): # note self.l_eval_current=nldb.l_eval_CURRENT self.l_out_current = l_out_current and self.l_eval_current self.X_order = X_order + SOLVE_MODES = ['', 'stnd', 'lstsq', 'svd'] + if solver not in SOLVE_MODES: + raise ValueError("Invalid solver mode. Expected one of: %s" % SOLVE_MODES) + self.solver = solver + SAMP_MODES = {'','linear', 'log', 'random'} + if samp_mod not in SAMP_MODES: + raise ValueError(f"Invalid sampling mode. Expected one of: {SAMP_MODES}") + self.samp_mod = samp_mod + self.nsamp = nsamp self.T_urange = T_urange self.freqs = np.array([efield["freq_range"][0] for efield in nldb.Efield], dtype=np.double) # Harmonic frequencies self.prefix = f'-{nldb.calc}' if nldb.calc != 'SAVE' else '' super().__init__() - def __str__(self): + def __str__(self): """ Print info of the class """ @@ -45,6 +56,12 @@ def __str__(self): s+="Number of runs : "+str(self.n_runs)+"\n" s+="Current evaluated : "+str(self.l_eval_current)+"\n" s+="Max harmonic order : "+str(self.X_order)+"\n" + if self.samp_mod != '': + s+="Sampling mode : "+str(self.samp_mode)+"\n" + if self.solver != '': + s+="Solver : "+str(self.solver)+"\n" + if self.nsamp > 0: + s+="Sampling points : "+str(self.nsamp)+"\n" if self.T_urange!=[-1, -1]: s+="User time range : "+str(self.T_urange)+"\n" s+="Frequency range: ["+str(self.freqs[0])+","+str(self.freqs[-1])+"] [au] \n" @@ -54,6 +71,10 @@ def __str__(self): s+="Output is set to polarization." return s + @abstractmethod + def set_defaults(self): + pass + @abstractmethod def get_sampling(self,idir,ifrq): pass @@ -63,41 +84,50 @@ def define_matrix(self,samp,ifrq): pass @abstractmethod - def update_time_range(self,T_period): + def update_time_range(self): pass @abstractmethod def get_Unit_of_Measure(self,i_order): pass - def solve_lin_system(self,mat,samp,solver): + def solve_lin_system(self,mat,samp): mat_dim = int(mat.shape[1]) out=np.zeros(mat_dim,dtype=np.cdouble) - SOLVE_MODES = ['stnd', 'lstsq', 'svd'] - if solver not in SOLVE_MODES: - raise ValueError("Invalid inversion mode. Expected one of: %s" % SOLVE_MODES) - if solver=="stdn" and ((not IsSquare(mat)) or (not IsWellDefined(mat))): - solver = "lstsq" - if solver=="stnd": + if self.solver=="stdn" and ((not IsSquare(mat)) or (not IsWellDefined(mat))): + print('WARNING: solver changed to least square') + self.solver = "lstsq" + if self.solver=="stnd": out = np.linalg.solve(mat,samp) - if solver=="lstsq": + if self.solver=="lstsq": out = np.linalg.lstsq(mat,samp,rcond=tol)[0] - if solver=="svd": - inv = np.linalg.pinv(M,rcond=tol) + if self.solver=="svd": + inv = np.linalg.pinv(mat,rcond=tol) for i_n in range(mat_dim): out[i_n]=out[i_n]+np.sum(inv[i_n,:]*samp[:]) return out - def perform_analysis(self,solver="stnd"): - out = np.zeros((self.X_order + 1, self.n_runs, 3), dtype=np.cdouble) + def perform_analysis(self): + self.set_defaults() + out = np.zeros((self.out_dim, self.n_runs, 3), dtype=np.cdouble) for i_f in tqdm(range(self.n_runs)): for i_d in range(3): samp_time, samp_sig= self.get_sampling(i_d,i_f) matrix = self.define_matrix(samp_time,i_f) - raw = self.solve_lin_system(matrix,samp_sig,solver) + raw = self.solve_lin_system(matrix,samp_sig) out[:, i_f, i_d] = raw[:self.X_order + 1] return out + def get_Unit_of_Measure(self,i_order): # not sure if this is a general or specific method - let it here for the moment + linear = 1.0 + ratio = SVCMm12VMm1 / AU2VMm1 # is there a better way than this? + if self.l_out_current: + linear = Junit/EFunit + ratio = 1.0/EFunit + if i_order == 0: + return np.power(ratio, 1, dtype=np.float64)*linear + return np.power(ratio, i_order - 1, dtype=np.float64)*linear + @abstractmethod def output_analysis(self,out,to_file): pass @@ -105,9 +135,9 @@ def output_analysis(self,out,to_file): @abstractmethod def reconstruct_signal(self,out,to_file): pass - - def IsSquare(self,m): +### some maths auxiliary functions + def IsSquare(m): return m.shape[0] == m.shape[1] - def IsWellConditioned(self,m): # with this I am trying to avoid inverting a matrix + def IsWellConditioned(m): # with this I am trying to avoid inverting a matrix return np.linalg.cond(m) < 1/sys.float_info.epsilon diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index 77be31a8..ccfed66c 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -13,26 +13,39 @@ import sys import os from abc import ABC,abstractmethod -from nl_analysis import * +from yambopy.nl.nl_analysis import Xn_from_signal # # # Derived class for monochromatic signal # class Xn_from_sine(Xn_from_signal): + @property + def set_defaults(self): + EFIELDS = ["SIN","SOFTSIN"] + if self.efield["name"] not in EFIELDS: + raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") + if any(name != '' for name in self.pumps[:]["name"]): + raise ValueError("This analysis is for one monochromatic field only.") + if self.solver == '': + self.solver = 'stnd' + self.out_dim = self.X_order + 1 + def get_sampling(self,idir,ifrq): - samp_dim = 2*self.X_order + 1 + samp_order = 2*self.X_order + 1 + if self.nsamp == -1: + self.nsamp = samp_order T_period = 2.0 * np.pi / self.freqs[ifrq] T_range, out_of_bounds = self.update_time_range(T_period) if (out_of_bounds): print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") i_t_start = int(np.round(T_range[0]/self.T_step)) - i_deltaT = int(np.round(T_period/self.T_step)/samp_dim) - T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(samp_dim)]) + i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) + T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) if self.l_out_current: - S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(samp_dim)]) # **CURRENT + S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) # **CURRENT else: - S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(samp_dim)]) + S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) return T_i,S_i def update_time_range(self,T_period): # not sure if this is a general or specific method - let it here for the moment @@ -62,16 +75,6 @@ def define_matrix(self,T_i,ifrq): M[:, i_n +self.X_order] = exp_pos return M - def get_Unit_of_Measure(self,i_order): # not sure if this is a general or specific method - let it here for the moment - linear = 1.0 - ratio = SVCMm12VMm1 / AU2VMm1 # is there a better way than this? - if self.l_out_current: - linear = Junit/EFunit - ratio = 1.0/EFunit - if i_order == 0: - return np.power(ratio, 1, dtype=np.float64)*linear - return np.power(ratio, i_order - 1, dtype=np.float64)*linear - def output_analysis(self,out,to_file=True): for i_order in range(self.X_order + 1): for i_f in range(self.n_runs): @@ -79,10 +82,10 @@ def output_analysis(self,out,to_file=True): out[i_order,:,:]*=self.get_Unit_of_Measure(i_order) if (to_file): output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_order}' - header = "[eV] " + " ".join([f"X/Im(z){i_order} X/Re(z){i_order}" for _ in range(3)]) + header = "E[eV] " + " ".join([f"X{i_order}/Im({d}) X{i_order}/Re({d})" for d in ('x','y','z')]) if self.l_out_current: output_file = f'o{self.prefix}.YamboPy-Sigma_probe_order_{i_order}' - header = "[eV] " + " ".join([f"S/Im(z){i_order} S/Re(z){i_order}" for _ in range(3)]) + header = "E[eV] " + " ".join([f"S{i_order}/Im({d}) S{i_order}/Re({d})" for d in ('x','y','z')]) values = np.column_stack((self.freqs * ha2ev, out[i_order, :, 0].imag, out[i_order, :, 0].real, out[i_order, :, 1].imag, out[i_order, :, 1].real, out[i_order, :, 2].imag, out[i_order, :, 2].real)) From 414b88c5049871fe1eb8989f64813761aac827b2 Mon Sep 17 00:00:00 2001 From: myrta75 Date: Thu, 11 Dec 2025 18:10:59 +0000 Subject: [PATCH 03/23] modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/nl_analysis.py modified: yambopy/nl/sin_analysis.py changes result from debug on a test case. --- yambopy/nl/freqmix_analysis.py | 89 +++++++++++++++++----------------- yambopy/nl/nl_analysis.py | 14 +++--- yambopy/nl/sin_analysis.py | 7 +-- 3 files changed, 57 insertions(+), 53 deletions(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index 638a5662..ebe4c081 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -12,10 +12,12 @@ from tqdm import tqdm import sys import os +import itertools from abc import ABC,abstractmethod from yambopy.nl.nl_analysis import Xn_from_signal # # +# # Derived class for frequency-mixed signal # class Xn_from_freqmix(Xn_from_signal): @@ -24,8 +26,8 @@ def set_defaults(self): self.l_out_current = False # there is no implementation for current yet though it can be easily introduced if self.solver == '': self.solver = 'svd' - if self.samp_mode == '': - self.samp_mode = 'log' + if self.samp_mod == '': + self.samp_mod = 'log' EFIELDS = ["SIN","SOFTSIN"] if self.efield["name"] not in EFIELDS: raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") @@ -33,9 +35,9 @@ def set_defaults(self): if self.pumps[0]["name"] not in EFIELDS: raise ValueError(f"Invalid pump field for frequency mixing analysis. Expected one of: {EFIELDS}") self.pump_freq = 0.0 - if self.pumps[0]["name"] != "": + if self.pumps[0]["name"] != "none": self.pump_freq = self.pumps[0]["freq_range"][0] - if self.pumps[1]["name"] != "": + if len(self.pumps)>1 and self.pumps[1]["name"] != "none": raise ValueError("Only mixing of two fields implemented.") if isinstance (self.X_order,int): # for frequency mixing it expects 2 orders, if one is given the same is given self.X_order = [self.X_order,self.X_order] @@ -48,7 +50,7 @@ def update_time_range(self): if T_range[0] <= 0.0: T_range[0] = self.T_deph if T_range[1] <= 0.0: - T_range[1]=time[-1] + T_range[1]=self.time[-1] return T_range # def get_sampling(self,idir,ifrq): @@ -56,32 +58,31 @@ def get_sampling(self,idir,ifrq): if self.nsamp == -1 or self.nsamp < samp_order: self.nsamp = int(samp_order**2) # - T_range = update_time_range() - i_t_start = int(np.round( T_range[0] / T_step)) - i_deltaT = int(np.round((T_range[1]-T_range[0])/T_step)/self.nsamp) - if self.samp_mode=='linear': - T_i = (i_t_start + i_deltaT * np.arange(self.nsamp))*T_step - P_i = P[i_t_start + i_deltaT * np.arange(self.nsamp)] - elif self.samp_mode=='log': - T_i = np.geomspace(i_t_start * T_step, T_range[1], self.nsamp, endpoint=False) - for i in range(len(T_i)): - i_t=int(np.round(T_i[i]/T_step)) - T_i[i]=i_t*T_step - P_i[i]=P[i_t] - elif self.samp_mode=='random': - T_i = np.random.uniform(i_t_start * T_step, T_range[1], self.nsamp) - P_i = [P[int(np.round(t / T_step))] for t in T_i] + T_range = self.update_time_range() + i_t_start = int(np.round( T_range[0] / self.T_step)) + i_deltaT = int(np.round((T_range[1]-T_range[0])/self.T_step)/self.nsamp) + if self.samp_mod=='linear': + T_i = (i_t_start + i_deltaT * np.arange(self.nsamp))*self.T_step + P_i = self.polarization[ifrq][idir,i_t_start + i_deltaT * np.arange(self.nsamp)] + elif self.samp_mod=='log': + T_i = np.geomspace(i_t_start * self.T_step, T_range[1], self.nsamp, endpoint=False) + i_t = np.array(np.round(T_i/self.T_step),dtype=int) + T_i = i_t*self.T_step + P_i = np.array([self.polarization[ifrq][idir,i] for i in i_t]) + elif self.samp_mod=='random': + T_i = np.random.uniform(i_t_start * self.T_step, T_range[1], self.nsamp) + P_i = np.array([self.polarization[ifrq][idir,int(np.round(t / self.T_step))] for t in T_i]) return T_i,P_i def define_matrix(self,T_i,ifrq): NX,MX = self.X_order[:] - W1 = self.freqs(ifrq) + W1 = self.freqs[ifrq] W2 = self.pump_freq - M_size = self.output_dim + M_size = self.out_dim M = np.zeros((self.nsamp, M_size), dtype=np.cdouble) for i_t in range(self.nsamp): for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): - M[i_t, i_c] = np.exp(-1j * (i_n*W1+i_n2*W2) * T_i[i_t],dtype=np.cdouble) + M[i_t, i_c] = np.exp(-1j * (i_n*W1+i_m*W2) * T_i[i_t],dtype=np.cdouble) return M def output_analysis(self,out,to_file=True): @@ -93,39 +94,39 @@ def output_analysis(self,out,to_file=True): field_fac = 1.0 if i_n !=0: field_fac *= Divide_by_the_Field(self.efields[i_f],abs(i_n)) #check this if self.pump_freq !=0: - if i_m !=0: field_fac *= Divide_by_the_Field(self.pumps,abs(i_m)) #check this! what is nldb.Efield2? + if i_m !=0: field_fac *= Divide_by_the_Field(self.pumps[0],abs(i_m)) #check this! what is nldb.Efield2? out[i_c,i_f,:] *= field_fac out[i_c,i_f,:] *= self.get_Unit_of_Measure(abs(i_n)+abs(i_m)) if (to_file): for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_n}_{i_m}' iorder = (i_n-i_m,i_n,i_m) - header = "E[eV] " + " ".join([f"X{i_order}/Im_{d} X{i_order}/Re_{d}" for d in ('x','y','z')]) + header = "E[eV] " + " ".join([f"X{iorder}/Im_{d} X{iorder}/Re_{d}" for d in ('x','y','z')]) values = np.column_stack((self.freqs * ha2ev, out[i_c, :, 0].imag, out[i_c, :, 0].real, out[i_c, :, 1].imag, out[i_c, :, 1].real, out[i_c, :, 2].imag, out[i_c, :, 2].real)) - np.savetxt(output_file, values, header=header, delimiter=' ', footer="Frequency mixing analysis with pump frequency {self.pump_freq*ha2ev}") + np.savetxt(output_file, values, header=header, delimiter=' ', footer=f"Frequency mixing analysis with pump frequency {self.pump_freq*ha2ev} eV") else: return (self.freqs, out) - def reconstruct_signal(self,out,to_file=True): - Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) + def reconstruct_signal(self,out,to_file=True): + Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) + for i_f in tqdm(range(self.n_runs)): + for i_d in range(3): + for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): + freq_term = np.exp(-1j * (i_n * self.freqs[i_f] + i_m*self.pump_freq) * self.time) + Seff[i_f, i_d, :] += out[i_c, i_f, i_d] * freq_term + if (to_file): for i_f in tqdm(range(self.n_runs)): - for i_d in range(3): - for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): - freq_term = np.exp(-1j * (i_n * self.freqs[i_f] + i_m*self.pump_freq) * self.time) - Seff[i_f, i_d, :] += out[i_c, i_f, i_d] * freq_term - if (to_file): - for i_f in tqdm(range(self.n_runs)): - values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) - output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' - header="[fs] Px Py Pz" - np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") - else: - return values + values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) + output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' + header="[fs] Px Py Pz" + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") + else: + return values - def spike_correction(X_eff):# Response function to check for spike - return 'Not implemented yet' + #def spike_correction(X_eff):# Response function to check for spike + # return 'Not implemented yet' - def residuals_func(x): - return 'Not implemented yet' + #def residuals_func(x): + # return 'Not implemented yet' diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index 0bedc998..2ceb3761 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -18,7 +18,7 @@ # Template class # class Xn_from_signal(ABC): - def __init__(self,nldb,X_order=4,T_urange=[-1, -1],l_out_current=False,nsamp=-1,samp_mod='',solver=''): # note I need to add user opportunity to set time-range + def __init__(self,nldb,X_order=4,T_range=[-1, -1],l_out_current=False,nsamp=-1,samp_mod='',solver='',tol=1e-10): # note I need to add user opportunity to set time-range self.time = nldb.IO_TIME_points # Time series self.T_step =self.time[1] - self.time[0] # Time step of the simulation self.T_deph =12.0/nldb.NL_damping @@ -40,9 +40,11 @@ def __init__(self,nldb,X_order=4,T_urange=[-1, -1],l_out_current=False,nsamp=-1, raise ValueError(f"Invalid sampling mode. Expected one of: {SAMP_MODES}") self.samp_mod = samp_mod self.nsamp = nsamp - self.T_urange = T_urange + self.T_urange = T_range self.freqs = np.array([efield["freq_range"][0] for efield in nldb.Efield], dtype=np.double) # Harmonic frequencies self.prefix = f'-{nldb.calc}' if nldb.calc != 'SAVE' else '' + self.out_dim = 0 + self.tol = tol super().__init__() def __str__(self): @@ -100,22 +102,22 @@ def solve_lin_system(self,mat,samp): if self.solver=="stnd": out = np.linalg.solve(mat,samp) if self.solver=="lstsq": - out = np.linalg.lstsq(mat,samp,rcond=tol)[0] + out = np.linalg.lstsq(mat,samp,rcond=self.tol)[0] if self.solver=="svd": - inv = np.linalg.pinv(mat,rcond=tol) + inv = np.linalg.pinv(mat,rcond=self.tol) for i_n in range(mat_dim): out[i_n]=out[i_n]+np.sum(inv[i_n,:]*samp[:]) return out def perform_analysis(self): - self.set_defaults() + _ = self.set_defaults() out = np.zeros((self.out_dim, self.n_runs, 3), dtype=np.cdouble) for i_f in tqdm(range(self.n_runs)): for i_d in range(3): samp_time, samp_sig= self.get_sampling(i_d,i_f) matrix = self.define_matrix(samp_time,i_f) raw = self.solve_lin_system(matrix,samp_sig) - out[:, i_f, i_d] = raw[:self.X_order + 1] + out[:, i_f, i_d] = raw[:self.out_dim] return out def get_Unit_of_Measure(self,i_order): # not sure if this is a general or specific method - let it here for the moment diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index ccfed66c..e7738af2 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -19,16 +19,17 @@ # Derived class for monochromatic signal # class Xn_from_sine(Xn_from_signal): - @property def set_defaults(self): EFIELDS = ["SIN","SOFTSIN"] if self.efield["name"] not in EFIELDS: raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") - if any(name != '' for name in self.pumps[:]["name"]): - raise ValueError("This analysis is for one monochromatic field only.") + for i_n in range(len(self.pumps)): + if self.pumps[i_n]["name"] != 'none': + raise ValueError("This analysis is for one monochromatic field only.") if self.solver == '': self.solver = 'stnd' self.out_dim = self.X_order + 1 + return def get_sampling(self,idir,ifrq): samp_order = 2*self.X_order + 1 From 3f6f80cfdbb1cbf6133e10ed2c733884a0d3ecf9 Mon Sep 17 00:00:00 2001 From: Myrta Date: Thu, 11 Dec 2025 18:42:17 +0000 Subject: [PATCH 04/23] Changes to be committed: modified: yambopy/nl/freqmix_analysis.py Changed header of output files --- yambopy/nl/freqmix_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index ebe4c081..ab05cbba 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -100,7 +100,7 @@ def output_analysis(self,out,to_file=True): if (to_file): for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_n}_{i_m}' - iorder = (i_n-i_m,i_n,i_m) + iorder = (i_n+i_m,i_n,i_m) header = "E[eV] " + " ".join([f"X{iorder}/Im_{d} X{iorder}/Re_{d}" for d in ('x','y','z')]) values = np.column_stack((self.freqs * ha2ev, out[i_c, :, 0].imag, out[i_c, :, 0].real, out[i_c, :, 1].imag, out[i_c, :, 1].real, From 92d02e856110c9f3b5d3d31917e28cb48da2e30e Mon Sep 17 00:00:00 2001 From: Myrta Date: Fri, 12 Dec 2025 11:45:12 +0000 Subject: [PATCH 05/23] Changes to be committed: modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/nl_analysis.py modified: yambopy/nl/sin_analysis.py Added an extra lstsq solver based on the scipy lstsq optimizer which takes initial values --- yambopy/nl/freqmix_analysis.py | 3 --- yambopy/nl/nl_analysis.py | 23 ++++++++++++++++++----- yambopy/nl/sin_analysis.py | 2 +- 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index ab05cbba..e3239fe8 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -127,6 +127,3 @@ def reconstruct_signal(self,out,to_file=True): #def spike_correction(X_eff):# Response function to check for spike # return 'Not implemented yet' - - #def residuals_func(x): - # return 'Not implemented yet' diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index 2ceb3761..0b199214 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -10,6 +10,7 @@ from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 from yambopy.nl.external_efield import Divide_by_the_Field from tqdm import tqdm +from scipy.optimize import least_squares import sys import os from abc import ABC,abstractmethod @@ -31,7 +32,7 @@ def __init__(self,nldb,X_order=4,T_range=[-1, -1],l_out_current=False,nsamp=-1,s self.l_eval_current=nldb.l_eval_CURRENT self.l_out_current = l_out_current and self.l_eval_current self.X_order = X_order - SOLVE_MODES = ['', 'stnd', 'lstsq', 'svd'] + SOLVE_MODES = ['', 'full', 'lstsq', 'lstsq_opt', 'svd'] if solver not in SOLVE_MODES: raise ValueError("Invalid solver mode. Expected one of: %s" % SOLVE_MODES) self.solver = solver @@ -93,13 +94,13 @@ def update_time_range(self): def get_Unit_of_Measure(self,i_order): pass - def solve_lin_system(self,mat,samp): + def solve_lin_system(self,mat,samp,init=None): mat_dim = int(mat.shape[1]) out=np.zeros(mat_dim,dtype=np.cdouble) - if self.solver=="stdn" and ((not IsSquare(mat)) or (not IsWellDefined(mat))): + if self.solver=="full" and ((not IsSquare(mat)) or (not IsWellDefined(mat))): print('WARNING: solver changed to least square') self.solver = "lstsq" - if self.solver=="stnd": + if self.solver=="full": out = np.linalg.solve(mat,samp) if self.solver=="lstsq": out = np.linalg.lstsq(mat,samp,rcond=self.tol)[0] @@ -107,6 +108,14 @@ def solve_lin_system(self,mat,samp): inv = np.linalg.pinv(mat,rcond=self.tol) for i_n in range(mat_dim): out[i_n]=out[i_n]+np.sum(inv[i_n,:]*samp[:]) + if self.solver=="lstsq_opt": + if(init is None): + x0_cmplx = np.linalg.lstsq(mat, samp, rcond=tol)[0] + else: + x0_cmplx = init + x0 = np.concatenate((x0_cmplx.real, x0_cmplx.imag)) + res = least_squares(residuals_func, x0, ftol=1e-11,gtol=1e-11,xtol=1e-11,verbose=0,x_scale='jac',args=(mat,samp)) + out = res.x[0:int(res.x.size/2)] + 1j * res.x[int(res.x.size/2):res.x.size] return out def perform_analysis(self): @@ -122,7 +131,7 @@ def perform_analysis(self): def get_Unit_of_Measure(self,i_order): # not sure if this is a general or specific method - let it here for the moment linear = 1.0 - ratio = SVCMm12VMm1 / AU2VMm1 # is there a better way than this? + ratio = SVCMm12VMm1 / AU2VMm1 # From AU to statVolt/cm ...is there a better way than this? if self.l_out_current: linear = Junit/EFunit ratio = 1.0/EFunit @@ -143,3 +152,7 @@ def IsSquare(m): def IsWellConditioned(m): # with this I am trying to avoid inverting a matrix return np.linalg.cond(m) < 1/sys.float_info.epsilon + + def residuals_func(x,M,S_i): + x_cmplx=x[0:int(x.size/2)] + 1j * x[int(x.size/2):x.size] + return np.linalg.norm(np.dot(M, x_cmplx) - S_i) diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index e7738af2..65245d4b 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -27,7 +27,7 @@ def set_defaults(self): if self.pumps[i_n]["name"] != 'none': raise ValueError("This analysis is for one monochromatic field only.") if self.solver == '': - self.solver = 'stnd' + self.solver = 'full' self.out_dim = self.X_order + 1 return From 4bf1f1e7cb1e64ecbd2b77cc3e111c710d8e100a Mon Sep 17 00:00:00 2001 From: Myrta Date: Fri, 12 Dec 2025 12:29:58 +0000 Subject: [PATCH 06/23] Changes to be committed: modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/nl_analysis.py coded the spike_correction function even though it is not used so far --- yambopy/nl/freqmix_analysis.py | 31 ++++++++++++++++++++++++++++--- yambopy/nl/nl_analysis.py | 2 +- 2 files changed, 29 insertions(+), 4 deletions(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index e3239fe8..d77ddf32 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -10,6 +10,7 @@ from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 from yambopy.nl.external_efield import Divide_by_the_Field from tqdm import tqdm +from scipy.andimage import uniform_filter1d import sys import os import itertools @@ -124,6 +125,30 @@ def reconstruct_signal(self,out,to_file=True): np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") else: return values - - #def spike_correction(X_eff):# Response function to check for spike - # return 'Not implemented yet' + # + # unused function so far + # check for spikes in the solution and recalculate the solution + # The function can be called after perform_analysis + # in a loop over cartesian directions + # Note that it can be taken into nl_analysis to be shared with other classes + # + def spike_correction(self,mat,samp,out): + window_size, threshold_local = (5,2.0) + out_corr = out + self.solver = "lstq_opt" + signal_im= abs(out[self.X_order[0]+1,self.X_order[1]+1,:].imag) + signal_re= abs(out[self.X_order[0]+1,self.X_order[1]+1,:].real) + smooth_signal_im = uniform_filter1d(signal_im, size=window_size) + smooth_signal_re = uniform_filter1d(signal_re, size=window_size) + spike_indices_local_im = np.where(np.abs(signal_im - smooth_signal_im)/smooth_signal_im > threshold_local)[0] + spike_indices_local_re = np.where(np.abs(signal_re - smooth_signal_re)/smooth_signal_re > threshold_local)[0] + spike_indices_local = np.unique(np.concatenate((spike_indices_local_im, spike_indices_local_re))) + for i_f in tqdm(spike_indices_local): + if(i_f==0 and i_f in spike_indices_local): + x0=out[:,i_f+1] + elif(i_f==n_frequencies-1 and i_f in spike_indices_local): + x0=out[:,i_f-1] + else: + x0=(out[:,i_f+1]+out[:,i_f-1])/2.0 + out_corr[:,i_f] = self.solve_lin_system(mat,samp,init=x0) + return out_corr, spike_indices_local diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index 0b199214..c2f94168 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -153,6 +153,6 @@ def IsSquare(m): def IsWellConditioned(m): # with this I am trying to avoid inverting a matrix return np.linalg.cond(m) < 1/sys.float_info.epsilon - def residuals_func(x,M,S_i): + def residuals_func(x,M,S_i): x_cmplx=x[0:int(x.size/2)] + 1j * x[int(x.size/2):x.size] return np.linalg.norm(np.dot(M, x_cmplx) - S_i) From 71d301edd7bd9bb616559bf6f2f447b5f0c02c79 Mon Sep 17 00:00:00 2001 From: Myrta Date: Fri, 12 Dec 2025 16:00:09 +0000 Subject: [PATCH 07/23] Changes to be committed: new file: yambopy/nl/Xn_from_signal_diagrams.md started a doc of the code suite for nl analysis ). use the md that can be then converted in whatevers --- yambopy/nl/Xn_from_signal_diagrams.md | 133 ++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 yambopy/nl/Xn_from_signal_diagrams.md diff --git a/yambopy/nl/Xn_from_signal_diagrams.md b/yambopy/nl/Xn_from_signal_diagrams.md new file mode 100644 index 00000000..ad679b7e --- /dev/null +++ b/yambopy/nl/Xn_from_signal_diagrams.md @@ -0,0 +1,133 @@ +# `Xn_from_signal` – Mermaid Diagrams + +This document contains Mermaid diagrams derived from `nl_analysis.py` for the abstract class **`Xn_from_signal`** and its key workflows. + +--- + +## 1) Class structure (`classDiagram`) +```mermaid +classDiagram + class Xn_from_signal { + <> + %% --- Attributes --- + +time : np.ndarray + +T_step : float + +T_deph : float + +efield : dict + +pumps : list + +efields : list + +n_runs : int + +polarization : np.ndarray + +current : np.ndarray + +l_eval_current : bool + +l_out_current : bool + +X_order : int + +solver : str + +samp_mod : str + +nsamp : int + +T_urange : list[float] + +freqs : np.ndarray + +prefix : str + +out_dim : int + +tol : float + + %% --- Constructor --- + +__init__(nldb, X_order=4, T_range=[-1,-1], l_out_current=False, nsamp=-1, samp_mod='', solver='', tol=1e-10) + + %% --- Special --- + +__str__() str + + %% --- Abstract hooks --- + #set_defaults() + #get_sampling(idir, ifrq) + #define_matrix(samp, ifrq) + #update_time_range() + #get_Unit_of_Measure(i_order) + #output_analysis(out, to_file) + #reconstruct_signal(out, to_file) + + %% --- Concrete methods --- + +solve_lin_system(mat, samp, init=None) np.ndarray + +perform_analysis() np.ndarray + +get_Unit_of_Measure(i_order) float + } + + %% Auxiliary (module-level) functions for completeness + class AuxMath { + +IsSquare(m) bool + +IsWellConditioned(m) bool + +residuals_func(x, M, S_i) float + } +``` + +--- + +## 2) Workflow: `perform_analysis` (`flowchart`) +```mermaid +graph TD + A[perform_analysis] --> B[set_defaults] + B --> C[allocate_out] + C --> D{loop_i_f} + D --> E{loop_i_d} + E --> F[get_sampling] + F --> G[define_matrix] + G --> H[solve_lin_system] + H --> I[assign_out] + I --> J{more_i_d} + J -- Yes --> E + J -- No --> K{more_i_f} + K -- Yes --> D + K -- No --> L[return_out] +``` + +--- + +## 3) Workflow: `solve_lin_system` (`flowchart`) +```mermaid +graph TD + A[solve_lin_system] --> B[init_out] + B --> C{solver_full} + C -- Yes --> D{square_and_well_conditioned} + D -- No --> E[set_solver_lstsq] + D -- Yes --> F[linalg_solve] + C -- No --> G{solver_lstsq} + G -- Yes --> H[lstsq] + G -- No --> I{solver_svd} + I -- Yes --> J[pinv] + J --> K[accumulate_inv_times_samp] + I -- No --> L{solver_lstsq_opt} + L -- Yes --> M{has_init} + M -- Yes --> O[use_init] + M -- No --> N[lstsq_for_x0] + N --> P[concat_real_and_imag] + O --> P + P --> Q[least_squares] + Q --> R[compose_complex_out] + L -- No --> S[return_out] + F --> T[return_out] + H --> T + K --> T + R --> T + E --> H +``` + +--- + +## 4) Sequence view (optional) +```mermaid +sequenceDiagram + participant Analyzer as Xn_from_signal + participant Impl as Subclass (implements abstract hooks) + + Analyzer->>Impl: set_defaults() + loop for each run i_f and direction i_d + Analyzer->>Impl: get_sampling(i_d, i_f) + Impl-->>Analyzer: (samp_time, samp_sig) + Analyzer->>Impl: define_matrix(samp_time, i_f) + Impl-->>Analyzer: matrix + Analyzer->>Analyzer: solve_lin_system(matrix, samp_sig) + Analyzer->>Analyzer: out[:, i_f, i_d] = raw[:out_dim] + end + Analyzer-->>Impl: output_analysis(out, to_file) + Analyzer-->>Impl: reconstruct_signal(out, to_file) +``` From d31d6211545199f6afebd5f349d62c2515511d4a Mon Sep 17 00:00:00 2001 From: Myrta Date: Fri, 12 Dec 2025 16:24:16 +0000 Subject: [PATCH 08/23] Changes to be committed: new file: yambopy/nl/pulse_analysis.py Started the new class for pulse analysiss --- yambopy/nl/pulse_analysis.py | 116 +++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 yambopy/nl/pulse_analysis.py diff --git a/yambopy/nl/pulse_analysis.py b/yambopy/nl/pulse_analysis.py new file mode 100644 index 00000000..99f8e433 --- /dev/null +++ b/yambopy/nl/pulse_analysis.py @@ -0,0 +1,116 @@ +# Copyright (c) 2023-2025, Claudio Attaccalite, +# Myrta Gruening, Anna Romani +# All rights reserved. +# +# This file is part of the yambopy project +# Calculate linear response from real-time calculations (yambo_nl) +# Modified to allow generalisation +# +import numpy as np +from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 +from yambopy.nl.external_efield import Divide_by_the_Field +from tqdm import tqdm +import sys +import os +from abc import ABC,abstractmethod +from yambopy.nl.nl_analysis import Xn_from_signal +# +# +# Derived class for monochromatic signal +# +class Xn_from_pulse(Xn_from_signal): + def set_defaults(self): + EFIELDS = ["QSSIN"] + if self.efield["name"] not in EFIELDS: + raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") + for i_n in range(len(self.pumps)): + if self.pumps[i_n]["name"] != 'none': + raise ValueError("This analysis is for one monochromatic field only.") + if self.solver == '': + self.solver = 'full' + self.out_dim = self.X_order + 1 # CHANGE THIS + return + + # def get_sampling(self,idir,ifrq): + # samp_order = 2*self.X_order + 1 + # if self.nsamp == -1: + # self.nsamp = samp_order + # T_period = 2.0 * np.pi / self.freqs[ifrq] + # T_range, out_of_bounds = self.update_time_range(T_period) + # if (out_of_bounds): + # print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') + # print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") + # i_t_start = int(np.round(T_range[0]/self.T_step)) + # i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) + # T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) + # if self.l_out_current: + # S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) # **CURRENT + # else: + # S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) + # return T_i,S_i + + # def update_time_range(self,T_period): + # T_range = self.T_urange + # out_of_bounds = False + # if T_range[0] <= 0.0: + # T_range[0] = self.time[-1] - T_period + # if T_range[1] > 0.0: + # T_range[1] = T_range[0] + T_period + # else: + # T_range[1] = self.time[-1] + # if T_range[1] > self.time[-1]: + # T_range[1] = slef.time[-1] + # T_range[0] = T_range[1] - T_period + # out_of_bound = True + # return T_range, out_of_bounds + + # def define_matrix(self,T_i,ifrq): + # M_size = len(T_i) + # M = np.zeros((M_size, M_size), dtype=np.cdouble) + # M[:, 0] = 1.0 + # W = self.freqs[ifrq] + # for i_n in range(1, self.X_order+1): + # exp_neg = np.exp(-1j * i_n* W * T_i, dtype=np.cdouble) + # exp_pos = np.exp(1j * i_n * W * T_i, dtype=np.cdouble) + # M[:, i_n] = exp_neg + # M[:, i_n +self.X_order] = exp_pos + # return M + + # def output_analysis(self,out,to_file=True): + # for i_order in range(self.X_order + 1): + # for i_f in range(self.n_runs): + # out[i_order, i_f, :] *= Divide_by_the_Field(self.efields[i_f], i_order) + # out[i_order,:,:]*=self.get_Unit_of_Measure(i_order) + # if (to_file): + # output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_order}' + # header = "E[eV] " + " ".join([f"X{i_order}/Im({d}) X{i_order}/Re({d})" for d in ('x','y','z')]) + # if self.l_out_current: + # output_file = f'o{self.prefix}.YamboPy-Sigma_probe_order_{i_order}' + # header = "E[eV] " + " ".join([f"S{i_order}/Im({d}) S{i_order}/Re({d})" for d in ('x','y','z')]) + # values = np.column_stack((self.freqs * ha2ev, out[i_order, :, 0].imag, out[i_order, :, 0].real, + # out[i_order, :, 1].imag, out[i_order, :, 1].real, + # out[i_order, :, 2].imag, out[i_order, :, 2].real)) + # np.savetxt(output_file, values, header=header, delimiter=' ', footer="Harmonic analysis results") + # else: + # return (self.freqs, out) + + # def reconstruct_signal(self,out,to_file=True): + # Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) + # for i_f in tqdm(range(self.n_runs)): + # for i_d in range(3): + # for i_order in range(self.X_order + 1): + # freq_term = np.exp(-1j * i_order * self.freqs[i_f] * self.time) + # Seff[i_f, i_d, :] += out[i_order, i_f, i_d] * freq_term + # Seff[i_f, i_d, :] += np.conj(out[i_order, i_f, i_d]) * np.conj(freq_term) + # if (to_file): + # for i_f in tqdm(range(self.n_runs)): + # values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) + # output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' + # header="[fs] Px Py Pz" + # if self.l_out_current: + # output_file = f'o{self.prefix}.YamboPy-curr_reconstructed_F{i_f + 1}' + # header="[fs] Jx Jy Jz" + # np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") + # else: + # return values + From 7c14fda8e62a941b419336f5892dedc383bddf09 Mon Sep 17 00:00:00 2001 From: Myrta Date: Tue, 23 Dec 2025 16:57:42 +0100 Subject: [PATCH 09/23] Changes to be committed: modified: yambopy/nl/external_efield.py modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/pulse_analysis.py modified: yambopy/nl/sin_analysis.py finished to code the pulse analysis and small chenges to the other files. pulse analysis to be tested. --- yambopy/nl/external_efield.py | 29 +++-- yambopy/nl/freqmix_analysis.py | 3 +- yambopy/nl/pulse_analysis.py | 227 +++++++++++++++++++-------------- yambopy/nl/sin_analysis.py | 24 ++-- 4 files changed, 163 insertions(+), 120 deletions(-) diff --git a/yambopy/nl/external_efield.py b/yambopy/nl/external_efield.py index 336b093f..9b3c7ba3 100644 --- a/yambopy/nl/external_efield.py +++ b/yambopy/nl/external_efield.py @@ -20,22 +20,27 @@ def Divide_by_the_Field(efield,order): divide_by_field=np.power(-2.0*1.0j/efield['amplitude'],order,dtype=np.cdouble) elif order==0: divide_by_field=4.0/np.power(efield['amplitude'],2.0,dtype=np.cdouble) - - elif efield['name'] == 'QSSIN': - # Approximate relations/does not work yet - sigma=efield['width'] - W_0=efield['freq_range'][0] - T_0= np.pi/W_0*float(round(W_0/np.pi*3.*sigma)) - T = 2*np.pi/W_0 - E_w= math.sqrt(np.pi/2)*sigma*np.exp(-1j*W_0*T_0)*(special.erf((T-T_0)/math.sqrt(2.0)/sigma)+special.erf(T_0/math.sqrt(2.0)/sigma)) + # elif efield['name'] == 'QSSIN': ! note that in class Xn_from_pulse there is a factor + # + # This part of code was copied from ypp but it was never used + # + # sigma=efield['width'] + # W_0=efield['freq_range'][0] + # T_0= np.pi/W_0*float(round(W_0/np.pi*3.*sigma)) + # T = 2*np.pi/W_0 + # E_w= math.sqrt(np.pi/2)*sigma*np.exp(-1j*W_0*T_0)*(special.erf((T-T_0)/math.sqrt(2.0)/sigma)+special.erf(T_0/math.sqrt(2.0)/sigma)) - if order!=0: - divide_by_field = (-2.0*1.0j/(E_w*efield['amplitude']))**order - elif order==0: - divide_by_field = 4.0/(E_w*efield['amplitude']*np.conj(E_w)) + # if order!=0: + # divide_by_field = (-2.0*1.0j/(E_w*efield['amplitude']))**order + # elif order==0: + # divide_by_field = 4.0/(E_w*efield['amplitude']*np.conj(E_w)) else: raise ValueError("Electric field not implemented in Divide_by_the_Field!") return divide_by_field +def Gaussian_centre(efield): + ratio=np.pi/efield['freq_range'][0] + sigma=efield["damping"]/(2.0*(2.0*np.log(2.0))**0.5) + return ratio * float(round(1.0 /ratio * sigma *efield["field_peak"])) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index d77ddf32..4baa4a55 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -79,8 +79,7 @@ def define_matrix(self,T_i,ifrq): NX,MX = self.X_order[:] W1 = self.freqs[ifrq] W2 = self.pump_freq - M_size = self.out_dim - M = np.zeros((self.nsamp, M_size), dtype=np.cdouble) + M = np.zeros((self.nsamp, self.out_dim), dtype=np.cdouble) for i_t in range(self.nsamp): for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): M[i_t, i_c] = np.exp(-1j * (i_n*W1+i_m*W2) * T_i[i_t],dtype=np.cdouble) diff --git a/yambopy/nl/pulse_analysis.py b/yambopy/nl/pulse_analysis.py index 99f8e433..9a7ed41b 100644 --- a/yambopy/nl/pulse_analysis.py +++ b/yambopy/nl/pulse_analysis.py @@ -7,110 +7,149 @@ # Modified to allow generalisation # import numpy as np -from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 -from yambopy.nl.external_efield import Divide_by_the_Field +from yambopy.units import ha2ev,fs2aut +from yambopy.nl.external_efield import Divide_by_the_Field, Gaussian_centre from tqdm import tqdm +from math import floor, ceil, factorial import sys import os from abc import ABC,abstractmethod from yambopy.nl.nl_analysis import Xn_from_signal # -# -# Derived class for monochromatic signal +# Derived class for pulse signal # class Xn_from_pulse(Xn_from_signal): - def set_defaults(self): - EFIELDS = ["QSSIN"] - if self.efield["name"] not in EFIELDS: - raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") - for i_n in range(len(self.pumps)): - if self.pumps[i_n]["name"] != 'none': - raise ValueError("This analysis is for one monochromatic field only.") - if self.solver == '': - self.solver = 'full' - self.out_dim = self.X_order + 1 # CHANGE THIS - return + # + def set_defaults(self): + EFIELDS = ["QSSIN"] + if self.efield["name"] not in EFIELDS: + raise ValueError(f"Invalid electric field for frequency mixing analysis. Expected one of: {EFIELDS}") + for i_n in range(len(self.pumps)): + if self.pumps[i_n]["name"] != 'none': + raise ValueError("This analysis is for one monochromatic field only.") + if self.solver == '': + self.solver = 'full' + dim_ = int(1+self.X_order+(self.X_order/2)**2) + D = (dim_,int(dim_ -0.25)) + self.out_dim = D[self.X_order%2] + self.T0 = Gaussian_centre(self.efield) + return + + def build_map(self): + M = 1+2*self.X_order + int(self.X_order*(self.X_order-1)/2) + L = (int(self.X_order/2*((self.X_order/2+1))),int(((self.X_order+1)/2)**2)) + I = np.zeros((M,2),dtype=int) + j = 0 + I[j,:] = [0,0] + for n in range(2,N+1,2): + j += 1 + I[j,:] = [n,floor(n/2)] + for n in range(1,N+1): + for m in range(0,ceil(n/2)): + j += 1 + I[j,:] = [n,m] + I[j+L[self.X_order%2],:] = [-n,m] + return I - # def get_sampling(self,idir,ifrq): - # samp_order = 2*self.X_order + 1 - # if self.nsamp == -1: - # self.nsamp = samp_order - # T_period = 2.0 * np.pi / self.freqs[ifrq] - # T_range, out_of_bounds = self.update_time_range(T_period) - # if (out_of_bounds): - # print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') - # print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") - # i_t_start = int(np.round(T_range[0]/self.T_step)) - # i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) - # T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) - # if self.l_out_current: - # S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) # **CURRENT - # else: - # S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) - # return T_i,S_i + def get_sampling(self,idir,ifrq): + M = 1+2*self.X_order + int(self.X_order*(self.X_order-1)/2) + if self.nsamp == -1 or self.nsamp < M: + self.nsamp = M + T_period = 2.0 * np.pi / self.freqs[ifrq] + T_range, out_of_bounds = self.update_time_range(T_period) + if (out_of_bounds): + print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') + print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") + i_t_start = int(np.round(T_range[0]/self.T_step)) + i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) + T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) + S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) + return T_i,S_i - # def update_time_range(self,T_period): - # T_range = self.T_urange - # out_of_bounds = False - # if T_range[0] <= 0.0: - # T_range[0] = self.time[-1] - T_period - # if T_range[1] > 0.0: - # T_range[1] = T_range[0] + T_period - # else: - # T_range[1] = self.time[-1] - # if T_range[1] > self.time[-1]: - # T_range[1] = slef.time[-1] - # T_range[0] = T_range[1] - T_period - # out_of_bound = True - # return T_range, out_of_bounds - - # def define_matrix(self,T_i,ifrq): - # M_size = len(T_i) - # M = np.zeros((M_size, M_size), dtype=np.cdouble) - # M[:, 0] = 1.0 - # W = self.freqs[ifrq] - # for i_n in range(1, self.X_order+1): - # exp_neg = np.exp(-1j * i_n* W * T_i, dtype=np.cdouble) - # exp_pos = np.exp(1j * i_n * W * T_i, dtype=np.cdouble) - # M[:, i_n] = exp_neg - # M[:, i_n +self.X_order] = exp_pos - # return M + def update_time_range(self,T_period): + T_range = self.T_urange + out_of_bounds = False + if T_range[0] <= 0.0: + T_range[0] = self.T0 - T_period/2 + T_range[1] = T_range[0] + T_period + if T_range[1] > self.time[-1]: + T_range[1] = slef.time[-1] + T_range[0] = T_range[1] - T_period + out_of_bound = True + return T_range, out_of_bounds - # def output_analysis(self,out,to_file=True): - # for i_order in range(self.X_order + 1): - # for i_f in range(self.n_runs): - # out[i_order, i_f, :] *= Divide_by_the_Field(self.efields[i_f], i_order) - # out[i_order,:,:]*=self.get_Unit_of_Measure(i_order) - # if (to_file): - # output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_order}' - # header = "E[eV] " + " ".join([f"X{i_order}/Im({d}) X{i_order}/Re({d})" for d in ('x','y','z')]) - # if self.l_out_current: - # output_file = f'o{self.prefix}.YamboPy-Sigma_probe_order_{i_order}' - # header = "E[eV] " + " ".join([f"S{i_order}/Im({d}) S{i_order}/Re({d})" for d in ('x','y','z')]) - # values = np.column_stack((self.freqs * ha2ev, out[i_order, :, 0].imag, out[i_order, :, 0].real, - # out[i_order, :, 1].imag, out[i_order, :, 1].real, - # out[i_order, :, 2].imag, out[i_order, :, 2].real)) - # np.savetxt(output_file, values, header=header, delimiter=' ', footer="Harmonic analysis results") - # else: - # return (self.freqs, out) + def define_matrix(self,T_i,ifrq): + dim_ = 1+2*self.X_order + int(self.X_order*(self.X_order-1)/2) + W = self.freqs[ifrq] + M = np.zeros((self.nsamp, dim_), dtype=np.cdouble) + i_map = self.build_map() + M[:,0] = 1.0 + for ii in range(1,len(i_map)): + i_n,i_m = i_map[ii] + n = abs(i_n) + M[:,ii] = np.exp(-n*(T_i[:]-self.T_0)**2/(2*self.sigma**2)) #sign + if (i_n%2 == 0 and i_m == int(i_n/2)): + M[:,ii] = M[:,ii] + elif (i_m < ceiling(i_n/2)): + if (i_n>0): + M[:,ii] *= np.exp(-1j * (n-2*i_m)*W * T_i[:]) + if (i_n<0): + M[:,ii] *= np.exp(1j * (n-2*i_m)*W * T_i[:]) + return M + + def divide_by_factor(self,f,n,m): # note that this must be generalised and moved to divide by field + # + # 1/2 (iEo)^n (-1)^m K (-p*omega; omega, ..., omega, -omega, ..., -omega) => specific for one freq. + # + # then K = n!/(m!(n-m)!) * 2**(1-n) if p /= 0, omega \=0 + # K = 2**(-n) if p = 0, omega \=0 + # K = 1 p = 0, omega =0 + p = int(n - 2*m) + factor = 0.5 + if f > 1e-4: + factor *=2**(-n) + if p != 0: + factor *= 2*factorial(n)/factorial(m)/factorial(n-m) + factor *= (-1)**m*np.power(1.0j*efield['amplitude'],n,dtype=np.cdouble) + return 1.0/factor + + def output_analysis(self,out,to_file=True): + i_map = self.build_map() + for ii in range(len(i_map)): + i_n,i_m = i_map[ii] + if (i_n < 0): continue + i_p = int(i_n - 2 * i_m) + for i_f in range(self.n_runs): + out[ii, i_f, :] *= divide_by_factor(self.efields[i_f],self.freqs[i_f],i_n,i_m) + out[ii,:,:]*=self.get_Unit_of_Measure(i_n) + if (to_file): + output_file = f'o{self.prefix}.YPy-X_order_{i_n}_{i_p}' + header = "E[eV] " + " ".join([f"X{i_n}({i_p}w)/Im({d}) X{i_n}({i_p}w)/Re({d})" for d in ('x','y','z')]) + values = np.column_stack((self.freqs * ha2ev, out[ii, :, 0].imag, out[ii, :, 0].real, + out[ii, :, 1].imag, out[ii, :, 1].real, + out[ii, :, 2].imag, out[ii, :, 2].real)) + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Pulse analysis results") + else: + return (self.freqs, out) - # def reconstruct_signal(self,out,to_file=True): - # Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) - # for i_f in tqdm(range(self.n_runs)): - # for i_d in range(3): - # for i_order in range(self.X_order + 1): - # freq_term = np.exp(-1j * i_order * self.freqs[i_f] * self.time) - # Seff[i_f, i_d, :] += out[i_order, i_f, i_d] * freq_term - # Seff[i_f, i_d, :] += np.conj(out[i_order, i_f, i_d]) * np.conj(freq_term) - # if (to_file): - # for i_f in tqdm(range(self.n_runs)): - # values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) - # output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' - # header="[fs] Px Py Pz" - # if self.l_out_current: - # output_file = f'o{self.prefix}.YamboPy-curr_reconstructed_F{i_f + 1}' - # header="[fs] Jx Jy Jz" - # np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") - # else: - # return values + def reconstruct_signal(self,out,to_file=True): + i_map = self.build_map() + Seff = np.zeros((self.n_runs, 3, len(self.time)), dtype=np.cdouble) + for i_f in tqdm(range(self.n_runs)): + for i_d in range(3): + for ii in range(len(i_map)): + i_n,i_m = i_map[ii] + if (i_n < 0): continue + i_p = int(i_n - 2 * i_m) + freq_term = (np.exp(-1j * i_p * self.freqs[i_f] * self.time))*np.exp(-i_n*(self.time-self.T_0)**2/(2*self.sigma**2)) + Seff[i_f, i_d, :] += out[i_order, i_f, i_d] * freq_term + Seff[i_f, i_d, :] += np.conj(out[i_order, i_f, i_d]) * np.conj(freq_term) + for i_f in tqdm(range(self.n_runs)): + values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) + output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' + header="[fs] Px Py Pz" + if (to_file): + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") + if (not to_file): + return values diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index 65245d4b..db776bec 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -7,7 +7,7 @@ # Modified to allow generalisation # import numpy as np -from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 +from yambopy.units import ha2ev,fs2aut from yambopy.nl.external_efield import Divide_by_the_Field from tqdm import tqdm import sys @@ -101,16 +101,16 @@ def reconstruct_signal(self,out,to_file=True): for i_order in range(self.X_order + 1): freq_term = np.exp(-1j * i_order * self.freqs[i_f] * self.time) Seff[i_f, i_d, :] += out[i_order, i_f, i_d] * freq_term - Seff[i_f, i_d, :] += np.conj(out[i_order, i_f, i_d]) * np.conj(freq_term) - if (to_file): - for i_f in tqdm(range(self.n_runs)): - values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) - output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' - header="[fs] Px Py Pz" - if self.l_out_current: - output_file = f'o{self.prefix}.YamboPy-curr_reconstructed_F{i_f + 1}' - header="[fs] Jx Jy Jz" + Seff[i_f, i_d, :] += np.conj(out[i_order, i_f, i_d]) * np.conj(freq_term) + for i_f in tqdm(range(self.n_runs)): + values = np.column_stack((self.time / fs2aut, Seff[i_f, 0, :].real, Seff[i_f, 1, :].real, Seff[i_f, 2, :].real)) + output_file = f'o{self.prefix}.YamboPy-pol_reconstructed_F{i_f + 1}' + header="[fs] Px Py Pz" + if self.l_out_current: + output_file = f'o{self.prefix}.YamboPy-curr_reconstructed_F{i_f + 1}' + header="[fs] Jx Jy Jz" + if (to_file): np.savetxt(output_file, values, header=header, delimiter=' ', footer="Reconstructed signal") - else: - return values + else: + return values From 93decbcf277afab45369e8b8827ade4b2d76d5b2 Mon Sep 17 00:00:00 2001 From: Myrta Date: Wed, 24 Dec 2025 12:42:00 +0100 Subject: [PATCH 10/23] Changes to be committed: modified: yambopy/nl/Xn_from_signal_diagrams.md Added a brief theory section to the documentation --- yambopy/nl/Xn_from_signal_diagrams.md | 34 +++++++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/yambopy/nl/Xn_from_signal_diagrams.md b/yambopy/nl/Xn_from_signal_diagrams.md index ad679b7e..a99e0f9a 100644 --- a/yambopy/nl/Xn_from_signal_diagrams.md +++ b/yambopy/nl/Xn_from_signal_diagrams.md @@ -1,6 +1,34 @@ -# `Xn_from_signal` – Mermaid Diagrams +# `Xn_from_signal`: documentation (version 1.0) -This document contains Mermaid diagrams derived from `nl_analysis.py` for the abstract class **`Xn_from_signal`** and its key workflows. +#### Myrta Grüning, Claudio Attaccalite, Mike Pointeck, Anna Romani, Mao Yuncheng + +This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” by Butcher and Cotter and the other sources listed in the [Bibliography](## Bibliography). A [minimal background](## 0. Theory ) is given in the next session to help understanding the code and facilitate further development. The rest of the code is dedicated to describe the code structure, key workflows, main functions and to provide an essential guide of the code use. + +## 0. Theory + +The problem solved is algebraic: + +$$ M_{kj} S_j = P_k,$$ + +where $P_k$ is the time-dependent polarization (or current) sampled on $N_t$ times $\{t_k\}$ which is output by the `yambo`/`lumen`code; the resulting $S_j$ is proportional to the susceptibility (conductivity) of nonlinear order $j$. The matrix of coefficients $M_{kj}$, of dimensions $N_t \times N_\text{nl}$ contains the time dependence to the applied electric field. So far three physical situations are implemented: +1. a single monochromatic electric field: $ {\bf E}_0 \sin(\omega_0 t)$ +2. two monochromatic electric fields: $ {\bf E}_0 (\sin(\omega_1 t) + \sin(\omega_2 t)) $ +3. a pulse-shaped electric field: $ {\bf E}(t) \sin(\omega_0 t)$. Here, it is assumed that the shape of the pulse ${\bf E}(t)$ varies slowly with respect to the period $2\pi/\omega_0$. So far, only a Gaussian pulse (${\bf E}(t) = {\bf E}_0 \exp(-(t-t_0)^2/(2\sigma^2))/(\sqrt{2}\sigma)$) has been implemented. + +Four solvers are available: + +1. the standard solver for full, well-determined matrix: calls [`numpy.linalg.solve`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html) +2. the least square solver, when $N_t \gg N_\text{nl}$ : calls [`numpy.linalg.lstsq`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsq) +3. the single value decomposition, using the Moore-Penrose pseudoinverse, when $N_t \gg N_\text{nl}$: calls [`numpy.linalg.pinv`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html#numpy.linalg.pinv) +4. the least square solver with an initial guess, when $N_t \gg N_\text{nl}$ : calls [`scipy.optimize.least_squares`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) + +From $S_j$ the susceptibilities (or conductivities) $\chi^{(n)}(-\omega_\sigma, \omega_1, \dots, \omega_n)$ are obtained using the following expression: + +$$ S_j = C_0 K (-\omega_\sigma, \omega_1, \dots, \omega_n)\chi^{(n)}(-\omega_\sigma, \omega_1, \dots, \omega_n) $$ + +where $K(-\omega_\sigma; \omega_1, \dots, \omega_n)$ is a numerical factor that accounts for the intrinsic permutation symmetry depending on the nonlinear order and frequency arguments of $\chi$. $C_0$ is a further numerical factor depending on the applied electric field. + +Details on the implementation can be found in the sources listed in the [Bibliography](## Bibliography) --- @@ -131,3 +159,5 @@ sequenceDiagram Analyzer-->>Impl: output_analysis(out, to_file) Analyzer-->>Impl: reconstruct_signal(out, to_file) ``` + +## Bibliography From 8d8402d3141030d17d9ee852ecb6295b14da2a9d Mon Sep 17 00:00:00 2001 From: Myrta Date: Fri, 26 Dec 2025 16:39:53 +0100 Subject: [PATCH 11/23] Changes to be committed: deleted: yambopy/nl/Xn_from_signal_diagrams.md new file: yambopy/nl/Xn_from_signal_doc.md First draft of documentation completed. Name of documentation file changed. --- yambopy/nl/Xn_from_signal_diagrams.md | 163 ------------------- yambopy/nl/Xn_from_signal_doc.md | 225 ++++++++++++++++++++++++++ 2 files changed, 225 insertions(+), 163 deletions(-) delete mode 100644 yambopy/nl/Xn_from_signal_diagrams.md create mode 100644 yambopy/nl/Xn_from_signal_doc.md diff --git a/yambopy/nl/Xn_from_signal_diagrams.md b/yambopy/nl/Xn_from_signal_diagrams.md deleted file mode 100644 index a99e0f9a..00000000 --- a/yambopy/nl/Xn_from_signal_diagrams.md +++ /dev/null @@ -1,163 +0,0 @@ -# `Xn_from_signal`: documentation (version 1.0) - -#### Myrta Grüning, Claudio Attaccalite, Mike Pointeck, Anna Romani, Mao Yuncheng - -This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” by Butcher and Cotter and the other sources listed in the [Bibliography](## Bibliography). A [minimal background](## 0. Theory ) is given in the next session to help understanding the code and facilitate further development. The rest of the code is dedicated to describe the code structure, key workflows, main functions and to provide an essential guide of the code use. - -## 0. Theory - -The problem solved is algebraic: - -$$ M_{kj} S_j = P_k,$$ - -where $P_k$ is the time-dependent polarization (or current) sampled on $N_t$ times $\{t_k\}$ which is output by the `yambo`/`lumen`code; the resulting $S_j$ is proportional to the susceptibility (conductivity) of nonlinear order $j$. The matrix of coefficients $M_{kj}$, of dimensions $N_t \times N_\text{nl}$ contains the time dependence to the applied electric field. So far three physical situations are implemented: -1. a single monochromatic electric field: $ {\bf E}_0 \sin(\omega_0 t)$ -2. two monochromatic electric fields: $ {\bf E}_0 (\sin(\omega_1 t) + \sin(\omega_2 t)) $ -3. a pulse-shaped electric field: $ {\bf E}(t) \sin(\omega_0 t)$. Here, it is assumed that the shape of the pulse ${\bf E}(t)$ varies slowly with respect to the period $2\pi/\omega_0$. So far, only a Gaussian pulse (${\bf E}(t) = {\bf E}_0 \exp(-(t-t_0)^2/(2\sigma^2))/(\sqrt{2}\sigma)$) has been implemented. - -Four solvers are available: - -1. the standard solver for full, well-determined matrix: calls [`numpy.linalg.solve`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html) -2. the least square solver, when $N_t \gg N_\text{nl}$ : calls [`numpy.linalg.lstsq`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsq) -3. the single value decomposition, using the Moore-Penrose pseudoinverse, when $N_t \gg N_\text{nl}$: calls [`numpy.linalg.pinv`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html#numpy.linalg.pinv) -4. the least square solver with an initial guess, when $N_t \gg N_\text{nl}$ : calls [`scipy.optimize.least_squares`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) - -From $S_j$ the susceptibilities (or conductivities) $\chi^{(n)}(-\omega_\sigma, \omega_1, \dots, \omega_n)$ are obtained using the following expression: - -$$ S_j = C_0 K (-\omega_\sigma, \omega_1, \dots, \omega_n)\chi^{(n)}(-\omega_\sigma, \omega_1, \dots, \omega_n) $$ - -where $K(-\omega_\sigma; \omega_1, \dots, \omega_n)$ is a numerical factor that accounts for the intrinsic permutation symmetry depending on the nonlinear order and frequency arguments of $\chi$. $C_0$ is a further numerical factor depending on the applied electric field. - -Details on the implementation can be found in the sources listed in the [Bibliography](## Bibliography) - ---- - -## 1) Class structure (`classDiagram`) -```mermaid -classDiagram - class Xn_from_signal { - <> - %% --- Attributes --- - +time : np.ndarray - +T_step : float - +T_deph : float - +efield : dict - +pumps : list - +efields : list - +n_runs : int - +polarization : np.ndarray - +current : np.ndarray - +l_eval_current : bool - +l_out_current : bool - +X_order : int - +solver : str - +samp_mod : str - +nsamp : int - +T_urange : list[float] - +freqs : np.ndarray - +prefix : str - +out_dim : int - +tol : float - - %% --- Constructor --- - +__init__(nldb, X_order=4, T_range=[-1,-1], l_out_current=False, nsamp=-1, samp_mod='', solver='', tol=1e-10) - - %% --- Special --- - +__str__() str - - %% --- Abstract hooks --- - #set_defaults() - #get_sampling(idir, ifrq) - #define_matrix(samp, ifrq) - #update_time_range() - #get_Unit_of_Measure(i_order) - #output_analysis(out, to_file) - #reconstruct_signal(out, to_file) - - %% --- Concrete methods --- - +solve_lin_system(mat, samp, init=None) np.ndarray - +perform_analysis() np.ndarray - +get_Unit_of_Measure(i_order) float - } - - %% Auxiliary (module-level) functions for completeness - class AuxMath { - +IsSquare(m) bool - +IsWellConditioned(m) bool - +residuals_func(x, M, S_i) float - } -``` - ---- - -## 2) Workflow: `perform_analysis` (`flowchart`) -```mermaid -graph TD - A[perform_analysis] --> B[set_defaults] - B --> C[allocate_out] - C --> D{loop_i_f} - D --> E{loop_i_d} - E --> F[get_sampling] - F --> G[define_matrix] - G --> H[solve_lin_system] - H --> I[assign_out] - I --> J{more_i_d} - J -- Yes --> E - J -- No --> K{more_i_f} - K -- Yes --> D - K -- No --> L[return_out] -``` - ---- - -## 3) Workflow: `solve_lin_system` (`flowchart`) -```mermaid -graph TD - A[solve_lin_system] --> B[init_out] - B --> C{solver_full} - C -- Yes --> D{square_and_well_conditioned} - D -- No --> E[set_solver_lstsq] - D -- Yes --> F[linalg_solve] - C -- No --> G{solver_lstsq} - G -- Yes --> H[lstsq] - G -- No --> I{solver_svd} - I -- Yes --> J[pinv] - J --> K[accumulate_inv_times_samp] - I -- No --> L{solver_lstsq_opt} - L -- Yes --> M{has_init} - M -- Yes --> O[use_init] - M -- No --> N[lstsq_for_x0] - N --> P[concat_real_and_imag] - O --> P - P --> Q[least_squares] - Q --> R[compose_complex_out] - L -- No --> S[return_out] - F --> T[return_out] - H --> T - K --> T - R --> T - E --> H -``` - ---- - -## 4) Sequence view (optional) -```mermaid -sequenceDiagram - participant Analyzer as Xn_from_signal - participant Impl as Subclass (implements abstract hooks) - - Analyzer->>Impl: set_defaults() - loop for each run i_f and direction i_d - Analyzer->>Impl: get_sampling(i_d, i_f) - Impl-->>Analyzer: (samp_time, samp_sig) - Analyzer->>Impl: define_matrix(samp_time, i_f) - Impl-->>Analyzer: matrix - Analyzer->>Analyzer: solve_lin_system(matrix, samp_sig) - Analyzer->>Analyzer: out[:, i_f, i_d] = raw[:out_dim] - end - Analyzer-->>Impl: output_analysis(out, to_file) - Analyzer-->>Impl: reconstruct_signal(out, to_file) -``` - -## Bibliography diff --git a/yambopy/nl/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md new file mode 100644 index 00000000..1d5b249b --- /dev/null +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -0,0 +1,225 @@ + + +# `Xn_from_signal`: documentation (version 1.0) + +#### Myrta Grüning, Claudio Attaccalite, Mike Pointeck, Anna Romani, Mao Yuncheng + +This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” by Butcher and Cotter and the other sources listed in the [Bibliography](## Bibliography). The [minimal background](## 0. Theory ) to understand the code and facilitate further development is given in the next session. The rest of the document is dedicated to describe the [code structure](## Code), key workflows, main functions and to provide an [essential guide](## How to use) of the code use. + +--- + +## 0. Theory + +The problem solved is algebraic: + +$$ M_{kj} S_j = P_k,$$ + +where $P_k$ is the time-dependent polarization (or current) sampled on $N_t$ times $\{t_k\}$ which is output by the `yambo`/`lumen`code; the resulting $S_j$ is proportional to the susceptibility (conductivity) of nonlinear order $j$. The matrix of coefficients $M_{kj}$, of dimensions $N_t \times N_\text{nl}$ contains the time dependence to the applied electric field. So far three physical situations are implemented: +1. a single monochromatic electric field: $ {\bf E}_0 \sin(\omega_0 t)$ +2. two monochromatic electric fields: $ {\bf E}_0 (\sin(\omega_1 t) + \sin(\omega_2 t)) $ +3. a pulse-shaped electric field: $ {\bf E}(t) \sin(\omega_0 t)$. Here, it is assumed that the shape of the pulse ${\bf E}(t)$ varies slowly with respect to the period $2\pi/\omega_0$. So far, only a Gaussian pulse, ${\bf E}(t) = {\bf E}_0 \exp(-(t-t_0)^2/(2\sigma^2))/(\sqrt{2}\sigma)$ has been implemented. + +Four solvers are available: + +1. the standard solver for full, well-determined matrix: calls [`numpy.linalg.solve`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html) +2. the least square solver, when $N_t \gg N_\text{nl}$ : calls [`numpy.linalg.lstsq`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsq) +3. the single value decomposition, using the Moore-Penrose pseudoinverse, when $N_t \gg N_\text{nl}$: calls [`numpy.linalg.pinv`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html#numpy.linalg.pinv) +4. the least square solver with an initial guess, when $N_t \gg N_\text{nl}$ : calls [`scipy.optimize.least_squares`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) + +From $S_j$ the susceptibilities (or conductivities) $\chi^{(n)}(-\omega_\sigma, \omega_1, \dots, \omega_n)$ are obtained using the following expression: + +$$ S_j = C_0 K (-\omega_\sigma, \omega_1, \dots, \omega_n)\chi^{(n)}(-\omega_\sigma, \omega_1, \dots, \omega_n) $$ + +where $K(-\omega_\sigma; \omega_1, \dots, \omega_n)$ is a numerical factor that accounts for the intrinsic permutation symmetry depending on the nonlinear order and frequency arguments of $\chi$. $C_0$ is a further numerical factor depending on the applied electric field (field strength, normalisation factor, nonlinear order). + +Details on the implementation can be found in the sources listed in the [Bibliography](## Bibliography) + +--- + +## 1. Code + +### 1.1 Input +It is supposed the `yambo_nl`code (part of the `yambo`/`lumen` suites) was run in the `nonlinear` mode with a loop on a number of frequencies and generated the `ndb.Nonlinear` (with the corresponding fragments). The nonlinear database `ndb.Nonlinear` is then read using the `nldb`class that outputs an object containing all information and data of the run. This object is the input. Further, the user can change in input the defaults of some analysis parameters (see below). + +### 1.2 Structure + +The code consists of the abstract class, `Xn_from_signal`, and three subclasses corresponding to the three physical situations: +1. `Xn_from_sine`: a single monochromatic electric field +2. `Xn_from_freqmix`: two monochromatic electric fields +3. `Xn_from_pulse`: a pulse-shaped electric field + +The main method in the abstract class `Xn_from_signal` is `perform_analysis` defining the sequence of operations to be performed. This is shown in the diagram below. First, defaults are set for each implementation (`set_defaults`). Then a double loop is entered on field frequencies and directions. For each implementation, the time-dependent signal is sampled (`get_sampling`) and the sampled-signal $P_k$ (`samp_sig`) together with the sampling times $\{t_k\}$ (`samp_time`) is returned. Next, for each implementation, the elements of the matrix $M_{kj}$ (`matrix`) is defined (`define_matrix`). Finally, the linear system is solved (`solve_lin_system`, common to all implementation) and the output passed to the `out` array. The latter is the input of `output_analysis` and `reconstruct_signal` that are implemented in each subclass. + +~~~mermaid +sequenceDiagram + participant Analyzer as Xn_from_signal + participant Impl as Subclass (implements abstract hooks) + +Analyzer->>Impl: set_defaults() +loop for each frequency i_f and direction i_d + Analyzer->>Impl: get_sampling(i_d, i_f) + Impl-->>Analyzer: (samp_time, samp_sig) + Analyzer->>Impl: define_matrix(samp_time, i_f) + Impl-->>Analyzer: matrix + Analyzer->>Analyzer: solve_lin_system(matrix, samp_sig) + Analyzer->>Analyzer: out[:, i_f, i_d] = raw[:out_dim] +end +~~~ +### 1.3 Abstract class diagram + +Attributes, constructor, abstract and concrete methods included in `Xn_from_signal`. + + +```mermaid +classDiagram + class Xn_from_signal { + <> + %% --- Attributes --- + +time : np.ndarray + +T_step : float + +T_deph : float + +efield : dict + +pumps : list + +efields : list + +n_runs : int + +polarization : np.ndarray + +current : np.ndarray + +l_eval_current : bool + +l_out_current : bool + +X_order : int + +solver : str + +samp_mod : str + +nsamp : int + +T_urange : list[float] + +freqs : np.ndarray + +prefix : str + +out_dim : int + +tol : float + + %% --- Constructor --- + +__init__(nldb, X_order=4, T_range=[-1,-1], l_out_current=False, nsamp=-1, samp_mod='', solver='', tol=1e-10) + + %% --- Special --- + +__str__() str + + %% --- Abstract hooks --- + #set_defaults() + #get_sampling(idir, ifrq) + #define_matrix(samp, ifrq) + #update_time_range() + #get_Unit_of_Measure(i_order) + #output_analysis(out, to_file) + #reconstruct_signal(out, to_file) + + %% --- Concrete methods --- + +solve_lin_system(mat, samp, init=None) np.ndarray + +perform_analysis() np.ndarray + +get_Unit_of_Measure(i_order) float + } + + %% Auxiliary (module-level) functions for completeness + class AuxMath { + +IsSquare(m) bool + +IsWellConditioned(m) bool + +residuals_func(x, M, S_i) float + } +``` + +### 1.4 Subclasses diagram + +The subclasses inherit the attributes and implement the abstract methods from `Xn_from_signal`. In addition: + +* `Xn_from_freqmix` has the `pump_freq` attribute, which is the frequency of the second electric field, and the `spike_correction` method which performs again the analysis for data points where the simple least square algebraic solution failed, using least square optimization starting from the averaged solution of the neighbouring data points. +* `Xn_from_pulse`has the `T0` attribute, centre of the Gaussian, and the `build_map` and `divide_by_factor`. The former maps the nonlinearity order $n$ and the number of negative frequencies $m$ into a single index. The latter is a modification of `Divide_by_Field` and is in this class until a proper generalisation of `Divide_by_Field` is available. + +```mermaid +classDiagram +note "Attributes and methods of Xn_from_signal in previous diagram" +Xn_from_signal <|-- Xn_from_sine +Xn_from_signal <|-- Xn_from_freqmix +Xn_from_signal <|-- Xn_from_pulse +class Xn_from_pulse{ ++ T0 : float ++ build_map() ++ divide_by_factor() +} +class Xn_from_freqmix{ ++ pump_freq : float ++ spike_correction() +} +class Xn_from_sine +``` +--- +## How to use + +The diagram below illustrates the general use of the code. + +Once the appropriate `ndb.Nonlinear` database has been created with `yambo_nl`, the `YamboNLDB` class is used to read the database and create the object containing all information. + +Depending on the external field used in `yambo_nl`, a different subclass is instantiated: + +* `Xn_from_sine`: a single monochromatic electric field +* `Xn_from_freqmix`: two monochromatic electric fields +* `Xn_from_pulse`: a pulse-shaped electric field + +All subclasses implement the `perform_analysis()` method (setting up and solving the algebraic problem in the [Theory section](## 0. Theory)). The `output_analysis` writes (by default) the susceptibilities (conductivities) on files. For checking the goodness of the analysis, one may output the reconstructed signal (to be compared with the input signal) with `reconstruct_signal`. + + +```mermaid + +graph LR + F[(ndb.Nonlinear)]--> A(Read database with YamboNLDB) -->|NLDB| B1(Xn_from_sine) -->|SIG| C(perform_analysis) -->|OUT| D(output_analysis) --> o.*@{shape: lean-l} + C --> |OUT|E(reconstruct_signal) + E --> o.*@{shape: lean-l} + A -->|NLDB| B2(Xn_from_freqmix) + A -->|NLDB| B3(Xn_from_pulse) + B2 -->|SIG| C + B3 -->|SIG| C + +``` + + + +### Example 1: Monochromatic external field. Harmonic generation from polarization. + +For one monochromatic external field, the `Xn_from_sine` class is instantiated. One can `print` the instance `SIG` to check the value of the class attributes read from `NLDB` and the defaults. Here, the default for `X_order` is overwritten and set to `5`. The output of `perform_analysis` is passed to `OUT`. This is passed as input to `output_analysis` that outputs ` o-DBs.YamboPy-X_probe_order_?`. Second harmonic is in `o-DBs.YamboPy-X_probe_order_2`. + +```python +NLDB=YamboNLDB(calc='DBs') +SIG = Xn_from_sine(NLDB,X_order=5) +print(SIG) +OUT = SIG.perform_analysis() +SIG.output_analysis(OUT) +``` +### Example 2: Monochromatic external field. Shift-current from current. + +For one monochromatic external field, the `Xn_from_sine` class is instantiated. Here, the default for `l_out_current` is overwritten and set to `True`. This means that the current, rather than the polarization is analysed. Note that the `yambo_nl` run must output the current together with polarization. The output of `perform_analysis` is passed to `output_analysis` that outputs ` o-DBs.YamboPy-Sigma_probe_order_?` (shift current is in the `o-DBs.YamboPy-X_probe_order_0` ), and to `reconstruct_signal` that outputs the files `o-DBs.YamboPy-curr_reconstructed_F*` + +```python +NLDB=YamboNLDB(calc="DBs") +SIG = Xn_from_sine(NLDB,l_out_current=True) +OUT = SIG.perform_analysis() +SIG.output_analysis(OUT) +SIG.reconstruct_signal(OUT) +``` +### Example 3: Frequency mixing. Sum/difference of frequencies. + +For two monochromatic external fields, the `Xn_from_freqmix` class is instantiated. Here, the default for the lower limit`T_range` is overwritten and set to `50 fs`. The output of `perform_analysis` is passed to `output_analysis` that in turn outputs `o-DBs.YamboPy-X_probe_order_?_?` each containing the $\chi^{(|n|+|m|)} (-\omega_\sigma; n\omega_1, m\omega_2 )$. The sum/difference of frequencies is given for $n,m=1$ (`o-DBs.YamboPy-X_probe_order_1_1`) and $n=1, m=-1$ (`o-DBs.YamboPy-X_probe_order_1_-1`). + +```python +NLDB=YamboNLDB(calc="DBs") +SIG = Xn_from_freqmix(NLDB,T_range=[50.0*fs2aut,-1.0]) +OUT = SIG.perform_analysis() +SIG.output_analysis(OUT) +``` + +Note: in all snippets one must add `from yambopy import *` + +--- + +## Bibliography + +1. Butcher PN, Cotter D. The constitutive relation. In: The Elements of Nonlinear Optics. Cambridge Studies in Modern Optics. Cambridge University Press; 1990:12-36. +2. Attaccalite C and Grüning M, [Phys. Rev. B 88, 235113 (2013)](https://doi.org/10.1103/PhysRevB.88.235113) +3. Pionteck MN, Grüning M, Sanna S, Attaccalite C, [SciPost Phys. 19, 129 (2025)](https://10.21468/SciPostPhys.19.5.129) +4. Romani A, Grüning M, 'Notes on nonlinear analysis from Gaussian pulses' (unpublished) \ No newline at end of file From 8d7275dab41fc1fcb93cc3d96af6467a202e8180 Mon Sep 17 00:00:00 2001 From: Myrta Date: Sat, 27 Dec 2025 12:56:59 +0100 Subject: [PATCH 12/23] Changes to be committed: modified: yambopy/nl/Xn_from_signal_doc.md Minor edits --- yambopy/nl/Xn_from_signal_doc.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yambopy/nl/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md index 1d5b249b..fa56cf61 100644 --- a/yambopy/nl/Xn_from_signal_doc.md +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -15,9 +15,9 @@ The problem solved is algebraic: $$ M_{kj} S_j = P_k,$$ where $P_k$ is the time-dependent polarization (or current) sampled on $N_t$ times $\{t_k\}$ which is output by the `yambo`/`lumen`code; the resulting $S_j$ is proportional to the susceptibility (conductivity) of nonlinear order $j$. The matrix of coefficients $M_{kj}$, of dimensions $N_t \times N_\text{nl}$ contains the time dependence to the applied electric field. So far three physical situations are implemented: -1. a single monochromatic electric field: $ {\bf E}_0 \sin(\omega_0 t)$ -2. two monochromatic electric fields: $ {\bf E}_0 (\sin(\omega_1 t) + \sin(\omega_2 t)) $ -3. a pulse-shaped electric field: $ {\bf E}(t) \sin(\omega_0 t)$. Here, it is assumed that the shape of the pulse ${\bf E}(t)$ varies slowly with respect to the period $2\pi/\omega_0$. So far, only a Gaussian pulse, ${\bf E}(t) = {\bf E}_0 \exp(-(t-t_0)^2/(2\sigma^2))/(\sqrt{2}\sigma)$ has been implemented. +1. a single monochromatic electric field: ${\bf E}_0 \sin(\omega_0 t)$ +2. two monochromatic electric fields: ${\bf E}_0 (\sin(\omega_1 t) + \sin(\omega_2 t))$ +3. a pulse-shaped electric field: ${\bf E}(t) \sin(\omega_0 t)$. Here, it is assumed that the shape of the pulse ${\bf E}(t)$ varies slowly with respect to the period $2\pi/\omega_0$. So far, only a Gaussian pulse, ${\bf E}(t) = {\bf E}_0 \exp(-(t-t_0)^2/(2\sigma^2))/(\sqrt{2}\sigma)$ has been implemented. Four solvers are available: @@ -222,4 +222,4 @@ Note: in all snippets one must add `from yambopy import *` 1. Butcher PN, Cotter D. The constitutive relation. In: The Elements of Nonlinear Optics. Cambridge Studies in Modern Optics. Cambridge University Press; 1990:12-36. 2. Attaccalite C and Grüning M, [Phys. Rev. B 88, 235113 (2013)](https://doi.org/10.1103/PhysRevB.88.235113) 3. Pionteck MN, Grüning M, Sanna S, Attaccalite C, [SciPost Phys. 19, 129 (2025)](https://10.21468/SciPostPhys.19.5.129) -4. Romani A, Grüning M, 'Notes on nonlinear analysis from Gaussian pulses' (unpublished) \ No newline at end of file +4. Romani A, Grüning M, 'Notes on nonlinear analysis from Gaussian pulses' (unpublished) From f40dcee812e2533d8ea7e08172f9487cfbd91841 Mon Sep 17 00:00:00 2001 From: Myrta Date: Sat, 27 Dec 2025 13:19:20 +0100 Subject: [PATCH 13/23] modified: yambopy/nl/Xn_from_signal_doc.md more minor edits --- yambopy/nl/Xn_from_signal_doc.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/yambopy/nl/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md index fa56cf61..ab8108cf 100644 --- a/yambopy/nl/Xn_from_signal_doc.md +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -4,11 +4,12 @@ #### Myrta Grüning, Claudio Attaccalite, Mike Pointeck, Anna Romani, Mao Yuncheng -This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” by Butcher and Cotter and the other sources listed in the [Bibliography](## Bibliography). The [minimal background](## 0. Theory ) to understand the code and facilitate further development is given in the next session. The rest of the document is dedicated to describe the [code structure](## Code), key workflows, main functions and to provide an [essential guide](## How to use) of the code use. +This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” and the other sources listed in the bibliography. The minimal background to understand the code and facilitate further development is given in the next session. The rest of the document is dedicated to describe the code structure, key workflows, main functions and to provide an essential guide of the code use. --- -## 0. Theory +## 0. Minimal theoretical compendium + The problem solved is algebraic: @@ -32,7 +33,7 @@ $$ S_j = C_0 K (-\omega_\sigma, \omega_1, \dots, \omega_n)\chi^{(n)}(-\omega_\si where $K(-\omega_\sigma; \omega_1, \dots, \omega_n)$ is a numerical factor that accounts for the intrinsic permutation symmetry depending on the nonlinear order and frequency arguments of $\chi$. $C_0$ is a further numerical factor depending on the applied electric field (field strength, normalisation factor, nonlinear order). -Details on the implementation can be found in the sources listed in the [Bibliography](## Bibliography) +Details on the implementation can be found in the sources listed in the bibliography. --- @@ -219,7 +220,7 @@ Note: in all snippets one must add `from yambopy import *` ## Bibliography -1. Butcher PN, Cotter D. The constitutive relation. In: The Elements of Nonlinear Optics. Cambridge Studies in Modern Optics. Cambridge University Press; 1990:12-36. +1. Butcher PN, Cotter D. The constitutive relation. In: The Elements of Nonlinear Optics. Cambridge Studies in Modern Optics. Cambridge University Press; 1990:12-36.[doi.org/10.1017/CBO9781139167994](https://doi.org/10.1017/CBO9781139167994) 2. Attaccalite C and Grüning M, [Phys. Rev. B 88, 235113 (2013)](https://doi.org/10.1103/PhysRevB.88.235113) 3. Pionteck MN, Grüning M, Sanna S, Attaccalite C, [SciPost Phys. 19, 129 (2025)](https://10.21468/SciPostPhys.19.5.129) 4. Romani A, Grüning M, 'Notes on nonlinear analysis from Gaussian pulses' (unpublished) From dd1afe34e49d4de49be5668d9bedbe9b589f4b10 Mon Sep 17 00:00:00 2001 From: Myrta Date: Sun, 28 Dec 2025 15:39:43 +0100 Subject: [PATCH 14/23] modified: yambopy/nl/Xn_from_signal_doc.md minor edit: fixed crossed references --- yambopy/nl/Xn_from_signal_doc.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/yambopy/nl/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md index ab8108cf..c8ef3cb5 100644 --- a/yambopy/nl/Xn_from_signal_doc.md +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -4,7 +4,7 @@ #### Myrta Grüning, Claudio Attaccalite, Mike Pointeck, Anna Romani, Mao Yuncheng -This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” and the other sources listed in the bibliography. The minimal background to understand the code and facilitate further development is given in the next session. The rest of the document is dedicated to describe the code structure, key workflows, main functions and to provide an essential guide of the code use. +This document describes the `Xn_from_signal` abstract python class and a set of derived classes, part of the `YamboPy` code, for extracting nonlinear susceptibilities and conductivities from the macroscopic time-dependent polarization $P$ and current $J$. An extended theoretical background can be found in Nonlinear Optics textbooks, see e.g. Sec. 2 of “The Elements of Nonlinear Optics” and the other sources listed in the [bibliography](# Bibliography). The [minimal background](# 0. Minimal theoretical compendium) to understand the code and facilitate further development is given in the next session. The rest of the document is dedicated to describe the [code structure, key workflows, main functions](# 1. Code) and to provide an essential guide of the [code use](# 2. How to use). --- @@ -151,7 +151,7 @@ class Xn_from_freqmix{ class Xn_from_sine ``` --- -## How to use +## 2. How to use The diagram below illustrates the general use of the code. @@ -224,3 +224,5 @@ Note: in all snippets one must add `from yambopy import *` 2. Attaccalite C and Grüning M, [Phys. Rev. B 88, 235113 (2013)](https://doi.org/10.1103/PhysRevB.88.235113) 3. Pionteck MN, Grüning M, Sanna S, Attaccalite C, [SciPost Phys. 19, 129 (2025)](https://10.21468/SciPostPhys.19.5.129) 4. Romani A, Grüning M, 'Notes on nonlinear analysis from Gaussian pulses' (unpublished) + +[# Bibliography]: From 6d2ac5714ddf52c5d30a4b7dc87b9fb7278b565c Mon Sep 17 00:00:00 2001 From: Myrta Date: Sun, 28 Dec 2025 15:48:06 +0100 Subject: [PATCH 15/23] Update Xn_from_signal_doc.md minor edits --- yambopy/nl/Xn_from_signal_doc.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/yambopy/nl/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md index c8ef3cb5..fdd734d2 100644 --- a/yambopy/nl/Xn_from_signal_doc.md +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -224,5 +224,3 @@ Note: in all snippets one must add `from yambopy import *` 2. Attaccalite C and Grüning M, [Phys. Rev. B 88, 235113 (2013)](https://doi.org/10.1103/PhysRevB.88.235113) 3. Pionteck MN, Grüning M, Sanna S, Attaccalite C, [SciPost Phys. 19, 129 (2025)](https://10.21468/SciPostPhys.19.5.129) 4. Romani A, Grüning M, 'Notes on nonlinear analysis from Gaussian pulses' (unpublished) - -[# Bibliography]: From e1ae6976837ba302d1dcb2f8f5bac8365efd7272 Mon Sep 17 00:00:00 2001 From: myrta75 Date: Sat, 3 Jan 2026 10:25:32 +0000 Subject: [PATCH 16/23] modified: yambopy/__init__.py modified: yambopy/dbs/nldb.py modified: yambopy/nl/external_efield.py modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/nl_analysis.py modified: yambopy/nl/pulse_analysis.py Corrected small bugs and added minor features --- yambopy/__init__.py | 1 + yambopy/dbs/nldb.py | 24 +++++++++++++++++---- yambopy/nl/external_efield.py | 3 ++- yambopy/nl/freqmix_analysis.py | 2 +- yambopy/nl/nl_analysis.py | 31 ++++++++++++++++++++++----- yambopy/nl/pulse_analysis.py | 39 ++++++++++++++++++++++------------ 6 files changed, 76 insertions(+), 24 deletions(-) diff --git a/yambopy/__init__.py b/yambopy/__init__.py index be033b6b..ca002da8 100644 --- a/yambopy/__init__.py +++ b/yambopy/__init__.py @@ -122,6 +122,7 @@ from yambopy.nl.nl_analysis import * from yambopy.nl.sin_analysis import * from yambopy.nl.freqmix_analysis import * +from yambopy.nl.pulse_analysis import * #doublegrid files from yambopy.double_grid.dg_convergence import * diff --git a/yambopy/dbs/nldb.py b/yambopy/dbs/nldb.py index e3f51c3e..802c80ba 100644 --- a/yambopy/dbs/nldb.py +++ b/yambopy/dbs/nldb.py @@ -43,11 +43,27 @@ def read_Efield(self,database,RT_step,n): efield["name"] =database.variables['Field_Name_'+str(n)][...].tobytes().decode().strip() efield["versor"] =database.variables['Field_Versor_'+str(n)][:].astype(np.double) efield["intensity"] =database.variables['Field_Intensity_'+str(n)][0].astype(np.double) - efield["damping"] =database.variables['Field_Damping_'+str(n)][0].astype(np.double) - efield["freq_range"] =database.variables['Field_Freq_range_'+str(n)][:].astype(np.double) - efield["freq_steps"] =database.variables['Field_Freq_steps_'+str(n)][:].astype(np.double) - efield["freq_step"] =database.variables['Field_Freq_step_'+str(n)][0].astype(np.double) + try: + efield["damping"] =database.variables['Field_FWHM_'+str(n)][0].astype(np.double) + except: + efield["damping"] =database.variables['Field_Damping_'+str(n)][0].astype(np.double) + try: + efield["freq_range"] =database.variables['Field_Freq_range_'+str(n)][:].astype(np.double) + except: + efield["freq_range"] =database.variables['Field_Freq_'+str(n)][:].astype(np.double) + try: + efield["freq_steps"] =database.variables['Field_Freq_steps_'+str(n)][:].astype(np.double) + except: + efield["freq_steps"] =1.0 + try: + efield["freq_step"] =database.variables['Field_Freq_step_'+str(n)][0].astype(np.double) + except: + efield["freq_step"] =0.0 efield["initial_time"] =database.variables['Field_Initial_time_'+str(n)][0].astype(np.double) + try: + efield["peak"] =database.variables['Field_peak_'+str(n)][0].astype(np.double) + except: + efield["peak"] =10.0 # # set t_initial according to Yambo # diff --git a/yambopy/nl/external_efield.py b/yambopy/nl/external_efield.py index 9b3c7ba3..7a6ca76b 100644 --- a/yambopy/nl/external_efield.py +++ b/yambopy/nl/external_efield.py @@ -42,5 +42,6 @@ def Divide_by_the_Field(efield,order): def Gaussian_centre(efield): ratio=np.pi/efield['freq_range'][0] sigma=efield["damping"]/(2.0*(2.0*np.log(2.0))**0.5) - return ratio * float(round(1.0 /ratio * sigma *efield["field_peak"])) + peak_ = float(1.0 /ratio * sigma *efield["peak"]) + return ratio * float(round(peak_)),sigma diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index 4baa4a55..fbc6ebef 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -10,7 +10,7 @@ from yambopy.units import ha2ev,fs2aut, Junit, EFunit,SVCMm12VMm1,AU2VMm1 from yambopy.nl.external_efield import Divide_by_the_Field from tqdm import tqdm -from scipy.andimage import uniform_filter1d +from scipy.ndimage import uniform_filter1d import sys import os import itertools diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index c2f94168..54fe7011 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -66,7 +66,7 @@ def __str__(self): if self.nsamp > 0: s+="Sampling points : "+str(self.nsamp)+"\n" if self.T_urange!=[-1, -1]: - s+="User time range : "+str(self.T_urange)+"\n" + s+="User time range : "+str(self.T_urange)+" [au] \n" s+="Frequency range: ["+str(self.freqs[0])+","+str(self.freqs[-1])+"] [au] \n" if self.l_out_current: s+="Output is set to current." @@ -97,8 +97,8 @@ def get_Unit_of_Measure(self,i_order): def solve_lin_system(self,mat,samp,init=None): mat_dim = int(mat.shape[1]) out=np.zeros(mat_dim,dtype=np.cdouble) - if self.solver=="full" and ((not IsSquare(mat)) or (not IsWellDefined(mat))): - print('WARNING: solver changed to least square') + if self.solver=="full" and ((not self.IsSquare(mat)) or (not self.IsWellConditioned(mat))): + print(f'WARNING: solver changed to least square since square:{self.IsSquare(mat)} and well-conditioned:{self.IsWellConditioned(mat)}') self.solver = "lstsq" if self.solver=="full": out = np.linalg.solve(mat,samp) @@ -124,8 +124,11 @@ def perform_analysis(self): for i_f in tqdm(range(self.n_runs)): for i_d in range(3): samp_time, samp_sig= self.get_sampling(i_d,i_f) + print(i_f,i_d,samp_time, samp_sig) matrix = self.define_matrix(samp_time,i_f) + print(i_f,i_d,matrix) raw = self.solve_lin_system(matrix,samp_sig) + print(i_f,i_d,raw) out[:, i_f, i_d] = raw[:self.out_dim] return out @@ -146,11 +149,29 @@ def output_analysis(self,out,to_file): @abstractmethod def reconstruct_signal(self,out,to_file): pass + + def append_runinfo(self,T): + s = "# Field details:" + s+= "# Type field "+str(self.efield["name"])+"\n" + s+= "# Field intensity "+str(self.efield["intensity"])+"\n" + s+= "# Field versor "+str(self.efield["versor"])+"\n" + if self.efield["name"]=='QSSIN': + s+= "Centre Gaussian "+str(self.T0/fs2aut)+"[fs] \n" + s+= "Gaussian sigma "+str(self.sigma/fs2aut)+"[fs] \n" + s+= "# Analysis details:" + s+="# Max harmonic order : "+str(self.X_order)+"\n" + s+="# Sampling mode : "+str(self.samp_mod) +"\n" + s+="# Solver : "+str(self.solver)+"\n" + s+="# Sampling points : "+str(self.nsamp)+"\n" + s+="# Start sampling time : "+str(T/fs2aut)+" [fs] \n" + return s + + ### some maths auxiliary functions - def IsSquare(m): + def IsSquare(self,m): return m.shape[0] == m.shape[1] - def IsWellConditioned(m): # with this I am trying to avoid inverting a matrix + def IsWellConditioned(self,m): # with this I am trying to avoid inverting a matrix return np.linalg.cond(m) < 1/sys.float_info.epsilon def residuals_func(x,M,S_i): diff --git a/yambopy/nl/pulse_analysis.py b/yambopy/nl/pulse_analysis.py index 9a7ed41b..db4c35d5 100644 --- a/yambopy/nl/pulse_analysis.py +++ b/yambopy/nl/pulse_analysis.py @@ -29,10 +29,13 @@ def set_defaults(self): raise ValueError("This analysis is for one monochromatic field only.") if self.solver == '': self.solver = 'full' - dim_ = int(1+self.X_order+(self.X_order/2)**2) - D = (dim_,int(dim_ -0.25)) + dim_ = 1+self.X_order+(self.X_order/2)**2 + D = (int(dim_),int(dim_ -0.25)) self.out_dim = D[self.X_order%2] - self.T0 = Gaussian_centre(self.efield) + self.T0 = Gaussian_centre(self.efield)[0] + self.sigma = Gaussian_centre(self.efield)[1] + if self.samp_mod== "": + self.samp_mod="linear" return def build_map(self): @@ -41,10 +44,10 @@ def build_map(self): I = np.zeros((M,2),dtype=int) j = 0 I[j,:] = [0,0] - for n in range(2,N+1,2): + for n in range(2,self.X_order+1,2): j += 1 I[j,:] = [n,floor(n/2)] - for n in range(1,N+1): + for n in range(1,self.X_order+1): for m in range(0,ceil(n/2)): j += 1 I[j,:] = [n,m] @@ -52,9 +55,10 @@ def build_map(self): return I def get_sampling(self,idir,ifrq): - M = 1+2*self.X_order + int(self.X_order*(self.X_order-1)/2) + M = int(1+2*self.X_order + int(self.X_order*(self.X_order-1)/2)) if self.nsamp == -1 or self.nsamp < M: self.nsamp = M + print(f'Sampling points set to {self.nsamp}') T_period = 2.0 * np.pi / self.freqs[ifrq] T_range, out_of_bounds = self.update_time_range(T_period) if (out_of_bounds): @@ -87,10 +91,10 @@ def define_matrix(self,T_i,ifrq): for ii in range(1,len(i_map)): i_n,i_m = i_map[ii] n = abs(i_n) - M[:,ii] = np.exp(-n*(T_i[:]-self.T_0)**2/(2*self.sigma**2)) #sign - if (i_n%2 == 0 and i_m == int(i_n/2)): + M[:,ii] = np.exp(-n*(T_i[:]-self.T0)**2/(2*self.sigma**2)) #sign + if (n%2 == 0 and i_m == int(n/2)): M[:,ii] = M[:,ii] - elif (i_m < ceiling(i_n/2)): + elif (i_m < ceil(n/2)): if (i_n>0): M[:,ii] *= np.exp(-1j * (n-2*i_m)*W * T_i[:]) if (i_n<0): @@ -106,11 +110,11 @@ def divide_by_factor(self,f,n,m): # note that this must be generalised and moved # K = 1 p = 0, omega =0 p = int(n - 2*m) factor = 0.5 - if f > 1e-4: - factor *=2**(-n) + if self.freqs[f] > 1e-4: + factor *=np.power(2,-n,dtype=np.cdouble) if p != 0: factor *= 2*factorial(n)/factorial(m)/factorial(n-m) - factor *= (-1)**m*np.power(1.0j*efield['amplitude'],n,dtype=np.cdouble) + factor *= (-1)**m*np.power(1.0j*self.efield['amplitude'],n,dtype=np.cdouble) return 1.0/factor def output_analysis(self,out,to_file=True): @@ -119,9 +123,14 @@ def output_analysis(self,out,to_file=True): i_n,i_m = i_map[ii] if (i_n < 0): continue i_p = int(i_n - 2 * i_m) + T = 10000.0 for i_f in range(self.n_runs): - out[ii, i_f, :] *= divide_by_factor(self.efields[i_f],self.freqs[i_f],i_n,i_m) + out[ii, i_f, :] *= self.divide_by_factor(i_f,i_n,i_m) + T_period = 2.0 * np.pi / self.freqs[i_f] + Trange, _ = self.update_time_range(T_period) + T = min(T,Trange[0]) out[ii,:,:]*=self.get_Unit_of_Measure(i_n) + run_info = self.append_runinfo(T) if (to_file): output_file = f'o{self.prefix}.YPy-X_order_{i_n}_{i_p}' header = "E[eV] " + " ".join([f"X{i_n}({i_p}w)/Im({d}) X{i_n}({i_p}w)/Re({d})" for d in ('x','y','z')]) @@ -129,7 +138,11 @@ def output_analysis(self,out,to_file=True): out[ii, :, 1].imag, out[ii, :, 1].real, out[ii, :, 2].imag, out[ii, :, 2].real)) np.savetxt(output_file, values, header=header, delimiter=' ', footer="Pulse analysis results") + outf = open(output_file, "a") + outf.writelines(run_info) + outf.close() else: + print(run_info) return (self.freqs, out) def reconstruct_signal(self,out,to_file=True): From b36ef2adc3e6be2321a87bf7f7d156ed9fd6424c Mon Sep 17 00:00:00 2001 From: myrta75 Date: Mon, 5 Jan 2026 09:40:37 +0000 Subject: [PATCH 17/23] modified: yambopy/nl/pulse_analysis.py Changed default for starting sampling time to T0 - sigma. --- yambopy/nl/pulse_analysis.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yambopy/nl/pulse_analysis.py b/yambopy/nl/pulse_analysis.py index db4c35d5..e6d1f5fd 100644 --- a/yambopy/nl/pulse_analysis.py +++ b/yambopy/nl/pulse_analysis.py @@ -74,10 +74,10 @@ def update_time_range(self,T_period): T_range = self.T_urange out_of_bounds = False if T_range[0] <= 0.0: - T_range[0] = self.T0 - T_period/2 + T_range[0] = self.T0 - self.sigma T_range[1] = T_range[0] + T_period if T_range[1] > self.time[-1]: - T_range[1] = slef.time[-1] + T_range[1] = self.time[-1] T_range[0] = T_range[1] - T_period out_of_bound = True return T_range, out_of_bounds From ac9d3d4ed5f895ae6ba7027875dabfc96a68a5e5 Mon Sep 17 00:00:00 2001 From: Myrta Date: Mon, 5 Jan 2026 18:34:20 +0100 Subject: [PATCH 18/23] modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/sin_analysis.py added run info at the end of output files --- yambopy/nl/freqmix_analysis.py | 8 +++++++- yambopy/nl/sin_analysis.py | 13 +++++++++++-- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index fbc6ebef..dfede9df 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -88,6 +88,8 @@ def define_matrix(self,T_i,ifrq): def output_analysis(self,out,to_file=True): # NX,MX = self.X_order[:] + T, _ = update_time_range() + run_info = self.append_runinfo(T) # for i_f in range(self.n_runs): for i_c,(i_n,i_m) in enumerate(itertools.product(range(-NX, NX+1),range(-MX, MX+1))): @@ -105,8 +107,12 @@ def output_analysis(self,out,to_file=True): values = np.column_stack((self.freqs * ha2ev, out[i_c, :, 0].imag, out[i_c, :, 0].real, out[i_c, :, 1].imag, out[i_c, :, 1].real, out[i_c, :, 2].imag, out[i_c, :, 2].real)) - np.savetxt(output_file, values, header=header, delimiter=' ', footer=f"Frequency mixing analysis with pump frequency {self.pump_freq*ha2ev} eV") + np.savetxt(output_file, values, header=header, delimiter=' ', footer=f"Frequency mixing analysis with pump frequency {self.pump_freq*ha2ev} eV") + outf = open(output_file, "a") + outf.writelines(run_info) + outf.close() else: + print(run_info) return (self.freqs, out) def reconstruct_signal(self,out,to_file=True): diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index db776bec..742c057a 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -78,9 +78,14 @@ def define_matrix(self,T_i,ifrq): def output_analysis(self,out,to_file=True): for i_order in range(self.X_order + 1): + T = 10000.0 for i_f in range(self.n_runs): out[i_order, i_f, :] *= Divide_by_the_Field(self.efields[i_f], i_order) - out[i_order,:,:]*=self.get_Unit_of_Measure(i_order) + T_period = 2.0 * np.pi / self.freqs[i_f] + Trange, _ = self.update_time_range(T_period) + T = min(T,Trange[0]) + out[i_order,:,:]*=self.get_Unit_of_Measure(i_order) + run_info = self.append_runinfo(T) if (to_file): output_file = f'o{self.prefix}.YamboPy-X_probe_order_{i_order}' header = "E[eV] " + " ".join([f"X{i_order}/Im({d}) X{i_order}/Re({d})" for d in ('x','y','z')]) @@ -90,8 +95,12 @@ def output_analysis(self,out,to_file=True): values = np.column_stack((self.freqs * ha2ev, out[i_order, :, 0].imag, out[i_order, :, 0].real, out[i_order, :, 1].imag, out[i_order, :, 1].real, out[i_order, :, 2].imag, out[i_order, :, 2].real)) - np.savetxt(output_file, values, header=header, delimiter=' ', footer="Harmonic analysis results") + np.savetxt(output_file, values, header=header, delimiter=' ', footer="Harmonic analysis results") + outf = open(output_file, "a") + outf.writelines(run_info) + outf.close() else: + print(run_info) return (self.freqs, out) def reconstruct_signal(self,out,to_file=True): From d764e47133abe177316d28e03bbbb2b3437e5efb Mon Sep 17 00:00:00 2001 From: myrta75 Date: Thu, 8 Jan 2026 16:22:43 +0000 Subject: [PATCH 19/23] modified: yambopy/nl/nl_analysis.py Added debug mode to hide comments --- yambopy/nl/nl_analysis.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index 54fe7011..e0b3ea60 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -1,4 +1,4 @@ -# Copyright (c) 2023-2025, Claudio Attaccalite, +# Copyright (c) 2023-2026, Claudio Attaccalite, # Myrta Gruening # All rights reserved. # @@ -19,7 +19,7 @@ # Template class # class Xn_from_signal(ABC): - def __init__(self,nldb,X_order=4,T_range=[-1, -1],l_out_current=False,nsamp=-1,samp_mod='',solver='',tol=1e-10): # note I need to add user opportunity to set time-range + def __init__(self,nldb,X_order=4,T_range=[-1, -1],l_out_current=False,nsamp=-1,samp_mod='',solver='',tol=1e-10,debug_mode=False): self.time = nldb.IO_TIME_points # Time series self.T_step =self.time[1] - self.time[0] # Time step of the simulation self.T_deph =12.0/nldb.NL_damping @@ -46,6 +46,7 @@ def __init__(self,nldb,X_order=4,T_range=[-1, -1],l_out_current=False,nsamp=-1,s self.prefix = f'-{nldb.calc}' if nldb.calc != 'SAVE' else '' self.out_dim = 0 self.tol = tol + self.debug_mode = debug_mode super().__init__() def __str__(self): @@ -124,12 +125,17 @@ def perform_analysis(self): for i_f in tqdm(range(self.n_runs)): for i_d in range(3): samp_time, samp_sig= self.get_sampling(i_d,i_f) - print(i_f,i_d,samp_time, samp_sig) matrix = self.define_matrix(samp_time,i_f) - print(i_f,i_d,matrix) raw = self.solve_lin_system(matrix,samp_sig) - print(i_f,i_d,raw) out[:, i_f, i_d] = raw[:self.out_dim] + if self.debug_mode: + print(f"Freq #{i_f}, direction: {i_d}:") + print("***Sampling:") + print(samp_time, samp_sig) + print("***Matrix:") + print(matrix) + print("***Solution:") + print(raw) return out def get_Unit_of_Measure(self,i_order): # not sure if this is a general or specific method - let it here for the moment From aebf75ca78064287f000029a4422e3bc795bc1a2 Mon Sep 17 00:00:00 2001 From: myrta75 Date: Thu, 8 Jan 2026 16:39:18 +0000 Subject: [PATCH 20/23] modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/sin_analysis.py modified: yambopy/nl/sum_frequencies.py Minor changes, corrections and additions --- yambopy/nl/freqmix_analysis.py | 2 +- yambopy/nl/sin_analysis.py | 2 ++ yambopy/nl/sum_frequencies.py | 1 - 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index dfede9df..72f5a0e7 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -88,7 +88,7 @@ def define_matrix(self,T_i,ifrq): def output_analysis(self,out,to_file=True): # NX,MX = self.X_order[:] - T, _ = update_time_range() + T, _ = self.update_time_range() run_info = self.append_runinfo(T) # for i_f in range(self.n_runs): diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index 742c057a..11c18b03 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -29,6 +29,8 @@ def set_defaults(self): if self.solver == '': self.solver = 'full' self.out_dim = self.X_order + 1 + if self.samp_mod== "": + self.samp_mod="linear" return def get_sampling(self,idir,ifrq): diff --git a/yambopy/nl/sum_frequencies.py b/yambopy/nl/sum_frequencies.py index 563b67ac..3d934d1b 100644 --- a/yambopy/nl/sum_frequencies.py +++ b/yambopy/nl/sum_frequencies.py @@ -8,7 +8,6 @@ import math from yambopy.units import ha2ev,fs2aut, SVCMm12VMm1,AU2VMm1 from yambopy.nl.external_efield import Divide_by_the_Field -from yambopy.nl.harmonic_analysis import update_T_range from scipy.optimize import least_squares from tqdm import tqdm from scipy.ndimage import uniform_filter1d From 54099eac0e8d6a88c94e4a5ecc7ac28fe37f8469 Mon Sep 17 00:00:00 2001 From: Myrta Date: Thu, 15 Jan 2026 11:05:19 +0000 Subject: [PATCH 21/23] Changes to be committed: modified: yambopy/nl/freqmix_analysis.py modified: yambopy/nl/nl_analysis.py modified: yambopy/nl/pulse_analysis.py modified: yambopy/nl/sin_analysis.py Moved get_sampling as common method in nl_analysis.py and created the abstract method set_sampling to get the specific range/sampling number --- yambopy/nl/freqmix_analysis.py | 17 ++--------- yambopy/nl/nl_analysis.py | 52 ++++++++++++++++++++++++++++++++-- yambopy/nl/pulse_analysis.py | 8 ++---- yambopy/nl/sin_analysis.py | 11 ++----- 4 files changed, 55 insertions(+), 33 deletions(-) diff --git a/yambopy/nl/freqmix_analysis.py b/yambopy/nl/freqmix_analysis.py index 72f5a0e7..c994c1d7 100644 --- a/yambopy/nl/freqmix_analysis.py +++ b/yambopy/nl/freqmix_analysis.py @@ -54,26 +54,13 @@ def update_time_range(self): T_range[1]=self.time[-1] return T_range # - def get_sampling(self,idir,ifrq): + def set_sampling(self,ifrq): samp_order = self.out_dim if self.nsamp == -1 or self.nsamp < samp_order: self.nsamp = int(samp_order**2) # T_range = self.update_time_range() - i_t_start = int(np.round( T_range[0] / self.T_step)) - i_deltaT = int(np.round((T_range[1]-T_range[0])/self.T_step)/self.nsamp) - if self.samp_mod=='linear': - T_i = (i_t_start + i_deltaT * np.arange(self.nsamp))*self.T_step - P_i = self.polarization[ifrq][idir,i_t_start + i_deltaT * np.arange(self.nsamp)] - elif self.samp_mod=='log': - T_i = np.geomspace(i_t_start * self.T_step, T_range[1], self.nsamp, endpoint=False) - i_t = np.array(np.round(T_i/self.T_step),dtype=int) - T_i = i_t*self.T_step - P_i = np.array([self.polarization[ifrq][idir,i] for i in i_t]) - elif self.samp_mod=='random': - T_i = np.random.uniform(i_t_start * self.T_step, T_range[1], self.nsamp) - P_i = np.array([self.polarization[ifrq][idir,int(np.round(t / self.T_step))] for t in T_i]) - return T_i,P_i + return T_range def define_matrix(self,T_i,ifrq): NX,MX = self.X_order[:] diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index e0b3ea60..81f0b8ef 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -80,8 +80,12 @@ def set_defaults(self): pass @abstractmethod - def get_sampling(self,idir,ifrq): - pass + def set_sampling(self,idir,ifrq): + pass + + @abstractmethod + def get_sampling(self): + pass @abstractmethod def define_matrix(self,samp,ifrq): @@ -94,6 +98,25 @@ def update_time_range(self): @abstractmethod def get_Unit_of_Measure(self,i_order): pass + + def get_sampling(self,T_range,idir,ifrq): # subtract self.efield["Initial_time"] from T_i CHECK!!!! + # + i_t_start = int(np.round( T_range[0] / self.T_step)) + i_deltaT = int(np.round((T_range[1]-T_range[0])/self.T_step)/self.nsamp) + if self.samp_mod=='linear': + i_t = i_t_start + i_deltaT * np.arange(self.nsamp) + elif self.samp_mod=='log': + T_i = np.geomspace(i_t_start * self.T_step, T_range[1], self.nsamp, endpoint=False) + i_t = np.array(np.round(T_i/self.T_step),dtype=int) + elif self.samp_mod=='random': + i_t = np.random.uniform(i_t_start,int(np.round(t / self.T_step)), self.nsamp) + T_i = i_t*self.T_step - self.efield["initial_time"] + if self.l_out_current: + S_i = np.array([self.current[ifrq][idir,i] for i in i_t]) + else: + S_i = np.array([self.polarization[ifrq][idir,i] for i in i_t]) + return T_i,S_i + def solve_lin_system(self,mat,samp,init=None): mat_dim = int(mat.shape[1]) @@ -123,8 +146,9 @@ def perform_analysis(self): _ = self.set_defaults() out = np.zeros((self.out_dim, self.n_runs, 3), dtype=np.cdouble) for i_f in tqdm(range(self.n_runs)): + T_r = self.set_sampling(i_f) for i_d in range(3): - samp_time, samp_sig= self.get_sampling(i_d,i_f) + samp_time, samp_sig= self.get_sampling(T_r,i_d,i_f) matrix = self.define_matrix(samp_time,i_f) raw = self.solve_lin_system(matrix,samp_sig) out[:, i_f, i_d] = raw[:self.out_dim] @@ -183,3 +207,25 @@ def IsWellConditioned(self,m): # with this I am trying to avoid inverting a matr def residuals_func(x,M,S_i): x_cmplx=x[0:int(x.size/2)] + 1j * x[int(x.size/2):x.size] return np.linalg.norm(np.dot(M, x_cmplx) - S_i) +#============================== + + + + def get_sampling(self,idir,ifrq): + # + i_t_start = int(np.round(T_range[0]/self.T_step)) + i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) + T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) + if self.l_out_current: + S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) # **CURRENT + else: + S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) + return T_i,S_i + + def get_sampling(self,idir,ifrq): + # + i_t_start = int(np.round(T_range[0]/self.T_step)) + i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) + T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) + S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) + return T_i,S_i diff --git a/yambopy/nl/pulse_analysis.py b/yambopy/nl/pulse_analysis.py index e6d1f5fd..f22a43c9 100644 --- a/yambopy/nl/pulse_analysis.py +++ b/yambopy/nl/pulse_analysis.py @@ -54,7 +54,7 @@ def build_map(self): I[j+L[self.X_order%2],:] = [-n,m] return I - def get_sampling(self,idir,ifrq): + def set_sampling(self,ifrq): M = int(1+2*self.X_order + int(self.X_order*(self.X_order-1)/2)) if self.nsamp == -1 or self.nsamp < M: self.nsamp = M @@ -64,11 +64,7 @@ def get_sampling(self,idir,ifrq): if (out_of_bounds): print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") - i_t_start = int(np.round(T_range[0]/self.T_step)) - i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) - T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) - S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) - return T_i,S_i + return T_range def update_time_range(self,T_period): T_range = self.T_urange diff --git a/yambopy/nl/sin_analysis.py b/yambopy/nl/sin_analysis.py index 11c18b03..0f996ca7 100644 --- a/yambopy/nl/sin_analysis.py +++ b/yambopy/nl/sin_analysis.py @@ -33,7 +33,7 @@ def set_defaults(self): self.samp_mod="linear" return - def get_sampling(self,idir,ifrq): + def set_sampling(self,ifrq): samp_order = 2*self.X_order + 1 if self.nsamp == -1: self.nsamp = samp_order @@ -42,14 +42,7 @@ def get_sampling(self,idir,ifrq): if (out_of_bounds): print(f'User range redifined for frequency {self.freqs[ifrq]* ha2ev:.3e} [eV]') print(f"Time range: {T_range[0] / fs2aut:.3f} - {T_range[1] / fs2aut:.3f} [fs]") - i_t_start = int(np.round(T_range[0]/self.T_step)) - i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) - T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) - if self.l_out_current: - S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) # **CURRENT - else: - S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) - return T_i,S_i + return T_range def update_time_range(self,T_period): # not sure if this is a general or specific method - let it here for the moment T_range = self.T_urange From 45bd9d0140aa8218535f04f7480042e0b6f9f85b Mon Sep 17 00:00:00 2001 From: myrta75 Date: Thu, 15 Jan 2026 12:11:06 +0000 Subject: [PATCH 22/23] modified: yambopy/nl/nl_analysis.py extra code deleted --- yambopy/nl/nl_analysis.py | 32 ++------------------------------ 1 file changed, 2 insertions(+), 30 deletions(-) diff --git a/yambopy/nl/nl_analysis.py b/yambopy/nl/nl_analysis.py index 81f0b8ef..46b59988 100644 --- a/yambopy/nl/nl_analysis.py +++ b/yambopy/nl/nl_analysis.py @@ -80,13 +80,9 @@ def set_defaults(self): pass @abstractmethod - def set_sampling(self,idir,ifrq): + def set_sampling(self,ifrq): pass - @abstractmethod - def get_sampling(self): - pass - @abstractmethod def define_matrix(self,samp,ifrq): pass @@ -99,8 +95,7 @@ def update_time_range(self): def get_Unit_of_Measure(self,i_order): pass - def get_sampling(self,T_range,idir,ifrq): # subtract self.efield["Initial_time"] from T_i CHECK!!!! - # + def get_sampling(self,T_range,idir,ifrq): i_t_start = int(np.round( T_range[0] / self.T_step)) i_deltaT = int(np.round((T_range[1]-T_range[0])/self.T_step)/self.nsamp) if self.samp_mod=='linear': @@ -117,7 +112,6 @@ def get_sampling(self,T_range,idir,ifrq): # subtract self.efield["Initial_time"] S_i = np.array([self.polarization[ifrq][idir,i] for i in i_t]) return T_i,S_i - def solve_lin_system(self,mat,samp,init=None): mat_dim = int(mat.shape[1]) out=np.zeros(mat_dim,dtype=np.cdouble) @@ -207,25 +201,3 @@ def IsWellConditioned(self,m): # with this I am trying to avoid inverting a matr def residuals_func(x,M,S_i): x_cmplx=x[0:int(x.size/2)] + 1j * x[int(x.size/2):x.size] return np.linalg.norm(np.dot(M, x_cmplx) - S_i) -#============================== - - - - def get_sampling(self,idir,ifrq): - # - i_t_start = int(np.round(T_range[0]/self.T_step)) - i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) - T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) - if self.l_out_current: - S_i = np.array([self.current[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) # **CURRENT - else: - S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) - return T_i,S_i - - def get_sampling(self,idir,ifrq): - # - i_t_start = int(np.round(T_range[0]/self.T_step)) - i_deltaT = int(np.round(T_period/self.T_step)/self.nsamp) - T_i = np.array([(i_t_start + i_deltaT * i) * self.T_step - self.efield["initial_time"] for i in range(self.nsamp)]) - S_i = np.array([self.polarization[ifrq][idir,i_t_start + i_deltaT * i] for i in range(self.nsamp)]) - return T_i,S_i From 9348ae795b79f7d92a52bee0170972e7dd96cd44 Mon Sep 17 00:00:00 2001 From: Myrta Date: Thu, 15 Jan 2026 12:45:41 +0000 Subject: [PATCH 23/23] modified: yambopy/nl/Xn_from_signal_doc.md Upadated documentation to be consistent with changes introduced in the code --- yambopy/nl/Xn_from_signal_doc.md | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/yambopy/nl/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md index fdd734d2..eecddb85 100644 --- a/yambopy/nl/Xn_from_signal_doc.md +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -1,6 +1,6 @@ -# `Xn_from_signal`: documentation (version 1.0) +# `Xn_from_signal`: documentation (version 1.1) #### Myrta Grüning, Claudio Attaccalite, Mike Pointeck, Anna Romani, Mao Yuncheng @@ -49,7 +49,7 @@ The code consists of the abstract class, `Xn_from_signal`, and three subclasses 2. `Xn_from_freqmix`: two monochromatic electric fields 3. `Xn_from_pulse`: a pulse-shaped electric field -The main method in the abstract class `Xn_from_signal` is `perform_analysis` defining the sequence of operations to be performed. This is shown in the diagram below. First, defaults are set for each implementation (`set_defaults`). Then a double loop is entered on field frequencies and directions. For each implementation, the time-dependent signal is sampled (`get_sampling`) and the sampled-signal $P_k$ (`samp_sig`) together with the sampling times $\{t_k\}$ (`samp_time`) is returned. Next, for each implementation, the elements of the matrix $M_{kj}$ (`matrix`) is defined (`define_matrix`). Finally, the linear system is solved (`solve_lin_system`, common to all implementation) and the output passed to the `out` array. The latter is the input of `output_analysis` and `reconstruct_signal` that are implemented in each subclass. +The main method in the abstract class `Xn_from_signal` is `perform_analysis` defining the sequence of operations to be performed. This is shown in the diagram below. First, defaults are set for each implementation (`set_defaults`). Then a loop is entered on field frequencies. For each implementation, the time range `range`is returned for a given frequency (`set_sampling`). Internally, this method also set the number of sampling points `nsamp`in case it is not user-defined. A second loop in entered on the field directions. The time-dependent signal is sampled (`get_sampling`, common to all implementations) and the sampled-signal $P_k$ (`samp_sig`) together with the sampling times $\{t_k\}$ (`samp_time`) is returned. Next, for each implementation, the elements of the matrix $M_{kj}$ (`matrix`) is defined (`define_matrix`). Finally, the linear system is solved (`solve_lin_system`, common to all implementation) and the output passed to the `out` array. The latter is the input of `output_analysis` and `reconstruct_signal` that are implemented in each subclass. ~~~mermaid sequenceDiagram @@ -57,13 +57,16 @@ sequenceDiagram participant Impl as Subclass (implements abstract hooks) Analyzer->>Impl: set_defaults() -loop for each frequency i_f and direction i_d - Analyzer->>Impl: get_sampling(i_d, i_f) - Impl-->>Analyzer: (samp_time, samp_sig) - Analyzer->>Impl: define_matrix(samp_time, i_f) - Impl-->>Analyzer: matrix - Analyzer->>Analyzer: solve_lin_system(matrix, samp_sig) - Analyzer->>Analyzer: out[:, i_f, i_d] = raw[:out_dim] +loop for each frequency i_f + Analyzer->>Impl: set_sampling(i_f) + Impl-->>Analyzer: range + loop for each direction i_d + Analyzer->>Analyzer: get_sampling(range, i_d, i_f) + Analyzer->>Impl: define_matrix(samp_time, i_f) + Impl-->>Analyzer: matrix + Analyzer->>Analyzer: solve_lin_system(matrix, samp_sig) + Analyzer->>Analyzer: out[:, i_f, i_d] = raw[:out_dim] + end end ~~~ ### 1.3 Abstract class diagram @@ -105,7 +108,7 @@ classDiagram %% --- Abstract hooks --- #set_defaults() - #get_sampling(idir, ifrq) + #set_sampling(ifrq) #define_matrix(samp, ifrq) #update_time_range() #get_Unit_of_Measure(i_order) @@ -116,6 +119,7 @@ classDiagram +solve_lin_system(mat, samp, init=None) np.ndarray +perform_analysis() np.ndarray +get_Unit_of_Measure(i_order) float + +get_sampling(range, idir, ifrq) } %% Auxiliary (module-level) functions for completeness