Sum frequency generation

This tutorial is for internal use only, these response functions are not implemented/tested in yambo/yambopy suite.

In this tutorial we will show how to calculate Sum Frequency Generation(SFG) and also Difference Frequency Generation (DFG) in bulk materials.
We suppose you are already familiar with the non-linear response using the Yambo code. If it is not the case please study the previous tutorials:
Linear response using Dynamical Berry Phase and Real time approach to non-linear response (SHG).
We thank Mike N. Pionteck for the script used in this tutorial.

DFT calculations

In this example, we will consider a single layer of hexagonal boron nitride (hBN). If you didn't before you can download input files and Yambo databases for this tutorial here: hBN-2D-RT.tar.gz. and/or follow the instructions to generate the databases here: Prerequisites for Real Time propagation with Yambo

Removing symmetries

In this tutorial we will calculate the SFG along when both fields are in the 'y' direction, therefore we remove symmetries not compatible with an external field along this direction, with the command ypp_nl -y:

fixsyms                          # [R] Remove symmetries not consistent with an 
external perturbation
% Efield1
  0.000000 | 1.000000 | 0.000000 |        # First external Electric Field
%
% Efield2
 0.000000 | 0.000000 | 0.000000 |        # Additional external Electric Field
%
BField= 0.000000           T     # [MAG] Magnetic field modulus
Bpsi= 0.000000             deg   # [MAG] Magnetic field psi angle [degree]
Btheta= 0.000000           deg   # [MAG] Magnetic field theta angle [degree]
#RmAllSymm                     # Remove all symmetries
RmTimeRev                     # Remove Time Reversal
#RmSpaceInv                    # Remove Spatial Inversion

Real-time simulation with two external fields

You go in the FixSymm folder and run again the setup. Then you can put the following input file, that has been generated with the command yambo_nl -u n , in the folder with the name yambo.in_sfg:

nloptics                         # [R] Non-linear spectroscopy
NLogCPUs=0                       # [PARALLEL] Live-timing CPU`s (0 for all)
PAR_def_mode= "balanced"         # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload"/"KQmemory")
NL_CPU= "10 1"                       # [PARALLEL] CPUs for each role
NL_ROLEs= "w k"                     # [PARALLEL] CPUs roles (w,k)
DIP_CPU= ""                      # [PARALLEL] CPUs for each role
DIP_ROLEs= ""                    # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0                    # [OPENMP/X] Number of threads for dipoles
NL_Threads=0                     # [OPENMP/NL] Number of threads for nl-optics
% NLBands
   3 |  6 |                           # [NL] Bands range
%
NLverbosity= "high"              # [NL] Verbosity level (low | high)
NLtime= 60.00000           fs    # [NL] Simulation Time
NLintegrator= "CRANKNIC"         # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC")
NLCorrelation= "IPA"             # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX")
NLLrcAlpha= 0.000000             # [NL] Long Range Correction
% NLEnRange
  2.000000 | 8.000000 |         eV    # [NL] Energy range (for loop on frequencies NLEnSteps/=0
%
 NLEnSteps= 30                  # [NL] Energy steps for the loop on frequencies
% NLrotaxis
 0.000000 | 0.000000 | 0.000000 |        # [NL] Rotation axis (for the loop on angles NLAngSteps/=0)
%
NLAngSteps=0                     # [NL] Angular steps (if NLAngSteps/=0 field versor will be ignored)
NLDamping= 0.200000        eV    # [NL] Damping (or dephasing)
RADLifeTime=-1.000000      fs    # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)
#EvalCurrent                   # [NL] Evaluate the current
#FrPolPerdic                   # [DIP] Force periodicity of polarization respect to the external field
HARRLvcs= 18475            RL    # [HA] Hartree     RL components
EXXRLvcs= 18475            RL    # [XX] Exchange    RL components
 % Field1_Freq
-0.100000 |-0.100000 |         eV    # [RT Field1] Frequency
%
Field1_NFreqs= 1                 # [RT Field1] Frequency
Field1_Int=   1000.00       kWLm2 # [RT Field1] Intensity
Field1_Width= 0.000000     fs    # [RT Field1] Width
Field1_kind= " SOFTSIN"           # [RT Field1] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field1_pol= "linear"             # [RT Field1] Pol(linear|circular)
% Field1_Dir
  0.000000 | 1.000000 | 0.000000 |        # [RT Field1] Versor
%
Field1_Tstart= 0.010000    fs    # [RT Field1] Initial Time
% Field2_Freq
3.000000 | 3.000000 |         eV    # [RT Field2] Frequency
%
Field2_NFreqs= 1                 # [RT Field2] Frequency
Field2_Int=   1000.00       kWLm2 # [RT Field2] Intensity
Field2_Width= 0.000000     fs    # [RT Field2] Width
Field2_kind= " SOFTSIN"           # [RT Field2] Kind(SIN|SOFTSIN| see more on src/modules/mod_fields.F)
Field2_pol= "linear"             # [RT Field2] Pol(linear|circular)
% Field2_Dir
  0.000000 | 1.000000 | 0.000000 |        # [RT Field2] Versor 
%
Field2_Tstart= 0.010000    fs    # [RT Field2] Initial Time

If you run this input file with the command yambo_nl -F yambo.in_sfg, the code will run 30 simulations for a first laser field with frequency between 2.0 and 8.0 while the frequency of the second field is fixed at 3.0 eV. For this reason we provide you a simple python script to change also the frequency of the second field (in blue in the input).
In the next section we will show how to analyse the non-linear response generated by the presence of two external fields.

Analysis of the single run

We will use an extension to two fields of the approach described in Ref. and used for the tutorials on second/third harmonic generation: Real time approach to non-linear response (SHG). These new subroutine are available in the new YamboPy package. Hereafter the pythonn script for the analysis:

from yambopy import *
from yambopy.nl.sum_frequencies import SF_Harmonic_Analysis
from yambopy.units import fs2aut
X_order=4
NLDB=YamboNLDB(calc='TEST')
SF_Harmonic_Analysis(NLDB,X_order=X_order,T_range=[50.0*fs2aut,-1.0],prn_Peff=True,prn_Xhi=True)

Notice that we set the beginning of time range where analysis is performed with the variable T_range=[50.0*fs2aut,-1.0], you can also omit this variable and the code automatically will set the optimal range for the analysis using the damping given in input in the real-time simulation.
This script will analyse the real-time response and produce all the different [math]\displaystyle{ \chi^{(n)}(\omega) }[/math] in the form:

o.YamboPy-SF_probe_order_1_0 ==> [math]\displaystyle{ \chi^{(1)}(\omega=\omega_1) }[/math]
o.YamboPy-SF_probe_order_1_1 ==> [math]\displaystyle{ \chi^{(2)}(\omega=\omega_1+\omega_2) }[/math]
o.YamboPy-SF_probe_order_2_0 ==> [math]\displaystyle{ \chi^{(2)}(\omega=2\omega_1) }[/math]
o.YamboPy-SF_probe_order_0_2 ==> [math]\displaystyle{ \chi^{(2)}(\omega=2\omega_2) }[/math]
o.YamboPy-SF_probe_order_1_2 ==> [math]\displaystyle{ \chi^{(3)}(\omega=\omega_1+2\omega_2) }[/math]
o.YamboPy-SF_probe_order_2_2 ==> [math]\displaystyle{ \chi^{(4)}(\omega=2\omega_1+2\omega_2) }[/math]
etc....

for example we can plot the 2 and 3 columns of the file o.YamboPy-SF_probe_order_1_1 that corresponds to the [math]\displaystyle{ \chi^{(2)}_{xyy} (\omega=\omega_1+\omega_2) }[/math] for [math]\displaystyle{ \omega_2 =3.0~eV }[/math]:

Here the python script used to generate the previous figure: plot_xhi2_monolayer.py

Checking the Fourier analysis

In the previous script we turned on the flag prn_Peff=True that automatically will produce a series of file to check if the Fourier was successful. The files o.YamboPy-sampling_FX contain the sampling points used to fix the real-time polarization, you can for example plot o.YamboPy-sampling_F1 vs o-TEST.polarization_F1 together to see which part of [math]\displaystyle{ P(t) }[/math] is used to extract the response functions:

Then you can check if real-time polarization is well reconstructed by the Fourier coefficients calculated by the script by plotting o.YamboPy-pol_reconstructed_F1 vs o-TEST.polarization_F1

As you see the polarization is well reconstructed starting from the extracted Fourier coefficients, expect for the beginning, but this is normal because the initial signal is affected by the internal eigen--modes of the systems that later are washed out by the de-phasing. For this reason we use only the last part of the signal for the analysis, for more details see Ref.

Analysis of all runs

For this reason we provide you a simple python script to change also the frequency of the second field (in blue in the input).

run_many_sfg.py

please modify the python script by set the correct path of your executable, the parallelization and the number of frequency steps you are interested in. In this example we will do 30 frequency steps for both the first (in the Yambo input) and the second field (in the python script).

Calculations will take some time, in the above example we parallelized on 10 cores to speed up them.
When the calculations are terminated you will find many folders called EF1, EF2, ... EF29. These folders are generated from the loop on the second field frequency(Field2_Freq), and each one contains the response real-time polarization at 30 different frequencies of the first external field (Field1_Freq).

Now we will use a script to read all the real-time polarization, Fourier analyze them as it was done in the previous section for a single one, and the write the [math]\displaystyle{ \chi^{(2)}(\omega=\omega_1 + \omega_2) }[/math] in a file in the form: [math]\displaystyle{ \omega_1 , \omega_2, Re \{\chi^{(2)}\} , Im \{\chi^{(2)}\} }[/math].

analysis_sfg.py

then we can interpolate the results [math]\displaystyle{ \chi^{(2)}(\omega) }[/math] and plot results in 3D or 2D: plot_xhi_2d.py, plot_xhi_3d.py

Comparison with SHG

Finally you can compare the diagonal of [math]\displaystyle{ \chi^{(2)}(\omega=\omega_1 +\omega_2) }[/math], when [math]\displaystyle{ \omega_1 =\omega_2 }[/math], with the standard second harmonic generation [math]\displaystyle{ \chi^{(2)}(\omega=2 \omega_1) }[/math] calculated as explained in the tutorial Real time approach to non-linear response (SHG):

Notice that there is factor two between the two response functions, this is the so-called degeneracy factor, for more details see Sec. 1.3 in Ref.. The small differences between the two spectra are given the different sampling strategies of the SF_Harmonic_Analysis function respect to the Harmonic_Analysis one.