Quadratic response#
First-order electric-dipole hyperpolarizability#
Tensor components of the first-order electric-dipole hyperpolarizability relate to quadratic response functions according to
where \(\omega_\sigma = \omega_1 + \omega_2\). In the limit of static fields, \(\beta(0;0,0)\) can also be determined from the third-order derivative of the energy with respect to electric field strengths.
Pertinent optical processes#
For dynamic fields, nonlinear optical processes of interest relating to this molecular property include:
Optical process |
Response function |
|---|---|
Second Harmonic Generation |
\(\beta^{\rm SHG}_{\mu\nu\eta}(-2\omega;\omega,\omega) =\langle\langle \hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta} \rangle\rangle_{\omega,\omega}\) |
Electro-Optic Pockels Effect |
\(\beta^{\rm EOPE}_{\mu\nu\eta}(-\omega;\omega,0) = \langle\langle\hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta}\rangle\rangle_{\omega,0}\) |
Optical Rectification |
\(\beta^{\rm OR}_{\mu\nu\eta}(0;\omega,-\omega) = \langle\langle\hat{\mu}_{\mu};\hat{\mu}_{\nu},\hat{\mu}_{\eta}\rangle\rangle_{\omega,-\omega}\) |
We note that the EOPE is also known under the name of the dc-Pockels effect.
Symmetry aspects#
The general definition of the quadratic response function, indicates its symmetry with respect to permutation of operators. Thus, for all the dipole hyperpolarizabilities we have
Also the relation
holds, and we have, for instance, \(\beta^{\rm OR}_{\mu\nu\eta}(0;\omega,-\omega)\) = \(\beta^{\rm EOPE}_{\eta\nu\mu}(-\omega;\omega,0)\). We note that this relation does not hold for macroscopic samples due to the differences in the local fields experienced by the molecule in these two different experimental set-ups.
In the case of the first hyperpolarizability, the measured quantity is the vector component of \(\beta\) in the direction of the permanent dipole moment \(\mu\), defining the molecular \(z\) axis. The relevant averages are given by
where the same sequence of optical frequencies (not given explicitly) is used for the laboratory axes and the molecular axes.
The number of independent non-zero tensor elements depends on the non-linear optical process and on the symmetry of the molecule. For example, \(\beta^{\rm SHG}\) is symmetric with respect to the permutation of the second and third indices and this can be used to simplify the equations for the parallel and perpendicular components. For non-zero frequencies, the number of independent tensor components to be computed decreases when Kleinman’s symmetry is assumed, that is, we assume that we can permute the indices of the incoming light without changing the corresponding frequencies,
Although for \(\omega_1 \neq \omega_2\) this is only an approximation and these two tensor components are not equal, Kleinman’s symmetry is often applied for low frequencies, where it it is usually found to be a reasonable approximation.
The basis set requirements for the calculation of first hyperpolarizabilities are much the same as for the linear polarizability. However, as the first hyperpolarizability probes even higher-order electric-field-perturbed densities of the molecule, care should be exercised to ensure that the basis set is sufficiently saturated with respect to diffuse polarizing functions. Special basis sets have been developed for the calculation of hyperpolarizabilities by Pluta and Sadlej, though the same basis sets that can be used for polarizabilities in most cases give reliable estimates also for first hyperpolarizabilities.
Due to the fact that the first hyperpolarizability involves multiple virtual excited states, the property has proven to be more sensitive to electron correlation effects than the polarizability. Molecules with push-pull conjugated structures are often characterized by large first hyperpolarizabilities. Due to the problems facing many commonly used local exchange-correlation functionals, DFT has been shown to have difficulties describing correctly nonlinear optical properties such as the first hyperpolarizability for extended systems}. The Coulomb-attenuated B3LYP functional (CAM-B3LYP) has, in contrast, shown quite good performance also for rather extended systems.
The atomic units of \(\beta\) are e\(^3\) a\(_0^3\) E\(_h^{-2}\). 1 au of \(\beta\) corresponds to 3.20636 \(\times 10^{-53}\) C\(^3\) m\(^3\) J\(^{-2}\) (SI) or to 8.6392 \(\times 10^{-33}\) statvolt\(^{-1}\) cm\(^4\) (esu).
Dispersion#
Similar to Cauchy expansions, expansions may be applied to describe the dispersions of hyperpolarizabilities
Since the expansion coefficient, \({\cal A}^{\beta}\), is frequency independent in itself, dispersions are transferable in between different nonlinear optical processes.
In VeloxChem#
As an illustration, let us calculate \(\beta_{zzz}(-2\omega; \omega, \omega)\) using damped response theory.
import veloxchem as vlx
First, we determine the reference state of the system.
mol_str = """3
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.740848095288 0.582094932012
H 0.000000000000 -0.740848095288 0.582094932012
"""
molecule = vlx.Molecule.read_xyz_string(mol_str)
basis = vlx.MolecularBasis.read(molecule, "def2-svpd", ostream=None)
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
scf_drv.xcfun = "b3lyp"
scf_results = scf_drv.compute(molecule, basis)
Second, we calculate the needed quadratic response functions with use of the QuadraticResponseDriver class.
qrf_drv = vlx.QuadraticResponseDriver()
qrf_drv.a_operator = "electric dipole"
qrf_drv.b_operator = "electric dipole"
qrf_drv.c_operator = "electric dipole"
# other available operators
#lrf.b_operator = "electric dipole"
#lrf.b_operator = "linear momentum"
#lrf.b_operator = "angular momentum"
#lrf.b_operator = "magnetic dipole"
qrf_drv.a_component = "z"
qrf_drv.b_component = "z"
qrf_drv.c_component = "z"
qrf_drv.b_frequencies = [0.0, 0.0656, 0.1312]
qrf_drv.c_frequencies = [0.0, 0.0656, 0.1312]
qrf_drv.damping = 0.004556 # 1000 cm-1
qrf_results = qrf_drv.compute(molecule, basis, scf_results)
Show code cell output
Quadratic Response Driver Setup
=================================
ERI Screening Threshold : 1.0e-12
Convergance Threshold : 1.0e-04
Max. Number of Iterations : 150
Damping Parameter : 4.556000e-03
Exchange-Correlation Functional : B3LYP
Molecular Grid Level : 4
Complex Response Solver Setup
===============================
Number of Frequencies : 1
Max. Number of Iterations : 150
Convergence Threshold : 1.0e-04
ERI Screening Scheme : Cauchy Schwarz + Density
ERI Screening Threshold : 1.0e-12
Exchange-Correlation Functional : B3LYP
Molecular Grid Level : 4
* Info * Using the Libxc library (version 6.2.2).
S. Lehtola, C. Steigemann, M. J.T. Oliveira, and M. A.L. Marques., SoftwareX 7, 1–5 (2018)
* Info * Using the B3LYP functional.
P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch., J. Phys. Chem. 98, 11623 (1994)
* Info * Molecular grid with 40720 points generated in 0.02 sec.
* Info * Processing Fock builds... (batch size: 8)
* Info * batch 1/1
* Info * 4 gerade trial vectors in reduced space
* Info * 4 ungerade trial vectors in reduced space
* Info * 94.37 kB of memory used for subspace procedure on the master node
* Info * 2.32 GB of memory available for the solver on the master node
*** Iteration: 1 * Residuals (Max,Min): 8.44e-01 and 5.68e-01
* Info * Processing Fock builds... (batch size: 8)
* Info * batch 1/1
* Info * 8 gerade trial vectors in reduced space
* Info * 8 ungerade trial vectors in reduced space
* Info * 116.13 kB of memory used for subspace procedure on the master node
* Info * 2.33 GB of memory available for the solver on the master node
*** Iteration: 2 * Residuals (Max,Min): 7.29e-02 and 5.16e-02
* Info * Processing Fock builds... (batch size: 8)
* Info * batch 1/1
* Info * 12 gerade trial vectors in reduced space
* Info * 12 ungerade trial vectors in reduced space
* Info * 137.89 kB of memory used for subspace procedure on the master node
* Info * 2.33 GB of memory available for the solver on the master node
*** Iteration: 3 * Residuals (Max,Min): 2.53e-03 and 1.87e-03
* Info * Processing Fock builds... (batch size: 10)
* Info * batch 1/1
* Info * 17 gerade trial vectors in reduced space
* Info * 17 ungerade trial vectors in reduced space
* Info * 165.09 kB of memory used for subspace procedure on the master node
* Info * 2.33 GB of memory available for the solver on the master node
*** Iteration: 4 * Residuals (Max,Min): 7.36e-05 and 5.73e-05
*** Complex response converged in 4 iterations. Time: 1.69 sec
Fock Matrix Computation
=========================
* Info * Processing Fock builds... (batch size: 6)
* Info * batch 1/1
* Info * Time spent in Fock matrices: 0.76 sec
Quadratic response function: << z;z,z >> (0.0,0.0)
=====================================================
Real Imaginary
----------------------------------------------------
QRF 6.73136841 0.00000000j
Quadratic response function: << z;z,z >> (0.0656,0.0656)
===========================================================
Real Imaginary
----------------------------------------------------
QRF 8.98788527 0.24601088j
Quadratic response function: << z;z,z >> (0.1312,0.1312)
===========================================================
Real Imaginary
----------------------------------------------------
QRF 28.04348163 2.22700923j
*** Time spent in quadratic response calculation: 2.50 sec ***
Third, we print the requested molecular properties.
print(f"{'w1':>8s}{'w2':>8s}{'beta_zzz':>24s}")
print("-"*40)
for key in qrf_results.keys():
print(f"{key[1]:8.4f}{key[2]:8.4f}{-qrf_results[key]:24.6f}")
w1 w2 beta_zzz
----------------------------------------
0.0000 0.0000 -6.731368-0.000000j
0.0656 0.0656 -8.987885-0.246011j
0.1312 0.1312 -28.043482-2.227009j
Note
The response functions in return from VeloxChem are determined for position operators, \(\hat{r}_\alpha\) and since we have \(\hat{\mu}_\alpha = - e \hat{r}_\alpha\), we need to change the sign of the response function values to get results for the first-order hyperpolarizability.