Machine Learning and Quantum Alchemy: Day 3

Part of the course Machine Learning and Quantum Alchemy.

Make sure to install pyscf first. If you use Google Colab, you can enter !pip install pyscf.

Exercises

1. Use the following code to calculate the total energy of a dimer given a distance. Make sure to understand the basic steps in the code. Plot the bonding potential of CO and N\(_2\) using "STO-3G" as basis.

   import pyscf.scf
   import pyscf.gto
   import numpy as np
   
   def dimer_energy(element1, element2, distance, basis):
       mol = pyscf.gto.M(atom=f"{element1} 0 0 0; {element2} 0 0 {distance}", basis=basis, verbose=0)
       calc = pyscf.scf.RHF(mol)
       return calc.kernel()
   

   element1 and element2 are the element labels of the two atoms as string, e.g. "N" for a nitrogen. distance is the distance between the atom in Angstrom and basis defines the quality of the quantum chemistry calculation.

   Example Solution

   xs = np.linspace(0.8, 2, 10)
   ys = [dimer_energy("N", "N", _, "STO-3G") for _ in xs]
   plt.plot(xs, ys, label="N2")
   ys = [dimer_energy("C", "O", _, "STO-3G") for _ in xs]
   plt.plot(xs, ys, label="CO")
   plt.legend()
   

2. A general strategy to calculate the derivative of a function is to use finite differences. From the definition of the (partial) derivative we can obtain a discretized simple expression for an approximate derivative:

   \[\left.\frac{\partial f(x)}{\partial x}\right|_{x=x_0} := \lim_{\delta\rightarrow 0} \frac{f(x_0+\delta)-f(x_0)}{\delta}\Rightarrow \left.\frac{\partial f(x)}{\partial x}\right|_{x=x_0} \approx \frac{f(x_0+\delta)-f(x_0)}{\delta}\]

   Use the function dimer_energy from above to calculate the first derivative using finite differences at \(x_0=1\) for N\(_2\) and the basis "STO-3G". Use different displacements \(\delta\) and make a plot showing the derivative value as a function of \(\delta\). What do you observe? How small can you make \(\delta\)? For comparison, the exact value is -1.3979936186235125.

   Example Solution

   def get_derivative(delta, center):
       a = dimer_energy("N", "N", 1 + delta, "STO-3G")
       return (a - center) / delta
   center = dimer_energy("N", "N", 1, "STO-3G")
   deltas = 10**np.linspace(-20, 0)
   derivs = [get_derivative(delta, center) for delta in deltas]
   refvalue  = -1.3979936186235125
   plt.loglog(deltas, np.abs(np.array(derivs)-refvalue))
   

3. Use the following code to calculate alchemical changes in a dimer. Plot the alchemical change from N\(_2\) to CO and the reverse once for basis="STO-3G" and once for basis="unc-def2-TZVP". Choose any reasonable distance for your plots. What do you observe?

   def dimer_energy(element1, element2, distance, alchemical_change, basis):
       mol = pyscf.gto.M(atom=f"{element1} 0 0 0; {element2} 0 0 {distance}", basis=basis, verbose=0)
       calc = pyscf.scf.RHF(mol)
       h1 = calc.get_hcore()
   
       s = 0
       for i, delta_Z in enumerate([1, -1]):
           mol.set_rinv_orig_(mol.atom_coords()[i])
           s -= delta_Z * alchemical_change * mol.intor("int1e_rinv")
   
       nn = (alchemical_change * np.diff(mol.atom_charges())[0] - alchemical_change**2)
       nn /= (distance / pyscf.data.nist.BOHR)
   
       calc.get_hcore = lambda *args, **kwargs: h1 + s
       return calc.kernel() + nn
   

   element1 and element2 are the element labels of the two atoms as string, e.g. "N" for a nitrogen. distance is the distance between the atom in Angstrom and basis defines the quality of the quantum chemistry calculation. The now additional parameter alchemical_change modifies the nuclear charges of the involved atoms such that their sum remains the same but one is incremented by alchemical_change and one is decremented by the same amount.

   Example Solution

   xs = np.linspace(0, 1, 10)
   basis = "def2-TZVP"
   N2_to_CO = [dimer_energy("N", "N", 1, _, basis) for _ in xs]
   CO_to_N2 = [dimer_energy("C", "O", 1, _, basis) for _ in xs]
   plt.plot(xs, N2_to_CO)
   plt.plot(xs, CO_to_N2[::-1])
   

4. Use first order finite differences to estimate the energy of N\(_2\) coming from CO. Why is it inaccurate?

5. For second order finite differences, one expression is

   \[\left.\frac{\partial^2 f(x)}{\partial x^2}\right|_{x=x_0} \approx \frac{f(x_0+\delta)-2f(x_0)+f(x_0-\delta)}{\delta^2}\]

   Use this to predict make a second-order Taylor approximation of the energy of N\(_2\) coming from CO:

   \[f(x+\delta)\approx f(x_0) + \delta f'(x_0) + \frac{\delta^2}{2}f''(x_0)\]

   Is that better now and why?