Example 2: PennyLane VQE calculation on water

Unlike the previous case, this example showcases the use of qc2 in combination with PennyLane to perform a oo-VQE run on water molecule. All circuit evaluations are conducted using the default.qubit state simulator device, which provides exact expectation values.

Import required packages

[1]:

import subprocess

# ASE molecule object
from ase.build import molecule

# PennyLane-related packages
import pennylane as qml
from pennylane import numpy as np

# ignore package import warnings
import warnings
warnings.filterwarnings('ignore')

# qc2 packages
from qc2.data import qc2Data
from qc2.ase import Psi4
from qc2.algorithms.utils import ActiveSpace
from qc2.algorithms.pennylane import oo_VQE

Instantiate `qc2Data` class and run qc2-ASE calculator

[2]:

# set Atoms object
mol = molecule('H2O')

# instantiate qc2Data class
qc2data = qc2Data(
    molecule=mol,
    filename='h2o.hdf5',
    schema='qcschema'
)

# specify the qchem calculator
qc2data.molecule.calc = Psi4(method='hf', basis='sto-3g')

# run calculation and save qchem data in the fcidump file
qc2data.run()

  Python driver attempted to set threads to 1.
  Psi4 was compiled without OpenMP, setting threads to 1.
  Python driver attempted to set threads to 1.
  Psi4 was compiled without OpenMP, setting threads to 1.
* Reference energy (Hartree): -74.96449224627995
* Saving qchem data in h2o.hdf5

This will generate a h2o.hdf5 file containing all relevant qchem data according to the QCSchema. Like before, let’s check its data structure via the shell command:

[3]:

!h5dump -n h2o.hdf5

HDF5 "h2o.hdf5" {
FILE_CONTENTS {
 group      /
 group      /keywords
 group      /model
 group      /molecule
 group      /properties
 group      /provenance
 group      /wavefunction
 dataset    /wavefunction/localized_orbitals_a
 dataset    /wavefunction/localized_orbitals_b
 dataset    /wavefunction/scf_eigenvalues_a
 dataset    /wavefunction/scf_eigenvalues_b
 dataset    /wavefunction/scf_eri
 dataset    /wavefunction/scf_eri_mo_aa
 dataset    /wavefunction/scf_eri_mo_ab
 dataset    /wavefunction/scf_eri_mo_ba
 dataset    /wavefunction/scf_eri_mo_bb
 dataset    /wavefunction/scf_fock_a
 dataset    /wavefunction/scf_fock_mo_a
 dataset    /wavefunction/scf_fock_mo_b
 dataset    /wavefunction/scf_orbitals_a
 dataset    /wavefunction/scf_orbitals_b
 }
}

You could also use !h5dump h2o.hdf5 to inspect the numerical data contained within each group and dataset.

Instantiate `qc2.algorithms.oo_VQE` class

[4]:

# set up activate space
active_space = ActiveSpace(
    num_active_electrons=(2, 2),  # => (n_alpha, n_beta)
    num_active_spatial_orbitals=4 # => active orbitals
)

# instantiate oo_VQE class
qc2data.algorithm = oo_VQE(
    active_space=active_space,
    mapper="jw",                                           # => use Jordan-Wigner mapper
    optimizer=qml.GradientDescentOptimizer(stepsize=0.5),  # => GradientDescentOptimizer optimizer from PennyLane
    device="default.qubit",
    max_iterations=100
)

Run oo-VQE

[5]:

# add extra options to `device` and `QNode` if needed...
results = qc2data.algorithm.run(
    device_kwargs={"shots": None},
    qnode_kwargs={"diff_method": "best"}
)

>>> Optimizing circuit and orbital parameters...
iter = 000, energy = -74.964404821247 Ha
iter = 001, energy = -74.974654732364 Ha
iter = 002, energy = -74.977136100996 Ha
iter = 003, energy = -74.977654273872 Ha
iter = 004, energy = -74.977797944648 Ha
iter = 005, energy = -74.977844599529 Ha
iter = 006, energy = -74.977860462172 Ha
iter = 007, energy = -74.977866022280 Ha
iter = 008, energy = -74.977867940065 Ha
iter = 009, energy = -74.977868636996 Ha
iter = 010, energy = -74.977868887832 Ha
iter = 011, energy = -74.977868978674 Ha
optimization finished.

=== PENNYLANE oo-VQE RESULTS ===
* Total ground state energy (Hartree): -74.977868978674

[6]:

# print optimized circuit parameters
print(f'* Optimized circuit parameters:')
print(results.optimal_circuit_params, "\n")

# orbital parameters
print(f'* Optimized orbital parameters:')
print(results.optimal_orbital_params, "\n")

print(f'* oo-VQE energy (Hartree):')
print(results.optimal_energy, "\n")

# vectors containing all intermediate results
# results.optimizer_evals, results.energy, results.circuit_parameters, results.orbital_parameters

* Optimized circuit parameters:
[ 2.07765036e-02 -1.11833627e-16 -2.07765036e-02 -1.57494321e-17
 -4.82754622e-17 -8.98925217e-17 -6.71880313e-17 -3.46826285e-17
  1.56019786e-01 -5.54938723e-17 -2.59957419e-17  4.35653437e-02
  1.45114627e-18 -9.95617852e-18  1.18705568e-17  1.02981524e-17
  1.21106275e-17  6.38136431e-18  9.59524869e-18  9.84554818e-18
  3.12784437e-19  3.55066164e-18  5.68534374e-02  3.11924957e-18
 -1.73185321e-19  3.13803385e-02]

* Optimized orbital parameters:
[0.0026712593299834616, -0.6190212888834982, 6.275542869417149e-10, -1.113516977506888e-16, 3.8213383832866345e-15, 5.586006655411282e-14, 7.555308746399609e-15, -2.7210673457436863e-05, -0.0019899144798891403, 9.447639902049476e-12, 0.0019118957526852315, 2.3253980286437208e-15, 2.0963033505657805e-15, 1.8613420600403597e-12, -0.002273727905803901, 2.9862067670493623e-12, -6.067518761511912e-16, 5.340727372484386e-10]

* oo-VQE energy (Hartree):
-74.97786897867395

Compare oo-VQE result with classical qchem calculations

Once again, let’s compare our oo-VQE energy with the one obtained from Psi4 calculations alone:

[7]:

import psi4

# Set Psi4 to run in serial mode
psi4.set_num_threads(1)

# Define the molecule in XYZ format
h2o = psi4.geometry("""
O 0.0000  0.000000000  0.1192622641
H 0.0000  0.763237638 -0.4770469398
H 0.0000 -0.763237638 -0.4770469398
symmetry c1
noreorient
nocom
""")

# Setup CAS space (Active space definition)
n_active_elec = 4    # Number of active electrons
n_active_orb = 4     # Number of active orbitals

# Set computation options
psi4.set_options({
    'basis': 'sto-3g',
    'reference': 'rhf'
})

# run reference Hartree-Fock
e_scf, wfn_scf = psi4.energy('scf', return_wfn=True)

# Perform FCI calculation
e_fci, wfn_fci = psi4.energy('fci/{}'.format(n_active_orb, n_active_elec), ref_wfn=wfn_scf, return_wfn=True)

# Set options for CASSCF
psi4.set_options({
    'icore': 1,
    'restricted_docc': 3
})

# Perform CASSCF calculation
e_casscf, wfn_casscf = psi4.energy('casscf/{}'.format(n_active_orb, n_active_elec), ref_wfn=wfn_scf, return_wfn=True)

# Print results
print('')
print("* Final oo-VQE energy:        {:.6f} hartree".format(results.optimal_energy))
print("* Final CASSCF/sto-3g energy: {:.6f} hartree".format(e_casscf))
print("* Final FCI/sto-3g energy:    {:.6f} hartree".format(e_fci))


* Final oo-VQE energy:        -74.977869 hartree
* Final CASSCF/sto-3g energy: -74.977870 hartree
* Final FCI/sto-3g energy:    -75.015429 hartree

As seen, our final oo-VQE energy is exactly the one expected when performing a CASSCF/sto-3g single-point calculation with the chosen active space.