The module
run_psi4.pyprovides a user-friendly way of running Psi4 calculations in OpenFermion. The basic idea is that once one generates aMolecularDatainstance, once can then call Psi4 with a specification of certain options (for instance, how muhc memory to use and what calculations to do_ in order to compute things about the molecule, update theMolecularDataobject, and save results of the calculation.The most common calculations users will want in OpenFermion would probabl be self-consistent field (aka Hartree-Fock calulations). Our code uses these "SCF" calculations to compute orbitals integrals, Hartree-Fock, energy, and more. Other common calculations are CISD and FCI calculations which also compute the 1-RDM and 2-RDM associated with their answers; CCSD calculations which also compute the CCSD amplitudes (useful for UCC) and MP2 calculations which can also be used to UCCSD initial guesses. Note that the
"delete_input"and"delete_output"options indicate whether to save automatically generated Psi4 input files or not.To use this plugin, you will need to personally download Psi4. Python code in
psi4_templateinstructs Psi4 (not python) to load theMolecularDataobject, populate it with information from calculations, and then save then save theMolecularDataobject as an HDF5. The modulerun_psi4uses subprocesses to actuallyrun_psi4and then loads the pickledMolecularData. Let's look at a simple example where we compute the energy of H$_2$ at various bond lengths.Warnings: electronic sructure calculations are finicky. They sometimes fail for surprising reasons. If a particular calculation is not converging it is probably necessary to change some of the SCF options set in
psi4_template. See the Psi4 documentation for more information or consult an electronic structure theory expert.
In [1]:
%matplotlib inline
In [2]:
from openfermion.hamiltonians import MolecularData
from openfermionpsi4 import run_psi4
# Set molecule parameters
basis = 'sto-3g'
multiplicity = 1
bond_length_interval = 0.2
n_points = 10
# Set calculation parameters
run_scf = 1
run_mp2 = 1
run_cisd = 0
run_ccsd = 0
run_fci = 1
delete_input = True
delete_output = True
# Generate molecule at different bond lengths
hf_energies = []
fci_energies = []
bond_lengths = []
for point in range(1, n_points + 1):
bond_length = bond_length_interval * float(point)
bond_lengths += [bond_length]
geometry = [('H', (0.,0.,0.)), ('H', (0.,0.,bond_length))]
molecule = MolecularData(
geometry, basis, multiplicity,
description=str(round(bond_length, 2)))
# Run Psi4
molecule = run_psi4(molecule,
run_scf=run_scf,
run_mp2=run_mp2,
run_cisd=run_cisd,
run_ccsd=run_ccsd,
run_fci=run_fci)
# Print out some results of calculation
print(f'\nAt bond length of {bond_length} ångstrom, molecular hydrogen has:')
print(f'Hartree-Fock energy of {molecule.hf_energy} Hartree.')
print(f'MP2 energy of {molecule.mp2_energy} Hartree.')
print(f'FCI energy of {molecule.fci_energy}')
print(f'Nuclear repulsion energy between protons is {molecule.nuclear_repulsion}')
for orbital in range(molecule.n_orbitals):
print(f'Spatial orbital {orbital} has energy of {molecule.orbital_energies[orbital]}')
hf_energies += [molecule.hf_energy]
fci_energies += [molecule.fci_energy]
# Plot
import matplotlib.pyplot as plt
plt.figure(0)
plt.plot(bond_lengths, fci_energies, 'x-'); plt.plot(bond_lengths, hf_energies, 'o-');
plt.ylabel('Energy in Hartree'); plt.xlabel('Bond length in ångstrom')
Out[2]:
I don't know what I did but this is cool. 🤓