Exercise 8

Part of the course Computational Chemistry.

Gradients are clearly relevant for many applications: e.g. to propagate a system using molecular dynamics or to minimize a geometry. Calculating them is a standard capability for quantum chemistry packages.

Note: these settings are not reliable and are only used for demonstration purposes. Please do not use them in actual research unless you know what you are doing.

Task 8.1: Slippery Slope

This time, our example system will be carbon monoxide. Write a function get_energy(distance: float) -> float to calculate the total energy of the system using restricted Hartree-Fock and the STO-3G basis set. Write another function get_force(distance: float) -> float which returns the deriative of the energy w.r.t. the distance of the two atoms using RHF.nuc_grad_method(). Finally, write a function get_finite_difference_force(distance: float) -> float which achieves the same using finite differences. Do the results agree?

Task 8.2: Time well spent

Extend the functions get_energy() and get_force() to accept a basis set get_force(distance: float, basis_set: str) -> float. Using time.time() measure the time it takes to evaluate just the energy and just the force for the basis sets cc-pVDZ, cc-pVTZ, cc-pVQZ. Make sure to average over a few runs to obtain reliable timings. What do you observe?

Task 8.3: Hellmann-Feynman

Hellmann-Feynman forces are obtained using the Hellmann-Feynman theorem which (subject to assumptions) states

\[\frac{\partial E}{\partial \lambda} = \left\langle \psi\left|\frac{\partial \hat{H}}{\partial \lambda}\right|\psi\right\rangle\]

For the electrostatic force acting on the atoms, this results in

\[\mathbf{F}_I = Z_I\int d\mathbf {r'}\rho(\mathbf{r'})\frac{(\mathbf{r'}-\mathbf{R'})}{|\mathbf{r'}-\mathbf{R}|^3}\]

We can use numerical integration. In PySCF, one way to achieve this is to use a DFT integration grid:

grid = pyscf.dft.gen_grid.Grids(mol)
grid.level = 3
grid.build()
dm1_ao = calc.make_rdm1()
ao_value = pyscf.dft.numint.eval_ao(mol, grid.coords, deriv=0)
rho = pyscf.dft.numint.eval_rho(mol, ao_value, dm1_ao, xctype="LDA")

which yields the electron density at all grid points. Note that the grid points carry different integration weights grid.weights. First, make sure that the integral in space of all density is the number of electrons, then implement Hellmann-Feynman forces for CO. Compare them to the actual forces for the basis sets STO-3G and cc-pVQZ. What do you observe?