1#   Licensed under the Apache License, Version 2.0 (the "License");
2#   you may not use this file except in compliance with the License.
3#   You may obtain a copy of the License at
4#
5#       http://www.apache.org/licenses/LICENSE-2.0
6#
7#   Unless required by applicable law or agreed to in writing, software
8#   distributed under the License is distributed on an "AS IS" BASIS,
9#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
10#   See the License for the specific language governing permissions and
11#   limitations under the License.
12""" Unit tests of remove_symmetry_qubits(). Test that it preserves the
13    correct ground state energy, and reduces the number of qubits required
14    by two.
15"""
16
17import unittest
18
19from openfermion.hamiltonians import fermi_hubbard
20from openfermion.chem import MolecularData
21from openfermion.transforms.opconversions import get_fermion_operator
22from openfermion.linalg.sparse_tools import (
23    get_sparse_operator, jw_get_ground_state_at_particle_number)
24from openfermion.linalg import eigenspectrum
25from openfermion.ops.operators import FermionOperator
26
27from openfermion.transforms.opconversions.remove_symmetry_qubits import (
28    symmetry_conserving_bravyi_kitaev)
29
30
31def LiH_sto3g():
32    """ Generates the Hamiltonian for LiH in
33        the STO-3G basis, at a distance of
34        1.45 A.
35    """
36    geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))]
37    molecule = MolecularData(geometry, 'sto-3g', 1, description="1.45")
38    molecule.load()
39    molecular_hamiltonian = molecule.get_molecular_hamiltonian()
40    hamiltonian = get_fermion_operator(molecular_hamiltonian)
41    num_electrons = molecule.n_electrons
42    num_orbitals = 2 * molecule.n_orbitals
43
44    return hamiltonian, num_orbitals, num_electrons
45
46
47def set_1D_hubbard(x_dim):
48    """ Returns a 1D Fermi-Hubbard Hamiltonian.
49    """
50    y_dim = 1
51    tunneling = 1.0
52    coulomb = 1.0
53    magnetic_field = 0.0
54    chemical_potential = 0.0
55    periodic = False
56    spinless = False
57
58    num_orbitals = 2 * x_dim * y_dim
59
60    hubbard_model = fermi_hubbard(x_dim, y_dim, tunneling, coulomb,
61                                  chemical_potential, magnetic_field, periodic,
62                                  spinless)
63
64    return hubbard_model, num_orbitals
65
66
67def number_of_qubits(qubit_hamiltonian, unreduced_orbitals):
68    """ Returns the number of qubits that a
69        qubit Hamiltonian acts upon.
70    """
71    max_orbital = 0
72    for i in range(unreduced_orbitals):
73        for term in qubit_hamiltonian.terms:
74            if (i, 'X') in term or (i, 'Y') in term or (i, 'Z') in term:
75                max_orbital = max_orbital + 1
76                break
77
78    return max_orbital
79
80
81class ReduceSymmetryQubitsTest(unittest.TestCase):
82
83    # Check whether fermionic and reduced qubit Hamiltonians
84    # have the same energy for LiH
85    def test_energy_reduce_symmetry_qubits(self):
86        # Generate the fermionic Hamiltonians,
87        # number of orbitals and number of electrons.
88        lih_sto_hamil, lih_sto_numorb, lih_sto_numel = LiH_sto3g()
89
90        # Use test function to reduce the qubits.
91        lih_sto_qbt = (symmetry_conserving_bravyi_kitaev(
92            lih_sto_hamil, lih_sto_numorb, lih_sto_numel))
93
94        self.assertAlmostEqual(
95            eigenspectrum(lih_sto_qbt)[0],
96            eigenspectrum(lih_sto_hamil)[0])
97
98    # Check that the qubit Hamiltonian acts on two fewer qubits
99    # for LiH.
100    def test_orbnum_reduce_symmetry_qubits(self):
101        # Generate the fermionic Hamiltonians,
102        # number of orbitals and number of electrons.
103        lih_sto_hamil, lih_sto_numorb, lih_sto_numel = LiH_sto3g()
104
105        # Use test function to reduce the qubits.
106        lih_sto_qbt = (symmetry_conserving_bravyi_kitaev(
107            lih_sto_hamil, lih_sto_numorb, lih_sto_numel))
108
109        self.assertEqual(number_of_qubits(lih_sto_qbt, lih_sto_numorb),
110                         lih_sto_numorb - 2)
111
112    # Check ValueErrors arise correctly.
113    def test_errors_reduce_symmetry_qubits(self):
114        # Generate the fermionic Hamiltonians,
115        # number of orbitals and number of electrons.
116        lih_sto_hamil, lih_sto_numorb, lih_sto_numel = LiH_sto3g()
117
118        with self.assertRaises(ValueError):
119            symmetry_conserving_bravyi_kitaev(0, lih_sto_numorb, lih_sto_numel)
120        with self.assertRaises(ValueError):
121            symmetry_conserving_bravyi_kitaev(lih_sto_hamil, 1.5, lih_sto_numel)
122        with self.assertRaises(ValueError):
123            symmetry_conserving_bravyi_kitaev(lih_sto_hamil, lih_sto_numorb,
124                                              3.6)
125
126    # Check energy is the same for Fermi-Hubbard model.
127    def test_hubbard_reduce_symmetry_qubits(self):
128        for i in range(4):
129            n_sites = i + 2
130            n_ferm = n_sites
131            hub_hamil, n_orb = set_1D_hubbard(n_sites)
132
133            # Use test function to reduce the qubits.
134            hub_qbt = (symmetry_conserving_bravyi_kitaev(
135                hub_hamil, n_orb, n_ferm))
136
137            sparse_op = get_sparse_operator(hub_hamil)
138            ground_energy, _ = jw_get_ground_state_at_particle_number(
139                sparse_op, n_ferm)
140
141            self.assertAlmostEqual(eigenspectrum(hub_qbt)[0], ground_energy)
142
143    # Check if single operator of a Hamiltonian is transformed consistently,
144    # e.g. system with 2 particles, 4 spin-orbitals, and want to determine
145    # elements of one- and two-particle reduced density matrix
146    def test_single_operator(self):
147        # Dummy operator acting only on 2 qubits of overall 4-qubit system
148        op = FermionOperator("0^ 1^ 1 0") + FermionOperator("1^ 0^ 0 1")
149        trafo_op = symmetry_conserving_bravyi_kitaev(op,
150                                                     active_fermions=2,
151                                                     active_orbitals=4)
152        # Check via eigenspectrum -- needs to stay the same
153        e_op = eigenspectrum(op)
154        e_trafo = eigenspectrum(trafo_op)
155        # Check eigenvalues
156        self.assertSequenceEqual(e_op.tolist(), e_trafo.tolist())
157