Alchemical perturbation

In this tutorial, we will show how to estimate properties of molecules using alchemical perturbation. Specifically, we will estimate the distance between two atoms in \(\mathrm{CO}\) and \(\mathrm{BF}\) molecules by only calculating the properties of \(\mathrm{N_2}\) and its alchemical perturbation.

Let’s denote the atomic number of the atoms in the diatomic molecule as

\[Z_{\pm}(\lambda) = 7 \pm \lambda\]

parameterized by a variable \(\lambda\). One of the atom takes the plus sign while another one takes the minus sign to keep the number of electrons constant.

The equilibrium distance between the atoms is defined as

\[s^*(\lambda) = \arg\min_s \mathcal{E}(s, \lambda)\]

where \(\mathcal{E}\) is the total energy as a function of the atomic distance \(s\) and \(\lambda\). What we will do is to estimate the equilibrium distance for \(\lambda = 1\) (for \(\mathrm{CO}\)) and \(\lambda = 2\) (for \(\mathrm{BF}\)) using Taylor expansion,

\[s^*(\lambda) \approx s^*(0) + \lambda \frac{\partial s^*}{\partial \lambda} + \frac{1}{2} \lambda^2 \frac{\partial^2 s^*}{\partial \lambda^2}\]

As a demonstration, we will use Hartree-Fock calculation with 3-21G basis set. First, we need to import modules and set up variables that we will need for the calculations.

import torch
import dqc
import xitorch.optimize  # for differentiable optimization
dtype = torch.double
basis = dqc.loadbasis("7:3-21G")

The 3-21G basis for atomz 7 does not exist, but we will download it
Downloaded to /home/docs/checkouts/readthedocs.org/user_builds/dqc/envs/latest/lib/python3.8/site-packages/dqc/api/.database/3-21g/07.gaussian94

xitorch is a great library that provides differentiable functionals that we will use in this tutorial. The last line with dqc.loadbasis() loads the basis 3-21G for atomic number 7. We will use the same basis for all values of \(\lambda\) to make sure there is no discontinuity in the properties.

Next, we need to define a function that calculates the energy given the distance \(s\) and \(\lambda\).

def get_energy(s, lmbda):
    atomzs = 7.0 + torch.tensor([1.0, -1.0], dtype=dtype) * lmbda
    atomposs = torch.tensor([[-0.5, 0, 0], [0.5, 0, 0]], dtype=dtype) * s
    mol = dqc.Mol((atomzs, atomposs), spin=0, basis=[basis, basis])
    qc = dqc.HF(mol).run()
    return qc.energy()

Once the function is defined, then we can calculate the equilibrium distance for \(\mathrm{N_2}\) molecule.

lmbda = torch.tensor(0.0, dtype=dtype).requires_grad_()
s0_n2 = torch.tensor(2.04, dtype=dtype)  # initial guess of the distance
smin_n2 = xitorch.optimize.minimize(
    get_energy, s0_n2, params=(lmbda,), method="gd", step=1e-2)
print(smin_n2)

tensor(2.0460, dtype=torch.float64, grad_fn=<_RootFinderBackward>)

xitorch.optimize.minimize finds the parameters s that minimizes the energy given the parameters lmbda. The output of xitorch.optimize.minimize is now differentiable with respect to the parameter lmbda.

grad_lmbda = torch.autograd.grad(smin_n2, lmbda, create_graph=True)[0]
grad2_lmbda = torch.autograd.grad(grad_lmbda, lmbda, create_graph=True)[0]
print(grad_lmbda.detach(), grad2_lmbda.detach())

tensor(1.3323e-15, dtype=torch.float64) tensor(0.1323, dtype=torch.float64)

Now, we can estimate the equilibrium distance of \(\mathrm{CO}\) and \(\mathrm{BF}\),

smin_co = smin_n2 + grad_lmbda + 0.5 * grad2_lmbda
smin_bf = smin_n2 + grad_lmbda * 2 + 0.5 * grad2_lmbda * 2 ** 2
print(smin_co.detach(), smin_bf.detach())

tensor(2.1121, dtype=torch.float64) tensor(2.3105, dtype=torch.float64)

For reference, the equilibrium distances for \(\mathrm{CO}\) and \(\mathrm{BF}\) by minimizing the energy are 2.1119 and 2.3103 Bohr, respectively, which are quite close to the values above.