diff --git a/yambopy/__init__.py b/yambopy/__init__.py index 5cf68556..ca002da8 100644 --- a/yambopy/__init__.py +++ b/yambopy/__init__.py @@ -119,7 +119,10 @@ 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 * +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/Xn_from_signal_doc.md b/yambopy/nl/Xn_from_signal_doc.md new file mode 100644 index 00000000..eecddb85 --- /dev/null +++ b/yambopy/nl/Xn_from_signal_doc.md @@ -0,0 +1,230 @@ + + +# `Xn_from_signal`: documentation (version 1.1) + +#### 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](# 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). + +--- + +## 0. Minimal theoretical compendium + + +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. + +--- + +## 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 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 + participant Analyzer as Xn_from_signal + participant Impl as Subclass (implements abstract hooks) + +Analyzer->>Impl: set_defaults() +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 + +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() + #set_sampling(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 + +get_sampling(range, idir, ifrq) + } + + %% 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 +``` +--- +## 2. 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.[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) diff --git a/yambopy/nl/external_efield.py b/yambopy/nl/external_efield.py index 336b093f..7a6ca76b 100644 --- a/yambopy/nl/external_efield.py +++ b/yambopy/nl/external_efield.py @@ -20,22 +20,28 @@ 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) + 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 new file mode 100644 index 00000000..c994c1d7 --- /dev/null +++ b/yambopy/nl/freqmix_analysis.py @@ -0,0 +1,146 @@ +# 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 +from scipy.ndimage import uniform_filter1d +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): + # + 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_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}") + 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"] != "none": + self.pump_freq = self.pumps[0]["freq_range"][0] + 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] + 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]=self.time[-1] + return T_range + # + 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() + return T_range + + def define_matrix(self,T_i,ifrq): + NX,MX = self.X_order[:] + W1 = self.freqs[ifrq] + W2 = self.pump_freq + 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) + return M + + def output_analysis(self,out,to_file=True): + # + NX,MX = self.X_order[:] + T, _ = self.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))): + 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[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{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=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): + 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 + # + # 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 new file mode 100644 index 00000000..46b59988 --- /dev/null +++ b/yambopy/nl/nl_analysis.py @@ -0,0 +1,203 @@ +# Copyright (c) 2023-2026, Claudio Attaccalite, +# 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 +from scipy.optimize import least_squares +import sys +import os +from abc import ABC,abstractmethod +# +# +# 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,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 + 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 + 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 + 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 + 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_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 + self.debug_mode = debug_mode + 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.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)+" [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." + else: + s+="Output is set to polarization." + return s + + @abstractmethod + def set_defaults(self): + pass + + @abstractmethod + def set_sampling(self,ifrq): + pass + + @abstractmethod + def define_matrix(self,samp,ifrq): + pass + + @abstractmethod + def update_time_range(self): + pass + + @abstractmethod + def get_Unit_of_Measure(self,i_order): + pass + + 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': + 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]) + out=np.zeros(mat_dim,dtype=np.cdouble) + 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) + if self.solver=="lstsq": + out = np.linalg.lstsq(mat,samp,rcond=self.tol)[0] + if self.solver=="svd": + 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): + _ = 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(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] + 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 + linear = 1.0 + 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 + 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 + + @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(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 + + 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/pulse_analysis.py b/yambopy/nl/pulse_analysis.py new file mode 100644 index 00000000..f22a43c9 --- /dev/null +++ b/yambopy/nl/pulse_analysis.py @@ -0,0 +1,164 @@ +# 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 +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 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' + 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)[0] + self.sigma = Gaussian_centre(self.efield)[1] + if self.samp_mod== "": + self.samp_mod="linear" + 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,self.X_order+1,2): + j += 1 + I[j,:] = [n,floor(n/2)] + for n in range(1,self.X_order+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 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 + 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): + 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]") + return T_range + + 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 - self.sigma + T_range[1] = T_range[0] + T_period + if T_range[1] > self.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 + + 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.T0)**2/(2*self.sigma**2)) #sign + if (n%2 == 0 and i_m == int(n/2)): + M[:,ii] = M[:,ii] + 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): + 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 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*self.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) + T = 10000.0 + for i_f in range(self.n_runs): + 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')]) + 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") + 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): + 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 new file mode 100644 index 00000000..0f996ca7 --- /dev/null +++ b/yambopy/nl/sin_analysis.py @@ -0,0 +1,120 @@ +# 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 +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_sine(Xn_from_signal): + 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}") + 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 + if self.samp_mod== "": + self.samp_mod="linear" + return + + def set_sampling(self,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]") + 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 + 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): + 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) + 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')]) + 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") + 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): + 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) + 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 + 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