Obtaining energy level estimates

Estimating the values of energy levels is one of the principal applications of quantum chemistry. This article outlines how you can perform this for the canonical example of molecular hydrogen. The sample referenced in this section is MolecularHydrogen in the chemistry samples repository. A more visual example that plots the output is the MolecularHydrogenGUI demo.

Estimating the energy values of molecular hydrogen

The first step is to construct the Hamiltonian representing molecular hydrogen. Although you can construct this using the NWChem tool, for brevity, this sample adds the Hamiltonian terms manually.

// The code snippets in this section require the following namespaces.
// Make sure to include these at the top of your file or namespace.
using Microsoft.Quantum.Chemistry.OrbitalIntegrals;
using Microsoft.Quantum.Chemistry.Fermion;
using Microsoft.Quantum.Chemistry.Paulis;
using Microsoft.Quantum.Chemistry.QSharpFormat;
using Microsoft.Quantum.Simulation.Simulators;

    // These orbital integrals are represented using the OrbitalIntegral
    // data structure.
    var energyOffset = 0.713776188; // This is the coulomb repulsion
    var nElectrons = 2; // Molecular hydrogen has two electrons
    var orbitalIntegrals = new OrbitalIntegral[]
    {
        new OrbitalIntegral(new[] { 0,0 }, -1.252477495),
        new OrbitalIntegral(new[] { 1,1 }, -0.475934275),
        new OrbitalIntegral(new[] { 0,0,0,0 }, 0.674493166),
        new OrbitalIntegral(new[] { 0,1,0,1 }, 0.181287518),
        new OrbitalIntegral(new[] { 0,1,1,0 }, 0.663472101),
        new OrbitalIntegral(new[] { 1,1,1,1 }, 0.697398010),
        // This line adds the identity term.
        new OrbitalIntegral(new int[] { }, energyOffset)
    };

    // Initialize a fermion Hamiltonian data structure and add terms to it.
    var fermionHamiltonian = new OrbitalIntegralHamiltonian(orbitalIntegrals).ToFermionHamiltonian();

Simulating the Hamiltonian requires converting the fermion operators to qubit operators. This conversion is performed through the Jordan-Wigner encoding as follows:

    // The Jordan-Wigner encoding converts the fermion Hamiltonian, 
    // expressed in terms of fermionic operators, to a qubit Hamiltonian,
    // expressed in terms of Pauli matrices. This is an essential step
    // for simulating our constructed Hamiltonians on a qubit quantum
    // computer.
    var jordanWignerEncoding = fermionHamiltonian.ToPauliHamiltonian(QubitEncoding.JordanWigner);

    // You also need to create an input quantum state to this Hamiltonian.
    // Use the Hartree-Fock state.
    var fermionWavefunction = fermionHamiltonian.CreateHartreeFockState(nElectrons);

    // This Jordan-Wigner data structure also contains a representation 
    // of the Hamiltonian and wavefunction made for consumption by the Q# operations.
    var qSharpHamiltonianData = jordanWignerEncoding.ToQSharpFormat();
    var qSharpWavefunctionData = fermionWavefunction.ToQSharpFormat();
    var qSharpData = Convert.ToQSharpFormat(qSharpHamiltonianData, qSharpWavefunctionData);

Next, pass qSharpData, which represents the Hamiltonian, to the TrotterStepOracle function (available in the Microsoft.Quantum.Chemistry.JordanWigner namespace). TrotterStepOracle returns a quantum operation that approximates the real-time evolution of the Hamiltonian. For more information, see Simulating Hamiltonian dynamics.

operation RunTrotterStep(qSharpData: JordanWignerEncodingData) : Unit {
    // Choose the integrator step size
    let stepSize = 1.0;

    // Choose the order of the Trotter—Suzuki integrator.
    let integratorOrder = 4;

    // `oracle` is an operation that applies a single time-step of evolution for duration `stepSize`.
    // `rescale` is just `1.0/stepSize` -- the number of steps required to simulate unit-time evolution.
    // `nQubits` is the number of qubits that must be allocated to run the `oracle` operation.
    let (nQubits, (rescale, oracle)) =  TrotterStepOracle (qSharpData, stepSize, integratorOrder);
}

At this point, you can use the standard library's phase estimation algorithms to learn the ground state energy using the previous simulation. This requires preparing a good approximation to the quantum ground state. Suggestions for such approximations are provided in the Broombridge schema. However, absent these suggestions, the default approach adds a number of hamiltonian.NElectrons electrons to greedily minimize the diagonal one-electron term energies. The phase estimation functions and operations are provided in DocFX notation in the Microsoft.Quantum.Characterization namespace, as well as the Microsoft.Quantum.Simulation namespace. Make sure to open these in your Q# file when using the GetEnergyByTrotterization code shown here.

The following snippet shows how the real-time evolution output by the chemistry simulation library integrates with quantum phase estimation.

open Microsoft.Quantum.Intrinsic;
open Microsoft.Quantum.Chemistry.JordanWigner;
open Microsoft.Quantum.Characterization;
open Microsoft.Quantum.Simulation;

operation GetEnergyByTrotterization (
    qSharpData : JordanWignerEncodingData,
    nBitsPrecision : Int,
    trotterStepSize : Double,
    trotterOrder : Int) : (Double, Double) {

    // The data describing the Hamiltonian for all these steps is contained in
    // `qSharpData`
    let (nSpinOrbitals, fermionTermData, statePrepData, energyOffset) = qSharpData!;

    // Using a Product formula, also known as `Trotterization`, to
    // simulate the Hamiltonian.
    let (nQubits, (rescaleFactor, oracle)) = 
        TrotterStepOracle(qSharpData, trotterStepSize, trotterOrder);

    // The operation that creates the trial state is defined here.
    // By default, greedy filling of spin-orbitals is used.
    let statePrep = PrepareTrialState(statePrepData, _);

    // Using the Robust Phase Estimation algorithm
    // of Kimmel, Low and Yoder.
    let phaseEstAlgorithm = RobustPhaseEstimation(nBitsPrecision, _, _);

    // This runs the quantum algorithm and returns a phase estimate.
    let estPhase = EstimateEnergy(nQubits, statePrep, oracle, phaseEstAlgorithm);

    // Now, obtain the energy estimate by rescaling the phase estimate
    // with the trotterStepSize. We also add the constant energy offset
    // to the estimated energy.
    let estEnergy = estPhase * rescaleFactor + energyOffset;

    // Return both the estimated phase and the estimated energy.
    return (estPhase, estEnergy);
}

You can now invoke the Q# code from the host program. The following C# code creates a full-state simulator and runs GetEnergyByTrotterization to obtain the ground state energy.

    using (var qsim = new QuantumSimulator())
    {
        // Specify the bits of precision desired in the phase estimation algorithm
        var bits = 7;

        // Specify the step size of the simulated time evolution. The step size needs to
        // be small enough to avoid aliasing of phases, and also to control the
        // error of simulation.
        var trotterStep = 0.4;

        // Choose the Trotter integrator order
        Int64 trotterOrder = 1;

        // As the quantum algorithm is probabilistic, run a few trials.

        // This may be compared to true value of
        Console.WriteLine("Exact molecular hydrogen ground state energy: -1.137260278.\n");
        Console.WriteLine("----- Performing quantum energy estimation by Trotter simulation algorithm");
        for (int i = 0; i < 5; i++)
        {
            // EstimateEnergyByTrotterization
            var (phaseEst, energyEst) = GetEnergyByTrotterization.Run(qsim, qSharpData, bits, trotterStep, trotterOrder).Result;
            Console.WriteLine("Estimated molecular hydrogen ground state energy: {0}", energyEst);
        }
    }

The operation returns two parameters:

• energyEst is the estimate of the ground state energy and should be close to -1.137 on average.
• phaseEst is the raw phase returned by the phase estimation algorithm. This useful for diagnosing aliasing when it occurs due to a trotterStep value that is too large.