Alchemical calculations allow to obtain electron densities and total electron energies of a set of isoelectronic target molecules without evaluating the electron density of the target molecule explicitly in a self-consistent fashion. While the theoretical background is detailed in aforementioned reference, this is more about the practical approach how to perform those calculations numerically. We will focus on Gaussian as quantum chemical program even though Horton, Orca and GPAW support similar calculations.
At first, we evaluate the electronic energy of N2 with the following minimal input file:
%Chk=run.chk #CCSD 6-31G(d) scf=tight run 0 1 N N 1 1.1
This calculates the total energy of a N2 molecule with a bond length of 1.1 Angstrom which is written to the output log file as
CCSD=-109.2552722
Exercise: Calculate the dissociation profile of N2 and plot the electronic energy as function of the bond length for five lengths between 0.9 and 1.5 Angstrom.
For alchemical transformations, we need to connect the states smoothly via the coupling parameter lambda. This requires a superimposed basis set, since the target molecule would be calculated with the basis set of the reference molecule otherwise. This in turn leads to an imprecise Taylor expansion where the approximated target (target molecule at reference basis set) is different from the actual target (target molecule at target basis set). Superimposing the basis sets makes the calculations numerically less stable (exercise: why?) and more expensive (exercise: why?) but is a requirement of the finite difference scheme that is employed in this example. Superpositioning basis sets in Gaussian is done with the Gen keyword (the other ones are technically needed for completion of the run):
%Chk=run.chk #CCSD(Full,MaxCyc=100) Gen scf=tight Massage guess=indo integral=NoXCTest Pop=None Density=CC run 0 1 N N 1 [Bond length here] N 0 [paste 6-31G(3df,3pd) basis sets in Gaussian94 format from https://bse.pnl.gov/bse/portal here for all involved atoms] **** 1 Nuc [Nuclear charge of site 1] 2 Nuc [Nuclear charge of site 2]
Once the self consistent electron density has been found, we need to evaluate it for the target molecule. This requires a new run of Gaussian where the density is not found self-consistently, but rather read in from the previous calculation. The keyword Density=Checkpoint does exactly this. In the following example, we calculate the electron-nuclei interaction of BF for the N2 density.
%Chk=run.chk #CCSD(T) sp scf=tight Massage guess=indo integral=NoXCTest Density=Checkpoint run 0 1 F B 1 [Bond length here]
When inspecting the output, it is advised to check the total charge output which can be found in a line similar to this:
Charge= 0.0000 electrons
If the total charge is not zero, the basis set has been read incorrectly from the checkpoint file. Now the log file contains energies as follows:
N-N= 1.984414532213D+01 E-N=-3.142289395254D+02
Where N-N is the Coulomb interaction between the nuclei and E-N is the interaction between the nuclei and the electron density. The energy prediction of the target molecule is then given as
E(t) = E(r) - EN(r) + EN(t) - NN(r) + NN(t)
All relevant energies are given in atomic units, i.e. Hartrees.
Exercise: Calculate the bond potential of N2 by performing self-consistent calculations for CO only.