Nonlinear optical response of an extended system¶
A short introduction¶
The nonlinear optical (NLO) response of an extended system can be obtained from its ground state electronic structure. Due to large computional cost, the local field effect effect is neglected. Hence, the NLO response is only frequencydependent. The NLO response can be descibed by nonlinear susceptibility tensors. There are numerious NLO processes, depending on the number of photons involved in the process.
The second harmonic generation (SHG) is a NLO process, in which two photons at the frequency of \(\omega\) generate a photon at the frequency of \(2\omega\). The SHG response is characterized by the secondorder (quadratic) susceptibility tensor, defined via
where \(\alpha,\beta={x,y,z}\) denotes the Cartezain coordinates, and \(E_{\alpha}\) and \(E_{\alpha}\) are the polarization and electric fields, respectivly. For bulk systems, \(\chi_{\gamma\alpha\beta}^{(2)}\) is expressed in the units of m/V. The details of computing \(\chi_{\gamma\alpha\beta}^{(2)}\) is documented in Ref. 1.
Example 1: SHG spectrum of semiconductor: Monolayer MoS2¶
To compute the SHG spectrum of given structure, 3 steps are performed:
Ground state (GS) calculation
import numpy as np from ase.build import mx2 from gpaw import GPAW, FermiDirac from gpaw.nlopt.matrixel import make_nlodata from gpaw.nlopt.shg import get_shg # Make the structure and add the vaccum around the layer atoms = mx2(formula='MoS2', a=3.16, thickness=3.17, vacuum=5.0) atoms.center(vacuum=15, axis=2) # GPAW parameters: nk = 40 params_gs = dict( mode='lcao', symmetry={'point_group': False, 'time_reversal': True}, nbands='nao', convergence={'bands': 10}, parallel={'domain': 1}, occupations=FermiDirac(width=0.05), kpts={'size': (nk, nk, 1), 'gamma': True}, xc='PBE', txt='gs.txt') # Ground state calculation: gs_name = 'gs.gpw' calc = GPAW(**params_gs) atoms.calc = calc atoms.get_potential_energy() atoms.calc.write(gs_name, mode='all')In this script a normal ground state calculation is performed with coarse kpoint grid. Note that LCAO basis is used here, but the PW basis set can also be used. For a smoother spectrum, a finer mesh should be employed.
Get the required matrix elements from the GS
Here, the matrix elements of momentum are computed. Then, all required quantities such as energies, occupations, and momentum matrix elements are saved in a file (‘mml.npz’). The GS file cane be removed after this step.
# Calculate momentum matrix: mml_name = 'mml.npz' make_nlodata(gs_name=gs_name, out_name=mml_name)
Note that, two optional paramters are available in the
gpaw.nlopt.matrixel.make_nlodata()
function:ni
andnf
as the first and last bands used for calculations of SHG.Compute the SHG spectrum
In this step, the SHG spectrum is calculated using the saved data. There are two wellknown gauges that can be used: length gauge or velocity gauge. Formally, they are equivalnent but they may generate different results. Here, SHG susceptibility is computed in both gauges and saved. The SHG susceptibility is a rank3 symmteric tensor with at most 18 independent components. In addition, the point group symmtery reduce the number of independent tensor elements. Monolayer MoS2 has only one independent tensor element: yyy. A broadening is necessary to obtain smooth graphs, and here 50 meV has been used.
# Shift parameters: eta = 0.05 # Broadening in eV w_ls = np.linspace(0, 6, 500) # in eV pol = 'yyy' # LG calculation shg_name1 = 'shg_' + pol + '_lg.npy' get_shg( freqs=w_ls, eta=eta, pol=pol, gauge='lg', out_name=shg_name1, mml_name=mml_name) # VG calculation shg_name2 = 'shg_' + pol + '_vg.npy' get_shg( freqs=w_ls, eta=eta, pol=pol, gauge='vg', out_name=shg_name2, mml_name=mml_name)
Result¶
Now the calculated SHG spectra are plotted at the end. Both real and imaginary parts of the computed SHG susceptibilities, obtained from two gauges are shown. The gauge invariance is confirmed from the calculation. Note that the bulk susceptibility (with SI units of m/V) is illdefined for 2D materials, since the volume cannot be defined without ambiguity in 2D systems. Instead, the sheet susceptibility, expressed in unit of m\(^2\)/V, is an unambiguous quantity for 2D materials. Hence, the bulk susceptibility is transformed to the unambiguous sheet susceptibility by multiplying with the width of the unit cell in the \(z\)direction.
# Creates: shg.png
import numpy as np
import matplotlib.pyplot as plt
from ase.io import read
# Plot and save both spectra
atoms = read('gs.txt')
cell = atoms.get_cell()
cellsize = atoms.get_cell_lengths_and_angles()
mult = cellsize[2] * 1e10 # make the sheet sus.
legls = []
res_name = ['shg_yyy_lg.npy', 'shg_yyy_vg.npy']
plt.figure(figsize=(6.0, 4.0), dpi=300)
for ii, name in enumerate(res_name):
# Load the data
shg = np.load(name)
w_l = shg[0]
plt.plot(np.real(w_l), np.real(mult * shg[1] * 1e18), '')
plt.plot(np.real(w_l), np.imag(mult * shg[1] * 1e18), '')
legls.append(f'{name}: Re')
legls.append(f'{name}: Im')
plt.xlabel(r'$\hbar\omega$ (eV)')
plt.ylabel(r'$\chi_yyy$ (nm$^2$/V)')
plt.legend(legls, ncol=2)
plt.tight_layout()
plt.savefig('shg.png', dpi=300)
The figure shown here is generated from scripts:
shg_MoS2.py
and shg_plot.py
.
It takes 30 minutes with 16 cpus on Intel Xeon X5570 2.93GHz.
Technical details¶
There are few points about the implementation that we emphasize:
The code is parallelized over kpoints.
The code employs only the time reversal symmtery to improve the speed.
Refer to 1 for the details on the NLO theory.
Useful tips¶
Use dry_run option to get an overview of a calculation (especially useful for heavy calculations!):
$ python3 filename.py gpaw=dfdryrun=8
It’s important to converge the results with respect to:
nbands
nkpt
(number of kpoints in gs calc.)
eta
ecut
ftol
vacuum (if there is)
API¶

gpaw.nlopt.matrixel.
make_nlodata
(gs_name: str = 'gs.gpw', out_name: str = 'mml.npz', spin: str = 'all', ni: int = 0, nf: int = 0) → None[source]¶ Get all required NLO data and store it in a file.
Writes NLO data to file: w_sk, f_skn, E_skn, p_skvnn.
Parameters:
 gs_name:
Ground state file name
 out_name:
Output filename
 spin:
Which spin channel (‘all’, ‘s0’ , ‘s1’)
 ni:
First band to compute the mml.
 nf:
Last band to compute the mml (relative to number of bands for nf <= 0).

gpaw.nlopt.shg.
get_shg
(freqs=[1.0], eta=0.05, pol='yyy', eshift=0.0, gauge='lg', ftol=0.0001, Etol=1e06, band_n=None, out_name='shg.npy', mml_name='mml.npz')[source]¶ Calculate RPA SHG spectrum for nonmagnetic semiconductors.
Output: shg.npy file with numpy array containing the spectrum and frequencies.
Parameters:
 freqs:
Excitation frequency array (a numpy array or list)
 eta:
Broadening, a number or an array (default 0.05 eV)
 pol:
Tensor element (default ‘yyy’)
 gauge:
Choose the gauge (lg or vg)
 Etol, ftol:
Tol. in energy and fermi to consider degeneracy
 band_n:
List of bands in the sum (default 0 to nb)
 out_name:
Output filename (default ‘shg.npy’)
 mml_name:
The momentum filename (default ‘mml.npz’)

gpaw.nlopt.linear.
get_chi_tensor
(freqs=[1.0], eta=0.05, ftol=0.0001, Etol=1e06, eshift=0.0, band_n=None, mml_name='mml.npz', out_name=None)[source]¶ Calculate full linear susceptibility tensor for nonmagnetic semiconductors.
 Input:
freqs Excitation frequency array (a numpy array or list) eta Broadening, a number or an array (default 0.05 eV) Etol, ftol Tol. in energy and fermi to consider degeneracy eshift Bandgap correction band_n List of bands in the sum (default 0 to nb) out_name If it is given: output filename mml_name The momentum filename (default ‘mml.npz’)
 Output:
chi_vvl The output tensor (3, 3, nw) chi.npy If specified: array containing the spectrum and freqs