Sum frequency generation

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

X_order=4
NLDB=YamboNLDB(calc='TEST')
SF_Harmonic_Analysis(NLDB,X_order=X_order,prn_Peff=True,prn_Xhi=True)

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....

notice that in the script we set 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:

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