Example 3: Computation of Energies with diabatization

In this example it is shown, how to run the basic energy computations with novel approach for diabatization of orbitals.

First, we specify the geometry of our system.

symbols = ['C', 'N', 'H', 'H', 'H']
coords = [[-0.000000000, 0.000000000, -0.000000000],
            [-0.000000000, 0.000000000, 1.412052000],
            [-0.000000000, -0.944215264, -0.525349226],
            [-0.000000000, 0.944215264, -0.525349226],
            [-0.176643782, 0.000000000, 2.413848671]]

Now, let’s define our system’s properties.

n_orbs_active = 2
n_elec_active = 2
charge = 0
multiplicity = 1
basis = 'sto-3g'

For the diabatization process, we use previously obtained wavefunction. Let’s load it now.

import psi4

load_wfn = psi4.core.Wavefunction.from_file(f'../tests/bruno_wfn_phi0_theta_180_sto3g')

Now, we create our problem using our previously computed wavefunction,

import saoovqe

problem = saoovqe.ProblemSet.with_dia_orbs_from_prev_wfn(symbols, coords, charge, multiplicity, n_elec_active,
                                                                 n_orbs_active, load_wfn, basis)

create initial circuits and ansatz,

n_states = 2
initial_circuits = saoovqe.OrthogonalCircuitSet.from_problem_set(n_states, problem)

repetitions = 1
ansatz = saoovqe.Ansatz.from_problem_set(saoovqe.AnsatzType.GUCCSD,
                                         problem,
                                         repetitions,
                                         qubit_mapper=problem.fermionic_mapper)

and finally, we create our solver.

from qiskit.primitives import Estimator

estimator = Estimator()
saoovqe_solver = saoovqe.SAOOVQE(estimator,
                                 initial_circuits,
                                 ansatz,
                                 problem,
                                 orbital_optimization_settings={})

As we now have all that is required, we calculate our energies.

from qiskit_algorithms.optimizers import SciPyOptimizer

energies = saoovqe_solver.get_energy(SciPyOptimizer('SLSQP', options={'maxiter': 500, 'ftol': 1e-8}))
energies

And now, we can compare them, with reference data:

import numpy as np

np.allclose(np.array([-92.69532984620506, -92.64558680739152]), energies)

And we can see, that the result is correct!