Finite field method

Finite field method#

The finite field method denotes the calculation of molecular properties of a given order by means of numerical differentiation of lower-order molecular properties. It is typically applied to the determination of electric properties since it is trivial to incorporate uniform static electric fields in the molecular Hamiltonian and carry out calculations in the presence of these fields.

Compared to analytic derivative techniques, finite field calculations are slower and less accurate but their appeal lies in the ease of implementation.

Energy expansion#

In the presence of a uniform static electric field, the molecular energy becomes equal to

\[ E(F) = E_0 - \mu_{\alpha} F_\alpha - \frac{1}{2} \alpha_{\alpha\beta} F_\alpha F_\beta - \frac{1}{6} \beta_{\alpha\beta\gamma} F_\alpha F_\beta F_\gamma - \frac{1}{24} \gamma_{\alpha\beta\gamma\delta} F_\alpha F_\beta F_\gamma F_\delta - \cdots \]

where the molecular property tensors \(\mu\), \(\alpha\), \(\beta\), and \(\gamma\) are the permanent dipole moment, polarizabilty, first-order hyperpolarizability, and second-order hyperpolarizability, respectively. These molecular properties refer to the isolated system and can thus be identified from derivatives of the energy in the limit of zero field strength:

\[\begin{align*} \mu_{\alpha} & = - \left. \frac{\partial E(F)}{\partial F_\alpha} \right|_{F = 0} \\ \alpha_{\alpha\beta}(0;0) & = - \left. \frac{\partial^2 E(F)}{\partial F_\alpha \partial F_\beta} \right|_{F = 0} \\ \beta_{\alpha\beta\gamma}(0;0,0) & = - \left. \frac{\partial^3 E(F)}{\partial F_\alpha \partial F_\beta \partial F_\gamma} \right|_{F = 0} \\ \gamma_{\alpha\beta\gamma\delta}(0;0,0,0) & = - \left. \frac{\partial^4 E(F)}{\partial F_\alpha \partial F_\beta \partial F_\gamma \partial F_\delta} \right|_{F = 0} \end{align*}\]

Dipole moment expansion#

It is also possible to define a dipole moment of the system in the presence of the external field. This field-dependent dipole moment takes the form

\[ \mu_\alpha^F = \mu_\alpha + \alpha_{\alpha\beta} F_\beta + \frac{1}{2} \beta_{\alpha\beta\gamma} F_\beta F_\gamma + \frac{1}{6} \gamma_{\alpha\beta\gamma\delta} F_\beta F_\gamma F_\delta - \cdots \]

We then get

\[\begin{align*} \alpha_{\alpha\beta}(0;0) & = \left. \frac{\partial \mu_\alpha^F}{\partial F_\beta } \right|_{F = 0} \\ \beta_{\alpha\beta\gamma}(0;0,0) & = \left. \frac{\partial^2 \mu_\alpha^F}{\partial F_\beta \partial F_\gamma} \right|_{F = 0} \\ \gamma_{\alpha\beta\gamma\delta}(0;0,0,0) & = \left. \frac{\partial^3 \mu_\alpha^F}{\partial F_\beta \partial F_\gamma \partial F_\delta} \right|_{F = 0} \end{align*}\]

Polarizability expansion#

It is also possible to define a polarizability of the system in the presence of the external field. This field-dependent polarizability takes the form

\[ \alpha_{\alpha\beta}^F(-\omega;\omega) = \alpha_{\alpha\beta}(-\omega;\omega) + \beta_{\alpha\beta\gamma}(-\omega;\omega,0) F_\gamma + \frac{1}{2} \gamma_{\alpha\beta\gamma\delta}(-\omega;\omega,0,0) F_\gamma F_\delta - \cdots \]

We then get

\[\begin{align*} \beta_{\alpha\beta\gamma}(-\omega;\omega,0) & = \left. \frac{\partial \alpha_{\alpha\beta}^F(-\omega;\omega)}{\partial F_\gamma} \right|_{F = 0} \\ \gamma_{\alpha\beta\gamma\delta}(-\omega;\omega,0,0) & = \left. \frac{\partial^2 \alpha_{\alpha\beta}^F(-\omega;\omega)}{\partial F_\gamma \partial F_\delta} \right|_{F = 0} \end{align*}\]

Hyperpolarizability expansion#

It is also possible to define a hyperpolarizability of the system in the presence of the external field. This field-dependent hyperpolarizability takes the form

\[ \beta_{\alpha\beta\gamma}^F(-\omega_\sigma;\omega_1, \omega_2) = \beta_{\alpha\beta\gamma}(-\omega_\sigma;\omega_1,\omega_2) + \gamma_{\alpha\beta\gamma\delta}(-\omega_\sigma;\omega_1,\omega_2,0) F_\delta - \cdots \]

We then get

\[\begin{align*} \gamma_{\alpha\beta\gamma\delta}(-\omega_\sigma;\omega_1,\omega_2,0) & = \left. \frac{\partial \beta_{\alpha\beta\gamma}^F(-\omega_\sigma;\omega_1,\omega_2)}{\partial F_\delta} \right|_{F = 0} \end{align*}\]

where \(\omega_\sigma = \omega_1+ \omega_2\).

Numerical example#

Let us illustrate the finite field method with practical example.

import numpy as np
import veloxchem as vlx

We define a water molecule with the dipole moment aligned along the positive z-axis.

h2o_xyz = """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(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "def2-svpd", ostream=None)

Setting up drivers for:

  • SCF optimization

  • first-order properties (for now an explicit function, will be changed later)

  • linear response function (real)

  • quadratic response function (complex)

  • cubic response function (complex)

Hide code cell source 
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()

scf_drv.xcfun = "b3lyp"
Hide code cell source 
element_charges = {"H": 1.0, "O": 8.0}

# electric-dipole
dipole_drv = vlx.ElectricDipoleIntegralsDriver()

dipole_mats = dipole_drv.compute(molecule, basis)

mu_x = -1.0 * dipole_mats.x_to_numpy()
mu_y = -1.0 * dipole_mats.y_to_numpy()
mu_z = -1.0 * dipole_mats.z_to_numpy()


def dipmom(D):
    mu_e = np.zeros(3)
    mu_n = np.zeros(3)

    # electronic part
    mu_e[0] = np.einsum("ab, ab", D, mu_x)
    mu_e[1] = np.einsum("ab, ab", D, mu_y)
    mu_e[2] = np.einsum("ab, ab", D, mu_z)

    # nuclear part
    for A, molecule_atom_label in enumerate(molecule.get_labels()):
        R_A = np.array(molecule.get_atom_coordinates(A))
        Z_A = element_charges[molecule_atom_label]

        mu_n += Z_A * R_A

    # total dipole moment
    mu = mu_e + mu_n

    return mu[2]
Hide code cell source 
lrf_drv = vlx.LinearResponseSolver()

lrf_drv.ostream.mute()

lrf_drv.a_operator = "electric dipole"
lrf_drv.b_operator = "electric dipole"

lrf_drv.a_components = "z"
lrf_drv.b_components = "z"

lrf_drv.frequencies = [0.0]
Hide code cell source 
qrf_drv = vlx.QuadraticResponseDriver()

qrf_drv.ostream.mute()

qrf_drv.a_component = "z"
qrf_drv.b_component = "z"
qrf_drv.c_component = "z"

qrf_drv.b_frequencies = [0.0]
qrf_drv.c_frequencies = [0.0]

qrf_drv.damping = 0.0
Hide code cell source 
crf_drv = vlx.CubicResponseDriver()

crf_drv.ostream.mute()

crf_drv.a_component = "z"
crf_drv.b_component = "z"
crf_drv.c_component = "z"
crf_drv.d_component = "z"

crf_drv.b_frequencies = [0.0]
crf_drv.c_frequencies = [0.0]
crf_drv.d_frequencies = [0.0]

crf_drv.damping = 0.0

Numerical differentiation#

Up to fourth-order numerical derivatives are calculated with use of a one-dimensional five-point stencil.

def numerical_derivatives(f):
    # numerical differentiation based on the five-point stencil

    d1 = (-f[2 * h] + 8 * f[h] - 8 * f[-h] + f[-2 * h]) / (12 * h)
    d2 = (-f[2 * h] + 16 * f[h] - 30 * f[0.0] + 16 * f[-h] - f[-2 * h]) / (12 * h**2)
    d3 = (f[2 * h] - 2 * f[h] + 2 * f[-h] - f[-2 * h]) / (2 * h**3)
    d4 = (f[2 * h] - 4 * f[h] + 6 * f[0.0] - 4 * f[-h] + f[-2 * h]) / h**4

    return (d1, d2, d3, d4)

Field-dependent molecular properties#

The following static molecular properties are determined as functions of the strength, \(F\), of the electric field applied along the positive z-axis:

  • electronic energy, \(E(F)\)

  • dipole moment, \(\mu_z^F\)

  • polarizability, \(\alpha_{zz}^F(0;0)\)

  • first-order hyperpolarizaiblity, \(\beta_{zzz}^F(0;0,0)\)

  • second-order hyperpolarizaiblity, \(\gamma_{zzzz}^F(0;0,0,0)\)

E = {}
mu = {}
alpha = {}
beta = {}
gamma = {}

h = 0.001
field_strengths = np.linspace(-2, 2, 5) * h # five-point stencil

for F in field_strengths:
    scf_drv.electric_field = [0.0, 0.0, F]
    scf_results = scf_drv.compute(molecule, basis)

    E[F] = scf_results["scf_energy"]

    mu[F] = dipmom(scf_results["D_alpha"] + scf_results["D_beta"])

    lrf_results = lrf_drv.compute(molecule, basis, scf_results)
    alpha[F] = -lrf_results["response_functions"][("z", "z", 0.0)]

    qrf_results = qrf_drv.compute(molecule, basis, scf_results)
    beta[F] = -qrf_results[('qrf', 0.0, 0.0)].real

    crf_results = crf_drv.compute(molecule, basis, scf_results)
    gamma[F] = -crf_results[('crf', 0.0, 0.0, 0.0)].real

Calculate derivatives#

dE = numerical_derivatives(E)
dmu = numerical_derivatives(mu)
dalpha = numerical_derivatives(alpha)
dbeta = numerical_derivatives(beta)

Results#

print("-"*50)
print(f"{'Molecular properties':>40s}")
print(" "*10 + "-"*40)
print(f"{'Method':10s}{'mu':>10s}{'alpha':>10s}{'beta':>10s}{'gamma':>10s}")
print("-"*50)
print("Analytic derivatives")
print(f"{'':10s}{mu[0.0]:10.6f}{alpha[0.0]:10.5f}{beta[0.0]:10.5f}{gamma[0.0]:10.3f}")
print("\nNumerical differentiation")
print(f"{' -energy':10s}{-dE[0]:10.6f}{-dE[1]:10.5f}{-dE[2]:10.5f}{-dE[3]:10.3f}")
print(f"{' -mu':10s}{'':10s}{dmu[0]:10.5f}{dmu[1]:10.5f}{dmu[2]:10.3f}")
print(f"{' -alpha':10s}{'':10s}{'':10s}{dalpha[0]:10.5f}{dalpha[1]:10.3f}")
print(f"{' -beta':10s}{'':10s}{'':10s}{'':10s}{dbeta[0]:10.3f}")
print("-"*50)

--------------------------------------------------
                    Molecular properties
          ----------------------------------------
Method            mu     alpha      beta     gamma
--------------------------------------------------
Analytic derivatives
            0.744762   9.34779  -6.73492   999.154

Numerical differentiation
 -energy    0.744763   9.34779  -6.71206  1047.880
 -mu                   9.34781  -6.73532   999.819
 -alpha                         -6.73499   999.151
 -beta                                     999.159
--------------------------------------------------

Using the results obtained with analytic derivatives as reference, we note that up to second-order derivatives are quite stable with the chosen thresholds of convergence (in the SCF optimization and response equation solvers) and strength of the applied field. Fourth-order derivatives are not reliable as seen in the result for \(\gamma_{zzzz}(0;0,0,0)\) based on differentiation of the energy.

Contents