 |
(git:7fc1fe1)
|
|
(git:7fc1fe1)
|
Functions/Subroutines | |
| subroutine, public | run_nlo (qs_env) |
| Calculate frequency-dependent χ^(2) to χ^(5) after the ground-state SCF. | |
| subroutine, public nlo::run_nlo | ( | type(qs_environment_type), pointer | qs_env | ) |
Calculate frequency-dependent χ^(2) to χ^(5) after the ground-state SCF.
| qs_env | Converged k-point SCF environment |
Method: length-gauge density-matrix response in the spirit of Aversa and Sipe, Phys. Rev. B 52, 14636 (1995), DOI: 10.1103/PhysRevB.52.14636; analytic covariant k derivatives follow Ventura et al., Phys. Rev. B 96, 035431 (2017), DOI: 10.1103/PhysRevB.96.035431. The finite AO projection, analytic position curvature, and numerical implementation are defined below and in nlo_current_tensor_derivatives.
Build lattice-image matrices for the AO Hamiltonian H, overlap S, dipole d. d_μν(R) = ⟨ϕ_(μ0)|r|ϕ_(νR)⟩ is a position matrix; the electron charge enters the later SI conversion. At each optical k point:
X(k) = Σ_R exp(2π i k_frac · R) X(R), with X = H, S, d.
Differentiating the phase gives every required k derivative. The AO overlap is removed, then H and position matrices are expressed in the eigenvector basis of H at the current k point.
At each optical k point, diagonalize the Hamiltonian once. The occupied bands define the equilibrium density matrix. A light field changes that density matrix: first order uses the equilibrium matrix, second order uses first-order responses, and so on through fifth order. Each new response is obtained from the previous one by an analytic k derivative and a band-energy denominator ε_m - ε_n - Ω - iΓ_mn. Here ε_m and ε_n are band energies, Ω is the sum of the incident photon energies, and Γ is the band-pair relaxation rate. For T₁ = T₂ it is the constant 1/T₁. For separate times, with γᵢ = 1/Tᵢ and Δε = ε_m - ε_n, it is
Γ_mn = γ₁ + (γ₂ - γ₁)[1 - exp(-(Δε/w)²)].
The positive width w sets the smooth transition at nearly degenerate bands. Its basis-covariant operator form is Γ_H(X) = γ₂X + (γ₁ - γ₂)exp[-(ad_H/w)²]X, where ad_H(X) = [H,X]. Its analytic k derivatives enter the response. All possible orders of the incident fields are averaged. No numerical derivatives of separately diagonalized bands or their phases are used.
The position operator in the band representation is
r_a = i ∂/∂k_a + A_a .
The subscript a denotes x, y or z. A_a is the matrix obtained from the AO dipole matrix after orthonormalizing and rotating into the central band basis. It is the finite-basis part of the position operator. When a field changes the current, an additional term from this matrix is required:
F_ia = ∂A_a/∂k_i - ∂A_i/∂k_a - i [A_i, A_a] .
Here i is the output-current direction. F_ia is the position curvature; its contribution to the current is always included. The other current contribution is the band-velocity matrix multiplied by the density response. Summing over optical k points and spin gives the cell-integrated current response I^(n) in atomic units, stored locally on each MPI rank.
Sum I^(n) across MPI ranks, convert it to SI units and divide by the simulation-cell volume V. The result is the nonlinear conductivity σ^(n) for each field order. For an output frequency Ω, the corresponding polarization obeys J(Ω) = -i Ω P(Ω). Hence:
χ²_iab = σ²_iab / [-i (ω₁+ω₂) ε₀] χ³_iabc = σ³_iabc / [-i (ω₁+ω₂+ω₃) ε₀] χ⁴_iabcd = σ⁴_iabcd / [-i (ω₁+ω₂+ω₃+ω₄) ε₀] χ⁵_iabcde = σ⁵_iabcde / [-i (ω₁+ω₂+ω₃+ω₄+ω₅) ε₀] .
The SI conversion is in nlo_write_susceptibility. For a two-dimensional slab of in-plane area A, the sheet response is χ_sheet = (V/A) χ_volume. Ω must be nonzero.
Definition at line 144 of file nlo.F.
1.15.0