Linear Response

Linear Response#

Electric-dipole polarizability#

The frequency-dependent electric-dipole polarizability is, next to the dipole moment, the most important effect characterizing the response of a molecule to electromagnetic radiation. It is a second-order symmetric tensor determined by the linear response function

\[\begin{equation*} \alpha_{\alpha\beta}(-\omega; \omega) = - \langle\langle \hat{\mu}_{\alpha};\hat{\mu}_{\beta} \rangle\rangle_{\omega} \end{equation*}\]

In the limit of static fields, \(\alpha(0;0)\) can also be determined from the second-order derivative of the energy with respect to electric field strengths.

Significance#

Several molecular properties which describe nonlinear effects, such as the Kerr effect or magnetic circular dichroism arising in the presence of radiation and additional electric or magnetic fields, are interpreted as derivatives of the dipole polarizability. They can be calculated as higher-order response functions. Similarly, relativistic corrections to the polarizabilities for heavy atoms can be estimated from higher-order response functions including the mass-velocity and Darwin operators, [KOJ95].

Isotropic and anisotropic components#

In isotropic fluids, the experimentally measured quantity is
usually the scalar component of the \(\alpha\) tensor given by the isotropic average, defined as

\[\begin{equation*} \alpha_{\rm iso} = \frac{1}{3} \alpha_{\eta\eta} \end{equation*}\]

The polarizability anisotropy may be defined as

\[\begin{eqnarray*} (\alpha_{\rm anis})^2 &=& (3 \alpha_{\mu\nu}\alpha_{\mu\nu} - \alpha_{\mu\mu}\alpha_{\nu\nu})/2, \end{eqnarray*}\]

which, using the components of the diagonalized tensor, leads to

\[\begin{eqnarray*} \alpha_{\rm anis} &=& 2^{-1/2} [ (\alpha_{xx} - \alpha_{yy})^2 + (\alpha_{yy} - \alpha_{zz})^2 + (\alpha_{zz} - \alpha_{xx})^2 ]^{1/2}. \end{eqnarray*}\]

For linear molecules, \(\alpha_{\rm anis}\) reduces to the difference between the parallel and perpendicular components,

\[\begin{eqnarray*} \alpha_{\rm anis} = \alpha_{\parallel} - \alpha_{\perp}. \end{eqnarray*}\]

Measurements#

The macroscopic property related to the molecular polarizability \(\alpha\) is the dielectric constant \(\epsilon\), defined by the ratio of the permittivity of the medium to the electric constant, \({\epsilon_0}\). It is also related to the refractive index \(n\) (\(n=\sqrt{\epsilon}\), assuming the vacuum magnetic permeability to be equal to 1). The quantity measured in a dielectric constant experiment is

\[\begin{equation*} \tilde{\alpha} = \left( \tilde{\mu}_{Z}^{F}/{F} \right)_{F \rightarrow 0} = \langle \alpha_{II} \rangle_\Omega + \mu_z^2/3kT \end{equation*}\]

where we use the notation \( <...>_\Omega \) for the average over all the orientations, \(k\) is the Boltzmann constant and \(T\) is the temperature. The temperature-dependent term is proportional to the square of the dipole moment, the temperature-independent term is the isotropic average of the polarizability (using the equation for the isotropic average given just above we can obtain a corresponding expression in terms of properties calculated in the molecular frame of reference).

The anisotropy of the electric dipole polarizability can be determined in Rayleigh scattering experiments, where the observed depolarization of the incident light is a function of \((\alpha_{\rm anis})^2\).

High-accuracy results#

The static electric dipole polarizability of the He atom has been studied in more detail than any other response property of any other atomic or molecular system. Not only the non-relativistic value is known accurately, and the relativistic corrections have been computed, but also the quantum-electrodynamic corrections have been analyzed and numerous effects related to the nuclear mass have been taken into account [LachJS04].

Basis set requirements#

Standard energy-optimized basis sets are not suitable for accurate calculations of electric polarizabilities. The simplest solution (adding the necessary polarization and diffuse functions) makes the basis sets too large to enable efficient calculations for large molecules. Significantly smaller basis sets, designed considering the electric-field dependence of the orbitals, provide results of similar or better quality at lower computational cost.

A practical problem that may be encountered in calculations using very extended basis sets is the appearance of linear dependencies in the basis set. This may make the response equations (as well as the equations determining the molecular energy/density) ill-conditioned and the calculations slowly convergent; in exceptional cases it may even be impossible to converge the equations. The problem is normally resolved by removing the linear dependencies in the basis set, or by removing manually the most diffuse basis functions from the basis set.

Electron correlation#

Electron correlation effects on the isotropic polarizability are in general moderate, being at the most 5-10%\ of the uncorrelated value and for nondipolar systems always leading to an increased polarizability. Correlation effects are usually larger on the individual tensor components than on the isotropic polarizability, and can therefore be expected to be more important for the polarizability anisotropies. In general, Hartree-Fock polarizability anisotropies cannot be considered accurate enough to allow for a quantitative comparison with experimental observations, or for the interpretation of the ellipticity arising in some birefringences, where the polarizability anisotropies often constitute an important (and dominant) orientational contribution.

DFT in general reproduces quite well correlation effects in the polarizability of small molecular systems. In particular, polarizability anisotropies calculated using DFT will usually be in better agreement with experiment. However, due to the local nature of the exchange-correlation functionals and the lack of general implementations of current density functionals, DFT does not perform so well for extended conjugated systems.

Cauchy moment expansion#

Response functions for a system in its electronic ground state\index{ground state} are analytic functions of the frequency arguments, except at the poles that occur when a frequency or a sum of frequencies is equal to an excitation energy. Thus, for frequencies below the first pole, the linear, quadratic or cubic response functions can be expanded in power series.

Let us consider in more detail the simplest case when only electric dipole operators are involved and the frequency dependence of the (hyper)polarizabilities is of interest. The Cauchy series - that is, the power series expansion of the frequency-dependent polarizability - is usually written as

\[ \alpha(-\omega ; \omega) = \sum_{k=0}^{\infty} \omega^{2k} S(-2k-2) = S(-2) + \omega^2 S(-4) + \omega^4 S(-6) + \cdots, \]

where \(S(-2) = \alpha(0;0)\). This expansion is valid also for purely imaginary frequency arguments.

In VeloxChem#

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)

Real response function#

We calculate the needed linear response functions with use of the LinearResponseSolver class.

lrf = vlx.LinearResponseSolver()

lrf.a_operator = "electric dipole"
lrf.b_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"

lrf.a_components = ["x", "y", "z"]
lrf.b_components = ["x", "y", "z"]

lrf.frequencies = [0.0, 0.0656]

lrf_results = lrf.compute(molecule, basis, scf_results)

Show code cell output Hide code cell output

                                                                                                                          
                                               Linear Response Solver Setup                                               
                                              ==============================                                              
                                                                                                                          
                               Number of Frequencies           : 2                                                        
                               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: 9)                                                                        
* Info *   batch 1/1

                                                                                                                          
* Info * 3 gerade trial vectors in reduced space                                                                          
* Info * 6 ungerade trial vectors in reduced space                                                                        
                                                                                                                          
* Info * 47.77 kB of memory used for subspace procedure on the master node                                                
* Info * 2.34 GB of memory available for the solver on the master node                                                    
                                                                                                                          
               *** Iteration:   1 * Residuals (Max,Min): 8.45e-01 and 1.89e-01                                            
                                                                                                                          
               <<x;x>>_0.0000 :     -9.01161246 Residual Norm: 0.19770107                                                 
               <<x;x>>_0.0656 :     -9.20248108 Residual Norm: 0.18863810                                                 
               <<y;y>>_0.0000 :     -9.53971511 Residual Norm: 0.31867382                                                 
               <<y;y>>_0.0656 :     -9.64430453 Residual Norm: 0.31469367                                                 
               <<z;z>>_0.0000 :     -8.46320379 Residual Norm: 0.84501162                                                 
               <<z;z>>_0.0656 :     -8.57203409 Residual Norm: 0.83255985                                                 
                                                                                                                          

* Info * Processing Fock builds... (batch size: 9)                                                                        
* Info *   batch 1/1                                                                                                      

                                                                                                                          
* Info * 6 gerade trial vectors in reduced space                                                                          
* Info * 12 ungerade trial vectors in reduced space                                                                       
                                                                                                                          
* Info * 72.25 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:   2 * Residuals (Max,Min): 1.12e-01 and 2.01e-02                                            
                                                                                                                          
               <<x;x>>_0.0000 :     -9.09615293 Residual Norm: 0.02138969                                                 
               <<x;x>>_0.0656 :     -9.29211162 Residual Norm: 0.02009095                                                 
               <<y;y>>_0.0000 :     -9.73585054 Residual Norm: 0.04219749                                                 
               <<y;y>>_0.0656 :     -9.84454457 Residual Norm: 0.04171911                                                 
               <<z;z>>_0.0000 :     -9.31931552 Residual Norm: 0.11162972                                                 
               <<z;z>>_0.0656 :     -9.45876238 Residual Norm: 0.11004217                                                 
                                                                                                                          

* Info * Processing Fock builds... (batch size: 8)                                                                        
* Info *   batch 1/1                                                                                                      

                                                                                                                          
* Info * 9 gerade trial vectors in reduced space                                                                          
* Info * 17 ungerade trial vectors in reduced space                                                                       
                                                                                                                          
* Info * 94.01 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:   3 * Residuals (Max,Min): 1.37e-02 and 7.82e-04                                            
                                                                                                                          
               <<x;x>>_0.0000 :     -9.09692448 Residual Norm: 0.00083876                                                 
               <<x;x>>_0.0656 :     -9.29287459 Residual Norm: 0.00078245                                                 
               <<y;y>>_0.0000 :     -9.73799055 Residual Norm: 0.00264362                                                 
               <<y;y>>_0.0656 :     -9.84673034 Residual Norm: 0.00268306                                                 
               <<z;z>>_0.0000 :     -9.34755032 Residual Norm: 0.01367302                                                 
               <<z;z>>_0.0656 :     -9.48835616 Residual Norm: 0.01341692                                                 
                                                                                                                          

* Info * Processing Fock builds... (batch size: 8)                                                                        
* Info *   batch 1/1                                                                                                      

                                                                                                                          
* Info * 12 gerade trial vectors in reduced space                                                                         
* Info * 22 ungerade trial vectors in reduced space                                                                       
                                                                                                                          
* Info * 115.81 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): 6.56e-04 and 5.14e-05                                            
                                                                                                                          
               <<x;x>>_0.0000 :     -9.09692566 Residual Norm: 0.00005136                                                 
               <<x;x>>_0.0656 :     -9.29287577 Residual Norm: 0.00005883                                                 
               <<y;y>>_0.0000 :     -9.73800472 Residual Norm: 0.00017941                                                 
               <<y;y>>_0.0656 :     -9.84674589 Residual Norm: 0.00018212                                                 
               <<z;z>>_0.0000 :     -9.34779283 Residual Norm: 0.00065565                                                 
               <<z;z>>_0.0656 :     -9.48860561 Residual Norm: 0.00064486                                                 
                                                                                                                          

* Info * Processing Fock builds... (batch size: 4)                                                                        
* Info *   batch 1/1                                                                                                      

                                                                                                                          
* Info * 14 gerade trial vectors in reduced space                                                                         
* Info * 24 ungerade trial vectors in reduced space                                                                       
                                                                                                                          
* Info * 126.65 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:   5 * Residuals (Max,Min): 5.88e-05 and 5.64e-06                                            
                                                                                                                          
               <<y;y>>_0.0000 :     -9.73800476 Residual Norm: 0.00000564                                                 
               <<y;y>>_0.0656 :     -9.84674594 Residual Norm: 0.00000865                                                 
               <<z;z>>_0.0000 :     -9.34779352 Residual Norm: 0.00003494                                                 
               <<z;z>>_0.0656 :     -9.48860634 Residual Norm: 0.00003797                                                 
                                                                                                                          

               *** Linear response converged in 5 iterations. Time: 1.93 sec

               Polarizability (w=0.0000)                                                                                  
               -------------------------                                                                                  
                                  X              Y              Z                                                         
               X         9.09692566     0.00000000    -0.00000000                                                         
               Y         0.00000000     9.73800476     0.00000000                                                         
               Z        -0.00000000     0.00000000     9.34779352                                                         
                                                                                                                          
               Polarizability (w=0.0656)                                                                                  
               -------------------------                                                                                  
                                  X              Y              Z                                                         
               X         9.29287577     0.00000000    -0.00000000                                                         
               Y         0.00000000     9.84674594    -0.00000000                                                         
               Z        -0.00000000     0.00000000     9.48860634                                                         
                                                                                                                          

Print the requested molecular properties.

print("Comp   Freq     alpha")
print("-"*21)
for key in lrf_results["response_functions"].keys():
    print(f" {key[0]}{key[1]}{key[2]:8.4f}{-lrf_results['response_functions'][key]:10.6f}")

Comp   Freq     alpha
---------------------
 xx  0.0000  9.096926
 yx  0.0000  0.000000
 zx  0.0000 -0.000000
 xx  0.0656  9.292876
 yx  0.0656  0.000000
 zx  0.0656 -0.000000
 xy  0.0000  0.000000
 yy  0.0000  9.738005
 zy  0.0000  0.000000
 xy  0.0656  0.000000
 yy  0.0656  9.846746
 zy  0.0656  0.000000
 xz  0.0000 -0.000000
 yz  0.0000  0.000000
 zz  0.0000  9.347794
 xz  0.0656 -0.000000
 yz  0.0656 -0.000000
 zz  0.0656  9.488606

Complex response function#

We calculate the needed linear response functions with use of the ComplexResponse class.

lrf = vlx.ComplexResponse() 

lrf.a_operator = "electric dipole"
lrf.b_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"

lrf.a_components = ["z"]
lrf.b_components = ["z"]

lrf.damping = 0.00454

lrf.frequencies = [0.0, 0.0656]

lrf_results = lrf.compute(molecule, basis, scf_results)