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"""Tests  _bravyi_kitaev_fast_test.py."""
13
14import os
15import unittest
16
17import numpy
18
19from openfermion.config import DATA_DIRECTORY
20from openfermion.chem import MolecularData
21from openfermion.ops.operators import FermionOperator, QubitOperator
22from openfermion.ops.representations import InteractionOperator
23from openfermion.transforms.opconversions import (bksf, get_fermion_operator,
24                                                  normal_ordered)
25from openfermion.transforms.opconversions.jordan_wigner import (
26    jordan_wigner, jordan_wigner_one_body)
27from openfermion.linalg import get_sparse_operator, eigenspectrum
28from openfermion.utils import count_qubits
29
30
31class bravyi_kitaev_fastTransformTest(unittest.TestCase):
32
33    def setUp(self):
34        geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))]
35        basis = 'sto-3g'
36        multiplicity = 1
37        filename = os.path.join(DATA_DIRECTORY, 'H2_sto-3g_singlet_0.7414')
38        self.molecule = MolecularData(geometry,
39                                      basis,
40                                      multiplicity,
41                                      filename=filename)
42        self.molecule.load()
43
44        # Get molecular Hamiltonian.
45        self.molecular_hamiltonian = self.molecule.get_molecular_hamiltonian()
46
47        # Get FCI RDM.
48        self.fci_rdm = self.molecule.get_molecular_rdm(use_fci=1)
49        # Get explicit coefficients.
50        self.nuclear_repulsion = self.molecular_hamiltonian.constant
51        self.one_body = self.molecular_hamiltonian.one_body_tensor
52        self.two_body = self.molecular_hamiltonian.two_body_tensor
53
54        # Get fermion Hamiltonian.
55        self.fermion_hamiltonian = normal_ordered(
56            get_fermion_operator(self.molecular_hamiltonian))
57
58        # Get qubit Hamiltonian.
59        self.qubit_hamiltonian = jordan_wigner(self.fermion_hamiltonian)
60
61        # Get the sparse matrix.
62        self.hamiltonian_matrix = get_sparse_operator(
63            self.molecular_hamiltonian)
64
65    def test_bad_input(self):
66        with self.assertRaises(TypeError):
67            bksf.bravyi_kitaev_fast(FermionOperator())
68
69    def test_bravyi_kitaev_fast_edgeoperator_Bi(self):
70        # checking the edge operators
71        edge_matrix = numpy.triu(numpy.ones((4, 4)))
72        edge_matrix_indices = numpy.array(
73            numpy.nonzero(
74                numpy.triu(edge_matrix) - numpy.diag(numpy.diag(edge_matrix))))
75
76        correct_operators_b0 = ((0, 'Z'), (1, 'Z'), (2, 'Z'))
77        correct_operators_b1 = ((0, 'Z'), (3, 'Z'), (4, 'Z'))
78        correct_operators_b2 = ((1, 'Z'), (3, 'Z'), (5, 'Z'))
79        correct_operators_b3 = ((2, 'Z'), (4, 'Z'), (5, 'Z'))
80
81        qterm_b0 = QubitOperator(correct_operators_b0, 1)
82        qterm_b1 = QubitOperator(correct_operators_b1, 1)
83        qterm_b2 = QubitOperator(correct_operators_b2, 1)
84        qterm_b3 = QubitOperator(correct_operators_b3, 1)
85        self.assertTrue(
86            qterm_b0 == bksf.edge_operator_b(edge_matrix_indices, 0))
87        self.assertTrue(
88            qterm_b1 == bksf.edge_operator_b(edge_matrix_indices, 1))
89        self.assertTrue(
90            qterm_b2 == bksf.edge_operator_b(edge_matrix_indices, 2))
91        self.assertTrue(
92            qterm_b3 == bksf.edge_operator_b(edge_matrix_indices, 3))
93
94    def test_bravyi_kitaev_fast_edgeoperator_Aij(self):
95        # checking the edge operators
96        edge_matrix = numpy.triu(numpy.ones((4, 4)))
97        edge_matrix_indices = numpy.array(
98            numpy.nonzero(
99                numpy.triu(edge_matrix) - numpy.diag(numpy.diag(edge_matrix))))
100        correct_operators_a01 = ((0, 'X'),)
101        correct_operators_a02 = ((0, 'Z'), (1, 'X'))
102        correct_operators_a03 = ((0, 'Z'), (1, 'Z'), (2, 'X'))
103        correct_operators_a12 = ((0, 'Z'), (1, 'Z'), (3, 'X'))
104        correct_operators_a13 = ((0, 'Z'), (2, 'Z'), (3, 'Z'), (4, 'X'))
105        correct_operators_a23 = ((1, 'Z'), (2, 'Z'), (3, 'Z'), (4, 'Z'), (5,
106                                                                          'X'))
107
108        qterm_a01 = QubitOperator(correct_operators_a01, 1)
109        qterm_a02 = QubitOperator(correct_operators_a02, 1)
110        qterm_a03 = QubitOperator(correct_operators_a03, 1)
111        qterm_a12 = QubitOperator(correct_operators_a12, 1)
112        qterm_a13 = QubitOperator(correct_operators_a13, 1)
113        qterm_a23 = QubitOperator(correct_operators_a23, 1)
114
115        self.assertTrue(
116            qterm_a01 == bksf.edge_operator_aij(edge_matrix_indices, 0, 1))
117        self.assertTrue(
118            qterm_a02 == bksf.edge_operator_aij(edge_matrix_indices, 0, 2))
119        self.assertTrue(
120            qterm_a03 == bksf.edge_operator_aij(edge_matrix_indices, 0, 3))
121        self.assertTrue(
122            qterm_a12 == bksf.edge_operator_aij(edge_matrix_indices, 1, 2))
123        self.assertTrue(
124            qterm_a13 == bksf.edge_operator_aij(edge_matrix_indices, 1, 3))
125        self.assertTrue(
126            qterm_a23 == bksf.edge_operator_aij(edge_matrix_indices, 2, 3))
127
128    def test_bravyi_kitaev_fast_jw_number_operator(self):
129        # bksf algorithm allows for even number of particles. So, compare the
130        # spectrum of number operator from jordan-wigner and bksf algorithm
131        # to make sure half of the jordan-wigner number operator spectrum
132        # can be found in bksf number operator spectrum.
133        bravyi_kitaev_fast_n = bksf.number_operator(self.molecular_hamiltonian)
134        jw_n = QubitOperator()
135        n_qubits = count_qubits(self.molecular_hamiltonian)
136        for i in range(n_qubits):
137            jw_n += jordan_wigner_one_body(i, i)
138        jw_eig_spec = eigenspectrum(jw_n)
139        bravyi_kitaev_fast_eig_spec = eigenspectrum(bravyi_kitaev_fast_n)
140        evensector = 0
141        for i in range(numpy.size(jw_eig_spec)):
142            if bool(
143                    numpy.size(
144                        numpy.where(
145                            jw_eig_spec[i] == bravyi_kitaev_fast_eig_spec))):
146                evensector += 1
147        self.assertEqual(evensector, 2**(n_qubits - 1))
148
149    def test_bravyi_kitaev_fast_jw_hamiltonian(self):
150        # make sure half of the jordan-wigner Hamiltonian eigenspectrum can
151        # be found in bksf Hamiltonian eigenspectrum.
152        n_qubits = count_qubits(self.molecular_hamiltonian)
153        bravyi_kitaev_fast_H = bksf.bravyi_kitaev_fast(
154            self.molecular_hamiltonian)
155        jw_H = jordan_wigner(self.molecular_hamiltonian)
156        bravyi_kitaev_fast_H_eig = eigenspectrum(bravyi_kitaev_fast_H)
157        jw_H_eig = eigenspectrum(jw_H)
158        bravyi_kitaev_fast_H_eig = bravyi_kitaev_fast_H_eig.round(5)
159        jw_H_eig = jw_H_eig.round(5)
160        evensector = 0
161        for i in range(numpy.size(jw_H_eig)):
162            if bool(
163                    numpy.size(
164                        numpy.where(jw_H_eig[i] == bravyi_kitaev_fast_H_eig))):
165                evensector += 1
166        self.assertEqual(evensector, 2**(n_qubits - 1))
167
168    def test_bravyi_kitaev_fast_generate_fermions(self):
169        n_qubits = count_qubits(self.molecular_hamiltonian)
170        # test for generating two fermions
171        edge_matrix = bksf.bravyi_kitaev_fast_edge_matrix(
172            self.molecular_hamiltonian)
173        edge_matrix_indices = numpy.array(
174            numpy.nonzero(
175                numpy.triu(edge_matrix) - numpy.diag(numpy.diag(edge_matrix))))
176        fermion_generation_operator = bksf.generate_fermions(
177            edge_matrix_indices, 2, 3)
178        fermion_generation_sp_matrix = get_sparse_operator(
179            fermion_generation_operator)
180        fermion_generation_matrix = fermion_generation_sp_matrix.toarray()
181        bksf_vacuum_state_operator = bksf.vacuum_operator(edge_matrix_indices)
182        bksf_vacuum_state_sp_matrix = get_sparse_operator(
183            bksf_vacuum_state_operator)
184        bksf_vacuum_state_matrix = bksf_vacuum_state_sp_matrix.toarray()
185        vacuum_state = numpy.zeros((2**(n_qubits), 1))
186        vacuum_state[0] = 1.
187        bksf_vacuum_state = numpy.dot(bksf_vacuum_state_matrix, vacuum_state)
188        two_fermion_state = numpy.dot(fermion_generation_matrix,
189                                      bksf_vacuum_state)
190        # using total number operator to check the number of fermions generated
191        tot_number_operator = bksf.number_operator(self.molecular_hamiltonian)
192        number_operator_sp_matrix = get_sparse_operator(tot_number_operator)
193        number_operator_matrix = number_operator_sp_matrix.toarray()
194        tot_fermions = numpy.dot(
195            two_fermion_state.conjugate().T,
196            numpy.dot(number_operator_matrix, two_fermion_state))
197        # checking the occupation number of site 2 and 3
198        number_operator_2 = bksf.number_operator(self.molecular_hamiltonian, 2)
199        number_operator_3 = bksf.number_operator(self.molecular_hamiltonian, 3)
200        number_operator_23 = number_operator_2 + number_operator_3
201        number_operator_23_sp_matrix = get_sparse_operator(number_operator_23)
202        number_operator_23_matrix = number_operator_23_sp_matrix.toarray()
203        tot_23_fermions = numpy.dot(
204            two_fermion_state.conjugate().T,
205            numpy.dot(number_operator_23_matrix, two_fermion_state))
206        self.assertTrue(2.0 - float(tot_fermions.real) < 1e-13)
207        self.assertTrue(2.0 - float(tot_23_fermions.real) < 1e-13)
208
209    def test_bravyi_kitaev_fast_excitation_terms(self):
210        # Testing on-site and excitation terms in Hamiltonian
211        constant = 0
212        one_body = numpy.array([[1, 2, 0, 3], [2, 1, 2, 0], [0, 2, 1, 2.5],
213                                [3, 0, 2.5, 1]])
214        # No Coloumb interaction
215        two_body = numpy.zeros((4, 4, 4, 4))
216        molecular_hamiltonian = InteractionOperator(constant, one_body,
217                                                    two_body)
218        n_qubits = count_qubits(molecular_hamiltonian)
219        # comparing the eigenspectrum of Hamiltonian
220        bravyi_kitaev_fast_H = bksf.bravyi_kitaev_fast(molecular_hamiltonian)
221        jw_H = jordan_wigner(molecular_hamiltonian)
222        bravyi_kitaev_fast_H_eig = eigenspectrum(bravyi_kitaev_fast_H)
223        jw_H_eig = eigenspectrum(jw_H)
224        bravyi_kitaev_fast_H_eig = bravyi_kitaev_fast_H_eig.round(5)
225        jw_H_eig = jw_H_eig.round(5)
226        evensector_H = 0
227        for i in range(numpy.size(jw_H_eig)):
228            if bool(
229                    numpy.size(
230                        numpy.where(jw_H_eig[i] == bravyi_kitaev_fast_H_eig))):
231                evensector_H += 1
232
233        # comparing eigenspectrum of number operator
234        bravyi_kitaev_fast_n = bksf.number_operator(molecular_hamiltonian)
235        jw_n = QubitOperator()
236        n_qubits = count_qubits(molecular_hamiltonian)
237        for i in range(n_qubits):
238            jw_n += jordan_wigner_one_body(i, i)
239        jw_eig_spec = eigenspectrum(jw_n)
240        bravyi_kitaev_fast_eig_spec = eigenspectrum(bravyi_kitaev_fast_n)
241        evensector_n = 0
242        for i in range(numpy.size(jw_eig_spec)):
243            if bool(
244                    numpy.size(
245                        numpy.where(
246                            jw_eig_spec[i] == bravyi_kitaev_fast_eig_spec))):
247                evensector_n += 1
248        self.assertEqual(evensector_H, 2**(n_qubits - 1))
249        self.assertEqual(evensector_n, 2**(n_qubits - 1))
250
251    def test_bravyi_kitaev_fast_number_excitation_operator(self):
252        # using hydrogen Hamiltonian and introducing some number operator terms
253        constant = 0
254        one_body = numpy.zeros((4, 4))
255        one_body[(0, 0)] = .4
256        one_body[(1, 1)] = .5
257        one_body[(2, 2)] = .6
258        one_body[(3, 3)] = .7
259        two_body = self.molecular_hamiltonian.two_body_tensor
260        # initiating number operator terms for all the possible cases
261        two_body[(1, 2, 3, 1)] = 0.1
262        two_body[(1, 3, 2, 1)] = 0.1
263        two_body[(1, 2, 1, 3)] = 0.15
264        two_body[(3, 1, 2, 1)] = 0.15
265        two_body[(0, 2, 2, 1)] = 0.09
266        two_body[(1, 2, 2, 0)] = 0.09
267        two_body[(1, 2, 3, 2)] = 0.11
268        two_body[(2, 3, 2, 1)] = 0.11
269        two_body[(2, 2, 2, 2)] = 0.1
270        molecular_hamiltonian = InteractionOperator(constant, one_body,
271                                                    two_body)
272        # comparing the eigenspectrum of Hamiltonian
273        n_qubits = count_qubits(molecular_hamiltonian)
274        bravyi_kitaev_fast_H = bksf.bravyi_kitaev_fast(molecular_hamiltonian)
275        jw_H = jordan_wigner(molecular_hamiltonian)
276        bravyi_kitaev_fast_H_eig = eigenspectrum(bravyi_kitaev_fast_H)
277        jw_H_eig = eigenspectrum(jw_H)
278        bravyi_kitaev_fast_H_eig = bravyi_kitaev_fast_H_eig.round(5)
279        jw_H_eig = jw_H_eig.round(5)
280        evensector_H = 0
281        for i in range(numpy.size(jw_H_eig)):
282            if bool(
283                    numpy.size(
284                        numpy.where(jw_H_eig[i] == bravyi_kitaev_fast_H_eig))):
285                evensector_H += 1
286
287        # comparing eigenspectrum of number operator
288        bravyi_kitaev_fast_n = bksf.number_operator(molecular_hamiltonian)
289        jw_n = QubitOperator()
290        n_qubits = count_qubits(molecular_hamiltonian)
291        for i in range(n_qubits):
292            jw_n += jordan_wigner_one_body(i, i)
293        jw_eig_spec = eigenspectrum(jw_n)
294        bravyi_kitaev_fast_eig_spec = eigenspectrum(bravyi_kitaev_fast_n)
295        evensector_n = 0
296        for i in range(numpy.size(jw_eig_spec)):
297            if bool(
298                    numpy.size(
299                        numpy.where(
300                            jw_eig_spec[i] == bravyi_kitaev_fast_eig_spec))):
301                evensector_n += 1
302        self.assertEqual(evensector_H, 2**(n_qubits - 1))
303        self.assertEqual(evensector_n, 2**(n_qubits - 1))
304