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 for sparse_tools.py.""" 13import os 14 15import unittest 16 17import itertools 18 19import numpy 20from numpy.linalg import multi_dot 21 22from scipy.linalg import eigh, norm 23from scipy.sparse import csc_matrix 24import scipy 25from scipy.special import comb 26 27from openfermion.config import DATA_DIRECTORY 28from openfermion.chem import MolecularData 29from openfermion.hamiltonians import (fermi_hubbard, jellium_model, 30 wigner_seitz_length_scale, 31 number_operator) 32from openfermion.ops.operators import (FermionOperator, BosonOperator, 33 QuadOperator, QubitOperator, 34 MajoranaOperator) 35from openfermion.ops.representations import DiagonalCoulombHamiltonian 36from openfermion.transforms.opconversions import (jordan_wigner, 37 get_fermion_operator, 38 normal_ordered) 39from openfermion.transforms.repconversions import (fourier_transform, 40 get_interaction_operator) 41 42from openfermion.testing.testing_utils import random_hermitian_matrix 43from openfermion.linalg.sparse_tools import ( 44 get_sparse_operator, get_ground_state, eigenspectrum, expectation, 45 jw_number_restrict_state, inner_product, jw_sz_restrict_state, 46 jw_get_ground_state_at_particle_number, jw_sparse_givens_rotation, 47 pauli_matrix_map, sparse_eigenspectrum, jordan_wigner_sparse, 48 kronecker_operators, identity_csc, pauli_x_csc, qubit_operator_sparse, 49 get_linear_qubit_operator_diagonal, jw_number_indices, 50 get_number_preserving_sparse_operator, jw_configuration_state, 51 jw_hartree_fock_state, jw_sz_restrict_operator, jw_sz_indices, 52 jw_number_restrict_operator, variance, 53 expectation_computational_basis_state, 54 expectation_db_operator_with_pw_basis_state, get_gap, boson_ladder_sparse, 55 boson_operator_sparse, single_quad_op_sparse, _iterate_basis_) 56 57from openfermion.utils.operator_utils import (is_hermitian, count_qubits, 58 hermitian_conjugated) 59from openfermion.utils.indexing import up_index, down_index 60from openfermion.utils.grid import Grid 61from openfermion.hamiltonians.jellium_hf_state import ( 62 lowest_single_particle_energy_states) 63from openfermion.linalg.linear_qubit_operator import LinearQubitOperator 64 65 66class EigenSpectrumTest(unittest.TestCase): 67 68 def setUp(self): 69 self.n_qubits = 5 70 self.majorana_operator = MajoranaOperator((1, 4, 9)) 71 self.fermion_term = FermionOperator('1^ 2^ 3 4', -3.17) 72 self.fermion_operator = self.fermion_term + hermitian_conjugated( 73 self.fermion_term) 74 self.qubit_operator = jordan_wigner(self.fermion_operator) 75 self.interaction_operator = get_interaction_operator( 76 self.fermion_operator) 77 78 def test_eigenspectrum(self): 79 fermion_eigenspectrum = eigenspectrum(self.fermion_operator) 80 qubit_eigenspectrum = eigenspectrum(self.qubit_operator) 81 interaction_eigenspectrum = eigenspectrum(self.interaction_operator) 82 for i in range(2**self.n_qubits): 83 self.assertAlmostEqual(fermion_eigenspectrum[i], 84 qubit_eigenspectrum[i]) 85 self.assertAlmostEqual(fermion_eigenspectrum[i], 86 interaction_eigenspectrum[i]) 87 88 with self.assertRaises(TypeError): 89 _ = eigenspectrum(BosonOperator()) 90 91 with self.assertRaises(TypeError): 92 _ = eigenspectrum(QuadOperator()) 93 94 95class SparseOperatorTest(unittest.TestCase): 96 97 def test_kronecker_operators(self): 98 99 self.assertAlmostEqual( 100 0, 101 numpy.amax( 102 numpy.absolute( 103 kronecker_operators(3 * [identity_csc]) - 104 kronecker_operators(3 * [pauli_x_csc])**2))) 105 106 def test_qubit_jw_fermion_integration(self): 107 108 # Initialize a random fermionic operator. 109 fermion_operator = FermionOperator(((3, 1), (2, 1), (1, 0), (0, 0)), 110 -4.3) 111 fermion_operator += FermionOperator(((3, 1), (1, 0)), 8.17) 112 fermion_operator += 3.2 * FermionOperator() 113 114 # Map to qubits and compare matrix versions. 115 qubit_operator = jordan_wigner(fermion_operator) 116 qubit_sparse = get_sparse_operator(qubit_operator) 117 qubit_spectrum = sparse_eigenspectrum(qubit_sparse) 118 fermion_sparse = jordan_wigner_sparse(fermion_operator) 119 fermion_spectrum = sparse_eigenspectrum(fermion_sparse) 120 self.assertAlmostEqual( 121 0., numpy.amax(numpy.absolute(fermion_spectrum - qubit_spectrum))) 122 123 124class JordanWignerSparseTest(unittest.TestCase): 125 126 def test_jw_sparse_0create(self): 127 expected = csc_matrix(([1], ([1], [0])), shape=(2, 2)) 128 self.assertTrue( 129 numpy.allclose( 130 jordan_wigner_sparse(FermionOperator('0^')).A, expected.A)) 131 132 def test_jw_sparse_1annihilate(self): 133 expected = csc_matrix(([1, -1], ([0, 2], [1, 3])), shape=(4, 4)) 134 self.assertTrue( 135 numpy.allclose( 136 jordan_wigner_sparse(FermionOperator('1')).A, expected.A)) 137 138 def test_jw_sparse_0create_2annihilate(self): 139 expected = csc_matrix(([-1j, 1j], ([4, 6], [1, 3])), 140 shape=(8, 8), 141 dtype=numpy.complex128) 142 self.assertTrue( 143 numpy.allclose( 144 jordan_wigner_sparse(FermionOperator('0^ 2', -1j)).A, 145 expected.A)) 146 147 def test_jw_sparse_0create_3annihilate(self): 148 expected = csc_matrix( 149 ([-1j, 1j, 1j, -1j], ([8, 10, 12, 14], [1, 3, 5, 7])), 150 shape=(16, 16), 151 dtype=numpy.complex128) 152 self.assertTrue( 153 numpy.allclose( 154 jordan_wigner_sparse(FermionOperator('0^ 3', -1j)).A, 155 expected.A)) 156 157 def test_jw_sparse_twobody(self): 158 expected = csc_matrix(([1, 1], ([6, 14], [5, 13])), shape=(16, 16)) 159 self.assertTrue( 160 numpy.allclose( 161 jordan_wigner_sparse(FermionOperator('2^ 1^ 1 3')).A, 162 expected.A)) 163 164 def test_qubit_operator_sparse_n_qubits_too_small(self): 165 with self.assertRaises(ValueError): 166 qubit_operator_sparse(QubitOperator('X3'), 1) 167 168 def test_qubit_operator_sparse_n_qubits_not_specified(self): 169 expected = csc_matrix(([1, 1, 1, 1], ([1, 0, 3, 2], [0, 1, 2, 3])), 170 shape=(4, 4)) 171 self.assertTrue( 172 numpy.allclose( 173 qubit_operator_sparse(QubitOperator('X1')).A, expected.A)) 174 175 def test_get_linear_qubit_operator_diagonal_wrong_n(self): 176 """Testing with wrong n_qubits.""" 177 with self.assertRaises(ValueError): 178 get_linear_qubit_operator_diagonal(QubitOperator('X3'), 1) 179 180 def test_get_linear_qubit_operator_diagonal_0(self): 181 """Testing with zero term.""" 182 qubit_operator = QubitOperator.zero() 183 vec_expected = numpy.zeros(8) 184 185 self.assertTrue( 186 numpy.allclose( 187 get_linear_qubit_operator_diagonal(qubit_operator, 3), 188 vec_expected)) 189 190 def test_get_linear_qubit_operator_diagonal_zero(self): 191 """Get zero diagonals from get_linear_qubit_operator_diagonal.""" 192 qubit_operator = QubitOperator('X0 Y1') 193 vec_expected = numpy.zeros(4) 194 195 self.assertTrue( 196 numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator), 197 vec_expected)) 198 199 def test_get_linear_qubit_operator_diagonal_non_zero(self): 200 """Get non zero diagonals from get_linear_qubit_operator_diagonal.""" 201 qubit_operator = QubitOperator('Z0 Z2') 202 vec_expected = numpy.array([1, -1, 1, -1, -1, 1, -1, 1]) 203 204 self.assertTrue( 205 numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator), 206 vec_expected)) 207 208 def test_get_linear_qubit_operator_diagonal_cmp_zero(self): 209 """Compare get_linear_qubit_operator_diagonal with 210 get_linear_qubit_operator.""" 211 qubit_operator = QubitOperator('Z1 X2 Y5') 212 vec_expected = numpy.diag( 213 LinearQubitOperator(qubit_operator) * numpy.eye(2**6)) 214 215 self.assertTrue( 216 numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator), 217 vec_expected)) 218 219 def test_get_linear_qubit_operator_diagonal_cmp_non_zero(self): 220 """Compare get_linear_qubit_operator_diagonal with 221 get_linear_qubit_operator.""" 222 qubit_operator = QubitOperator('Z1 Z2 Z5') 223 vec_expected = numpy.diag( 224 LinearQubitOperator(qubit_operator) * numpy.eye(2**6)) 225 226 self.assertTrue( 227 numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator), 228 vec_expected)) 229 230 231class ComputationalBasisStateTest(unittest.TestCase): 232 233 def test_computational_basis_state(self): 234 comp_basis_state = jw_configuration_state([0, 2, 5], 7) 235 self.assertAlmostEqual(comp_basis_state[82], 1.) 236 self.assertAlmostEqual(sum(comp_basis_state), 1.) 237 238 239class JWHartreeFockStateTest(unittest.TestCase): 240 241 def test_jw_hartree_fock_state(self): 242 hartree_fock_state = jw_hartree_fock_state(3, 7) 243 self.assertAlmostEqual(hartree_fock_state[112], 1.) 244 self.assertAlmostEqual(sum(hartree_fock_state), 1.) 245 246 247class JWNumberIndicesTest(unittest.TestCase): 248 249 def test_jw_sparse_index(self): 250 """Test the indexing scheme for selecting specific particle numbers""" 251 expected = [1, 2] 252 calculated_indices = jw_number_indices(1, 2) 253 self.assertEqual(expected, calculated_indices) 254 255 expected = [3] 256 calculated_indices = jw_number_indices(2, 2) 257 self.assertEqual(expected, calculated_indices) 258 259 def test_jw_number_indices(self): 260 n_qubits = numpy.random.randint(1, 12) 261 n_particles = numpy.random.randint(n_qubits + 1) 262 263 number_indices = jw_number_indices(n_particles, n_qubits) 264 subspace_dimension = len(number_indices) 265 266 self.assertEqual(subspace_dimension, comb(n_qubits, n_particles)) 267 268 for index in number_indices: 269 binary_string = bin(index)[2:].zfill(n_qubits) 270 n_ones = binary_string.count('1') 271 self.assertEqual(n_ones, n_particles) 272 273 274class JWSzIndicesTest(unittest.TestCase): 275 276 def test_jw_sz_indices(self): 277 """Test the indexing scheme for selecting specific sz value""" 278 279 def sz_integer(bitstring): 280 """Computes the total number of occupied up sites 281 minus the total number of occupied down sites.""" 282 n_sites = len(bitstring) // 2 283 284 n_up = len([ 285 site for site in range(n_sites) 286 if bitstring[up_index(site)] == '1' 287 ]) 288 n_down = len([ 289 site for site in range(n_sites) 290 if bitstring[down_index(site)] == '1' 291 ]) 292 293 return n_up - n_down 294 295 def jw_sz_indices_brute_force(sz_value, n_qubits): 296 """Computes the correct indices by brute force.""" 297 indices = [] 298 for bitstring in itertools.product(['0', '1'], repeat=n_qubits): 299 if (sz_integer(bitstring) == int(2 * sz_value)): 300 indices.append(int(''.join(bitstring), 2)) 301 302 return indices 303 304 # General test 305 n_sites = numpy.random.randint(1, 10) 306 n_qubits = 2 * n_sites 307 sz_int = ((-1)**numpy.random.randint(2) * 308 numpy.random.randint(n_sites + 1)) 309 sz_value = sz_int / 2. 310 311 correct_indices = jw_sz_indices_brute_force(sz_value, n_qubits) 312 subspace_dimension = len(correct_indices) 313 314 calculated_indices = jw_sz_indices(sz_value, n_qubits) 315 316 self.assertEqual(len(calculated_indices), subspace_dimension) 317 318 for index in calculated_indices: 319 binary_string = bin(index)[2:].zfill(n_qubits) 320 self.assertEqual(sz_integer(binary_string), sz_int) 321 322 # Test fixing particle number 323 n_particles = abs(sz_int) 324 325 correct_indices = [ 326 index for index in correct_indices 327 if bin(index)[2:].count('1') == n_particles 328 ] 329 subspace_dimension = len(correct_indices) 330 331 calculated_indices = jw_sz_indices(sz_value, 332 n_qubits, 333 n_electrons=n_particles) 334 335 self.assertEqual(len(calculated_indices), subspace_dimension) 336 337 for index in calculated_indices: 338 binary_string = bin(index)[2:].zfill(n_qubits) 339 self.assertEqual(sz_integer(binary_string), sz_int) 340 self.assertEqual(binary_string.count('1'), n_particles) 341 342 # Test exceptions 343 with self.assertRaises(ValueError): 344 jw_sz_indices(3, 3) 345 346 with self.assertRaises(ValueError): 347 jw_sz_indices(3.1, 4) 348 349 with self.assertRaises(ValueError): 350 jw_sz_indices(1.5, 8, n_electrons=6) 351 352 with self.assertRaises(ValueError): 353 jw_sz_indices(1.5, 8, n_electrons=1) 354 355 356class JWNumberRestrictOperatorTest(unittest.TestCase): 357 358 def test_jw_restrict_operator(self): 359 """Test the scheme for restricting JW encoded operators to number""" 360 # Make a Hamiltonian that cares mostly about number of electrons 361 n_qubits = 4 362 target_electrons = 2 363 penalty_const = 10. 364 number_sparse = jordan_wigner_sparse(number_operator(n_qubits)) 365 bias_sparse = jordan_wigner_sparse( 366 sum([ 367 FermionOperator(((i, 1), (i, 0)), 1.0) for i in range(n_qubits) 368 ], FermionOperator())) 369 hamiltonian_sparse = penalty_const * ( 370 number_sparse - 371 target_electrons * scipy.sparse.identity(2**n_qubits) 372 ).dot(number_sparse - target_electrons * 373 scipy.sparse.identity(2**n_qubits)) + bias_sparse 374 375 restricted_hamiltonian = jw_number_restrict_operator( 376 hamiltonian_sparse, target_electrons, n_qubits) 377 true_eigvals, _ = eigh(hamiltonian_sparse.A) 378 test_eigvals, _ = eigh(restricted_hamiltonian.A) 379 380 self.assertAlmostEqual(norm(true_eigvals[:6] - test_eigvals[:6]), 0.0) 381 382 def test_jw_restrict_operator_hopping_to_1_particle(self): 383 hop = FermionOperator('3^ 1') + FermionOperator('1^ 3') 384 hop_sparse = jordan_wigner_sparse(hop, n_qubits=4) 385 hop_restrict = jw_number_restrict_operator(hop_sparse, 1, n_qubits=4) 386 expected = csc_matrix(([1, 1], ([0, 2], [2, 0])), shape=(4, 4)) 387 388 self.assertTrue(numpy.allclose(hop_restrict.A, expected.A)) 389 390 def test_jw_restrict_operator_interaction_to_1_particle(self): 391 interaction = FermionOperator('3^ 2^ 4 1') 392 interaction_sparse = jordan_wigner_sparse(interaction, n_qubits=6) 393 interaction_restrict = jw_number_restrict_operator(interaction_sparse, 394 1, 395 n_qubits=6) 396 expected = csc_matrix(([], ([], [])), shape=(6, 6)) 397 398 self.assertTrue(numpy.allclose(interaction_restrict.A, expected.A)) 399 400 def test_jw_restrict_operator_interaction_to_2_particles(self): 401 interaction = (FermionOperator('3^ 2^ 4 1') + 402 FermionOperator('4^ 1^ 3 2')) 403 interaction_sparse = jordan_wigner_sparse(interaction, n_qubits=6) 404 interaction_restrict = jw_number_restrict_operator(interaction_sparse, 405 2, 406 n_qubits=6) 407 408 dim = 6 * 5 // 2 # shape of new sparse array 409 410 # 3^ 2^ 4 1 maps 2**4 + 2 = 18 to 2**3 + 2**2 = 12 and vice versa; 411 # in the 2-particle subspace (1, 4) and (2, 3) are 7th and 9th. 412 expected = csc_matrix(([-1, -1], ([7, 9], [9, 7])), shape=(dim, dim)) 413 414 self.assertTrue(numpy.allclose(interaction_restrict.A, expected.A)) 415 416 def test_jw_restrict_operator_hopping_to_1_particle_default_nqubits(self): 417 interaction = (FermionOperator('3^ 2^ 4 1') + 418 FermionOperator('4^ 1^ 3 2')) 419 interaction_sparse = jordan_wigner_sparse(interaction, n_qubits=6) 420 # n_qubits should default to 6 421 interaction_restrict = jw_number_restrict_operator( 422 interaction_sparse, 2) 423 424 dim = 6 * 5 // 2 # shape of new sparse array 425 426 # 3^ 2^ 4 1 maps 2**4 + 2 = 18 to 2**3 + 2**2 = 12 and vice versa; 427 # in the 2-particle subspace (1, 4) and (2, 3) are 7th and 9th. 428 expected = csc_matrix(([-1, -1], ([7, 9], [9, 7])), shape=(dim, dim)) 429 430 self.assertTrue(numpy.allclose(interaction_restrict.A, expected.A)) 431 432 def test_jw_restrict_jellium_ground_state_integration(self): 433 n_qubits = 4 434 grid = Grid(dimensions=1, length=n_qubits, scale=1.0) 435 jellium_hamiltonian = jordan_wigner_sparse( 436 jellium_model(grid, spinless=False)) 437 438 # 2 * n_qubits because of spin 439 number_sparse = jordan_wigner_sparse(number_operator(2 * n_qubits)) 440 441 restricted_number = jw_number_restrict_operator(number_sparse, 2) 442 restricted_jellium_hamiltonian = jw_number_restrict_operator( 443 jellium_hamiltonian, 2) 444 445 _, ground_state = get_ground_state(restricted_jellium_hamiltonian) 446 447 number_expectation = expectation(restricted_number, ground_state) 448 self.assertAlmostEqual(number_expectation, 2) 449 450 451class JWSzRestrictOperatorTest(unittest.TestCase): 452 453 def test_restrict_interaction_hamiltonian(self): 454 """Test restricting a coulomb repulsion Hamiltonian to a specified 455 Sz manifold.""" 456 x_dim = 3 457 y_dim = 2 458 459 interaction_term = fermi_hubbard(x_dim, y_dim, 0., 1.) 460 interaction_sparse = get_sparse_operator(interaction_term) 461 sz_value = 2 462 interaction_restricted = jw_sz_restrict_operator( 463 interaction_sparse, sz_value) 464 restricted_interaction_values = set( 465 [int(value.real) for value in interaction_restricted.diagonal()]) 466 # Originally the eigenvalues run from 0 to 6 but after restricting, 467 # they should run from 0 to 2 468 self.assertEqual(restricted_interaction_values, {0, 1, 2}) 469 470 471class JWNumberRestrictStateTest(unittest.TestCase): 472 473 def test_jw_number_restrict_state(self): 474 n_qubits = numpy.random.randint(1, 12) 475 n_particles = numpy.random.randint(0, n_qubits) 476 477 number_indices = jw_number_indices(n_particles, n_qubits) 478 subspace_dimension = len(number_indices) 479 480 # Create a vector that has entry 1 for every coordinate with 481 # the specified particle number, and 0 everywhere else 482 vector = numpy.zeros(2**n_qubits, dtype=float) 483 vector[number_indices] = 1 484 485 # Restrict the vector 486 restricted_vector = jw_number_restrict_state(vector, n_particles) 487 488 # Check that it has the correct shape 489 self.assertEqual(restricted_vector.shape[0], subspace_dimension) 490 491 # Check that it has the same norm as the original vector 492 self.assertAlmostEqual( 493 inner_product(vector, vector), 494 inner_product(restricted_vector, restricted_vector)) 495 496 497class JWSzRestrictStateTest(unittest.TestCase): 498 499 def test_jw_sz_restrict_state(self): 500 n_sites = numpy.random.randint(1, 10) 501 n_qubits = 2 * n_sites 502 sz_int = ((-1)**numpy.random.randint(2) * 503 numpy.random.randint(n_sites + 1)) 504 sz_value = sz_int / 2 505 506 sz_indices = jw_sz_indices(sz_value, n_qubits) 507 subspace_dimension = len(sz_indices) 508 509 # Create a vector that has entry 1 for every coordinate in 510 # the specified subspace, and 0 everywhere else 511 vector = numpy.zeros(2**n_qubits, dtype=float) 512 vector[sz_indices] = 1 513 514 # Restrict the vector 515 restricted_vector = jw_sz_restrict_state(vector, sz_value) 516 517 # Check that it has the correct shape 518 self.assertEqual(restricted_vector.shape[0], subspace_dimension) 519 520 # Check that it has the same norm as the original vector 521 self.assertAlmostEqual( 522 inner_product(vector, vector), 523 inner_product(restricted_vector, restricted_vector)) 524 525 526class JWGetGroundStatesByParticleNumberTest(unittest.TestCase): 527 528 def test_jw_get_ground_state_at_particle_number_herm_conserving(self): 529 # Initialize a particle-number-conserving Hermitian operator 530 ferm_op = FermionOperator('0^ 1') + FermionOperator('1^ 0') + \ 531 FermionOperator('1^ 2') + FermionOperator('2^ 1') + \ 532 FermionOperator('1^ 3', -.4) + FermionOperator('3^ 1', -.4) 533 jw_hamiltonian = jordan_wigner(ferm_op) 534 sparse_operator = get_sparse_operator(jw_hamiltonian) 535 n_qubits = 4 536 537 num_op = get_sparse_operator(number_operator(n_qubits)) 538 539 # Test each possible particle number 540 for particle_number in range(n_qubits): 541 # Get the ground energy and ground state at this particle number 542 energy, state = jw_get_ground_state_at_particle_number( 543 sparse_operator, particle_number) 544 545 # Check that it's an eigenvector with the correct eigenvalue 546 self.assertTrue( 547 numpy.allclose(sparse_operator.dot(state), energy * state)) 548 549 # Check that it has the correct particle number 550 num = expectation(num_op, state) 551 self.assertAlmostEqual(num, particle_number) 552 553 def test_jw_get_ground_state_at_particle_number_hubbard(self): 554 555 model = fermi_hubbard(2, 2, 1.0, 4.0) 556 sparse_operator = get_sparse_operator(model) 557 n_qubits = count_qubits(model) 558 num_op = get_sparse_operator(number_operator(n_qubits)) 559 560 # Test each possible particle number 561 for particle_number in range(n_qubits): 562 # Get the ground energy and ground state at this particle number 563 energy, state = jw_get_ground_state_at_particle_number( 564 sparse_operator, particle_number) 565 566 # Check that it's an eigenvector with the correct eigenvalue 567 self.assertTrue( 568 numpy.allclose(sparse_operator.dot(state), energy * state)) 569 570 # Check that it has the correct particle number 571 num = expectation(num_op, state) 572 self.assertAlmostEqual(num, particle_number) 573 574 def test_jw_get_ground_state_at_particle_number_jellium(self): 575 576 grid = Grid(2, 2, 1.0) 577 model = jellium_model(grid, spinless=True, plane_wave=False) 578 sparse_operator = get_sparse_operator(model) 579 n_qubits = count_qubits(model) 580 num_op = get_sparse_operator(number_operator(n_qubits)) 581 582 # Test each possible particle number 583 for particle_number in range(n_qubits): 584 # Get the ground energy and ground state at this particle number 585 energy, state = jw_get_ground_state_at_particle_number( 586 sparse_operator, particle_number) 587 588 # Check that it's an eigenvector with the correct eigenvalue 589 self.assertTrue( 590 numpy.allclose(sparse_operator.dot(state), energy * state)) 591 592 # Check that it has the correct particle number 593 num = expectation(num_op, state) 594 self.assertAlmostEqual(num, particle_number) 595 596 597class JWSparseGivensRotationTest(unittest.TestCase): 598 599 def test_bad_input(self): 600 with self.assertRaises(ValueError): 601 jw_sparse_givens_rotation(0, 2, 1., 1., 5) 602 with self.assertRaises(ValueError): 603 jw_sparse_givens_rotation(4, 5, 1., 1., 5) 604 605 606class GroundStateTest(unittest.TestCase): 607 608 def test_get_ground_state_hermitian(self): 609 ground = get_ground_state( 610 get_sparse_operator( 611 QubitOperator('Y0 X1') + QubitOperator('Z0 Z1'))) 612 expected_state = csc_matrix(([1j, 1], ([1, 2], [0, 0])), 613 shape=(4, 1), 614 dtype=numpy.complex128).A 615 expected_state /= numpy.sqrt(2.0) 616 617 self.assertAlmostEqual(ground[0], -2) 618 self.assertAlmostEqual( 619 numpy.absolute(expected_state.T.conj().dot(ground[1]))[0], 1.) 620 621 622class ExpectationTest(unittest.TestCase): 623 624 def test_expectation_correct_sparse_matrix(self): 625 operator = get_sparse_operator(QubitOperator('X0'), n_qubits=2) 626 vector = numpy.array([0., 1.j, 0., 1.j]) 627 self.assertAlmostEqual(expectation(operator, vector), 2.0) 628 629 density_matrix = scipy.sparse.csc_matrix( 630 numpy.outer(vector, numpy.conjugate(vector))) 631 self.assertAlmostEqual(expectation(operator, density_matrix), 2.0) 632 633 def test_expectation_correct_linear_operator(self): 634 operator = LinearQubitOperator(QubitOperator('X0'), n_qubits=2) 635 vector = numpy.array([0., 1.j, 0., 1.j]) 636 self.assertAlmostEqual(expectation(operator, vector), 2.0) 637 638 def test_expectation_handles_column_vector(self): 639 operator = get_sparse_operator(QubitOperator('X0'), n_qubits=2) 640 vector = numpy.array([[0.], [1.j], [0.], [1.j]]) 641 self.assertAlmostEqual(expectation(operator, vector), 2.0) 642 643 def test_expectation_correct_zero(self): 644 operator = get_sparse_operator(QubitOperator('X0'), n_qubits=2) 645 vector = numpy.array([1j, -1j, -1j, -1j]) 646 self.assertAlmostEqual(expectation(operator, vector), 0.0) 647 648 def test_execptions(self): 649 operator = LinearQubitOperator(QubitOperator('X0'), n_qubits=2) 650 vector = scipy.sparse.csc_matrix(numpy.array([1, 0, 0, 0])) 651 with self.assertRaises(ValueError): 652 _ = expectation(operator, vector) 653 with self.assertRaises(ValueError): 654 _ = expectation(vector, operator) 655 656 657class VarianceTest(unittest.TestCase): 658 659 def test_variance_row_vector(self): 660 661 X = pauli_matrix_map['X'] 662 Z = pauli_matrix_map['Z'] 663 zero = numpy.array([1., 0.]) 664 plus = numpy.array([1., 1.]) / numpy.sqrt(2) 665 minus = numpy.array([1., -1.]) / numpy.sqrt(2) 666 667 self.assertAlmostEqual(variance(Z, zero), 0.) 668 self.assertAlmostEqual(variance(X, zero), 1.) 669 670 self.assertAlmostEqual(variance(Z, plus), 1.) 671 self.assertAlmostEqual(variance(X, plus), 0.) 672 673 self.assertAlmostEqual(variance(Z, minus), 1.) 674 self.assertAlmostEqual(variance(X, minus), 0.) 675 676 def test_variance_column_vector(self): 677 678 X = pauli_matrix_map['X'] 679 Z = pauli_matrix_map['Z'] 680 zero = numpy.array([[1.], [0.]]) 681 plus = numpy.array([[1.], [1.]]) / numpy.sqrt(2) 682 minus = numpy.array([[1.], [-1.]]) / numpy.sqrt(2) 683 684 self.assertAlmostEqual(variance(Z, zero), 0.) 685 self.assertAlmostEqual(variance(X, zero), 1.) 686 687 self.assertAlmostEqual(variance(Z, plus), 1.) 688 self.assertAlmostEqual(variance(X, plus), 0.) 689 690 self.assertAlmostEqual(variance(Z, minus), 1.) 691 self.assertAlmostEqual(variance(X, minus), 0.) 692 693 694class ExpectationComputationalBasisStateTest(unittest.TestCase): 695 696 def test_expectation_fermion_operator_single_number_terms(self): 697 operator = FermionOperator('3^ 3', 1.9) + FermionOperator('2^ 1') 698 state = csc_matrix(([1], ([15], [0])), shape=(16, 1)) 699 700 self.assertAlmostEqual( 701 expectation_computational_basis_state(operator, state), 1.9) 702 703 def test_expectation_fermion_operator_two_number_terms(self): 704 operator = (FermionOperator('2^ 2', 1.9) + FermionOperator('2^ 1') + 705 FermionOperator('2^ 1^ 2 1', -1.7)) 706 state = csc_matrix(([1], ([6], [0])), shape=(16, 1)) 707 708 self.assertAlmostEqual( 709 expectation_computational_basis_state(operator, state), 3.6) 710 711 def test_expectation_identity_fermion_operator(self): 712 operator = FermionOperator.identity() * 1.1 713 state = csc_matrix(([1], ([6], [0])), shape=(16, 1)) 714 715 self.assertAlmostEqual( 716 expectation_computational_basis_state(operator, state), 1.1) 717 718 def test_expectation_state_is_list_single_number_terms(self): 719 operator = FermionOperator('3^ 3', 1.9) + FermionOperator('2^ 1') 720 state = [1, 1, 1, 1] 721 722 self.assertAlmostEqual( 723 expectation_computational_basis_state(operator, state), 1.9) 724 725 def test_expectation_state_is_list_fermion_operator_two_number_terms(self): 726 operator = (FermionOperator('2^ 2', 1.9) + FermionOperator('2^ 1') + 727 FermionOperator('2^ 1^ 2 1', -1.7)) 728 state = [0, 1, 1] 729 730 self.assertAlmostEqual( 731 expectation_computational_basis_state(operator, state), 3.6) 732 733 def test_expectation_state_is_list_identity_fermion_operator(self): 734 operator = FermionOperator.identity() * 1.1 735 state = [0, 1, 1] 736 737 self.assertAlmostEqual( 738 expectation_computational_basis_state(operator, state), 1.1) 739 740 def test_expectation_bad_operator_type(self): 741 with self.assertRaises(TypeError): 742 expectation_computational_basis_state( 743 'never', csc_matrix(([1], ([6], [0])), shape=(16, 1))) 744 745 def test_expectation_qubit_operator_not_implemented(self): 746 with self.assertRaises(NotImplementedError): 747 expectation_computational_basis_state( 748 QubitOperator(), csc_matrix(([1], ([6], [0])), shape=(16, 1))) 749 750 751class ExpectationDualBasisOperatorWithPlaneWaveBasisState(unittest.TestCase): 752 753 def setUp(self): 754 grid_length = 4 755 dimension = 1 756 wigner_seitz_radius = 10. 757 self.spinless = True 758 self.n_spatial_orbitals = grid_length**dimension 759 760 n_qubits = self.n_spatial_orbitals 761 self.n_particles = 3 762 763 # Compute appropriate length scale and the corresponding grid. 764 length_scale = wigner_seitz_length_scale(wigner_seitz_radius, 765 self.n_particles, dimension) 766 767 self.grid1 = Grid(dimension, grid_length, length_scale) 768 # Get the occupied orbitals of the plane-wave basis Hartree-Fock state. 769 hamiltonian = jellium_model(self.grid1, self.spinless, plane_wave=True) 770 hamiltonian = normal_ordered(hamiltonian) 771 hamiltonian.compress() 772 773 occupied_states = numpy.array( 774 lowest_single_particle_energy_states(hamiltonian, self.n_particles)) 775 self.hf_state_index1 = numpy.sum(2**occupied_states) 776 777 self.hf_state1 = numpy.zeros(2**n_qubits) 778 self.hf_state1[self.hf_state_index1] = 1.0 779 780 self.orbital_occupations1 = [ 781 digit == '1' for digit in bin(self.hf_state_index1)[2:] 782 ][::-1] 783 self.occupied_orbitals1 = [ 784 index for index, occupied in enumerate(self.orbital_occupations1) 785 if occupied 786 ] 787 788 self.reversed_occupied_orbitals1 = list(self.occupied_orbitals1) 789 for i in range(len(self.reversed_occupied_orbitals1)): 790 self.reversed_occupied_orbitals1[i] = -1 + int( 791 numpy.log2(self.hf_state1.shape[0]) 792 ) - self.reversed_occupied_orbitals1[i] 793 794 self.reversed_hf_state_index1 = sum( 795 2**index for index in self.reversed_occupied_orbitals1) 796 797 def test_1body_hopping_operator_1D(self): 798 operator = FermionOperator('2^ 0') 799 operator = normal_ordered(operator) 800 transformed_operator = normal_ordered( 801 fourier_transform(operator, self.grid1, self.spinless)) 802 803 expected = expectation(get_sparse_operator(transformed_operator), 804 self.hf_state1) 805 actual = expectation_db_operator_with_pw_basis_state( 806 operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals, 807 self.grid1, self.spinless) 808 self.assertAlmostEqual(expected, actual) 809 810 def test_1body_number_operator_1D(self): 811 operator = FermionOperator('2^ 2') 812 operator = normal_ordered(operator) 813 transformed_operator = normal_ordered( 814 fourier_transform(operator, self.grid1, self.spinless)) 815 816 expected = expectation(get_sparse_operator(transformed_operator), 817 self.hf_state1) 818 actual = expectation_db_operator_with_pw_basis_state( 819 operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals, 820 self.grid1, self.spinless) 821 self.assertAlmostEqual(expected, actual) 822 823 def test_2body_partial_number_operator_high_1D(self): 824 operator = FermionOperator('2^ 1^ 2 0') 825 operator = normal_ordered(operator) 826 transformed_operator = normal_ordered( 827 fourier_transform(operator, self.grid1, self.spinless)) 828 829 expected = expectation(get_sparse_operator(transformed_operator), 830 self.hf_state1) 831 actual = expectation_db_operator_with_pw_basis_state( 832 operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals, 833 self.grid1, self.spinless) 834 self.assertAlmostEqual(expected, actual) 835 836 def test_2body_partial_number_operator_mid_1D(self): 837 operator = FermionOperator('1^ 0^ 1 2') 838 operator = normal_ordered(operator) 839 transformed_operator = normal_ordered( 840 fourier_transform(operator, self.grid1, self.spinless)) 841 842 expected = expectation(get_sparse_operator(transformed_operator), 843 self.hf_state1) 844 actual = expectation_db_operator_with_pw_basis_state( 845 operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals, 846 self.grid1, self.spinless) 847 self.assertAlmostEqual(expected, actual) 848 849 def test_3body_double_number_operator_1D(self): 850 operator = FermionOperator('3^ 2^ 1^ 3 1 0') 851 operator = normal_ordered(operator) 852 transformed_operator = normal_ordered( 853 fourier_transform(operator, self.grid1, self.spinless)) 854 855 expected = expectation(get_sparse_operator(transformed_operator), 856 self.hf_state1) 857 actual = expectation_db_operator_with_pw_basis_state( 858 operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals, 859 self.grid1, self.spinless) 860 self.assertAlmostEqual(expected, actual) 861 862 def test_2body_adjacent_number_operator_1D(self): 863 operator = FermionOperator('3^ 2^ 2 1') 864 operator = normal_ordered(operator) 865 transformed_operator = normal_ordered( 866 fourier_transform(operator, self.grid1, self.spinless)) 867 868 expected = expectation(get_sparse_operator(transformed_operator), 869 self.hf_state1) 870 actual = expectation_db_operator_with_pw_basis_state( 871 operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals, 872 self.grid1, self.spinless) 873 self.assertAlmostEqual(expected, actual) 874 875 def test_1d5_with_spin_10particles(self): 876 dimension = 1 877 grid_length = 5 878 n_spatial_orbitals = grid_length**dimension 879 wigner_seitz_radius = 9.3 880 881 spinless = False 882 n_qubits = n_spatial_orbitals 883 if not spinless: 884 n_qubits *= 2 885 n_particles_big = 10 886 887 length_scale = wigner_seitz_length_scale(wigner_seitz_radius, 888 n_particles_big, dimension) 889 890 self.grid3 = Grid(dimension, grid_length, length_scale) 891 # Get the occupied orbitals of the plane-wave basis Hartree-Fock state. 892 hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True) 893 hamiltonian = normal_ordered(hamiltonian) 894 hamiltonian.compress() 895 896 occupied_states = numpy.array( 897 lowest_single_particle_energy_states(hamiltonian, n_particles_big)) 898 self.hf_state_index3 = numpy.sum(2**occupied_states) 899 900 self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])), 901 shape=(2**n_qubits, 1)) 902 903 self.orbital_occupations3 = [ 904 digit == '1' for digit in bin(self.hf_state_index3)[2:] 905 ][::-1] 906 self.occupied_orbitals3 = [ 907 index for index, occupied in enumerate(self.orbital_occupations3) 908 if occupied 909 ] 910 911 self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3) 912 for i in range(len(self.reversed_occupied_orbitals3)): 913 self.reversed_occupied_orbitals3[i] = -1 + int( 914 numpy.log2(self.hf_state3.shape[0]) 915 ) - self.reversed_occupied_orbitals3[i] 916 917 self.reversed_hf_state_index3 = sum( 918 2**index for index in self.reversed_occupied_orbitals3) 919 920 operator = (FermionOperator('6^ 0^ 1^ 3 5 4', 2) + 921 FermionOperator('7^ 6^ 5 4', -3.7j) + 922 FermionOperator('3^ 3', 2.1) + FermionOperator('3^ 2', 1.7)) 923 operator = normal_ordered(operator) 924 normal_ordered(fourier_transform(operator, self.grid3, spinless)) 925 926 expected = 2.1 927 # Calculated from expectation(get_sparse_operator( 928 # transformed_operator), self.hf_state3) 929 actual = expectation_db_operator_with_pw_basis_state( 930 operator, self.reversed_occupied_orbitals3, n_spatial_orbitals, 931 self.grid3, spinless) 932 933 self.assertAlmostEqual(expected, actual) 934 935 def test_1d5_with_spin_7particles(self): 936 dimension = 1 937 grid_length = 5 938 n_spatial_orbitals = grid_length**dimension 939 wigner_seitz_radius = 9.3 940 941 spinless = False 942 n_qubits = n_spatial_orbitals 943 if not spinless: 944 n_qubits *= 2 945 n_particles_big = 7 946 947 length_scale = wigner_seitz_length_scale(wigner_seitz_radius, 948 n_particles_big, dimension) 949 950 self.grid3 = Grid(dimension, grid_length, length_scale) 951 # Get the occupied orbitals of the plane-wave basis Hartree-Fock state. 952 hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True) 953 hamiltonian = normal_ordered(hamiltonian) 954 hamiltonian.compress() 955 956 occupied_states = numpy.array( 957 lowest_single_particle_energy_states(hamiltonian, n_particles_big)) 958 self.hf_state_index3 = numpy.sum(2**occupied_states) 959 960 self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])), 961 shape=(2**n_qubits, 1)) 962 963 self.orbital_occupations3 = [ 964 digit == '1' for digit in bin(self.hf_state_index3)[2:] 965 ][::-1] 966 self.occupied_orbitals3 = [ 967 index for index, occupied in enumerate(self.orbital_occupations3) 968 if occupied 969 ] 970 971 self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3) 972 for i in range(len(self.reversed_occupied_orbitals3)): 973 self.reversed_occupied_orbitals3[i] = -1 + int( 974 numpy.log2(self.hf_state3.shape[0]) 975 ) - self.reversed_occupied_orbitals3[i] 976 977 self.reversed_hf_state_index3 = sum( 978 2**index for index in self.reversed_occupied_orbitals3) 979 980 operator = (FermionOperator('6^ 0^ 1^ 3 5 4', 2) + 981 FermionOperator('7^ 2^ 4 1') + 982 FermionOperator('3^ 3', 2.1) + 983 FermionOperator('5^ 3^ 1 0', 7.3)) 984 operator = normal_ordered(operator) 985 normal_ordered(fourier_transform(operator, self.grid3, spinless)) 986 987 expected = 1.66 - 0.0615536707435j 988 # Calculated with expected = expectation(get_sparse_operator( 989 # transformed_operator), self.hf_state3) 990 actual = expectation_db_operator_with_pw_basis_state( 991 operator, self.reversed_occupied_orbitals3, n_spatial_orbitals, 992 self.grid3, spinless) 993 994 self.assertAlmostEqual(expected, actual) 995 996 def test_3d2_spinless(self): 997 dimension = 3 998 grid_length = 2 999 n_spatial_orbitals = grid_length**dimension 1000 wigner_seitz_radius = 9.3 1001 1002 spinless = True 1003 n_qubits = n_spatial_orbitals 1004 n_particles_big = 5 1005 1006 length_scale = wigner_seitz_length_scale(wigner_seitz_radius, 1007 n_particles_big, dimension) 1008 1009 self.grid3 = Grid(dimension, grid_length, length_scale) 1010 # Get the occupied orbitals of the plane-wave basis Hartree-Fock state. 1011 hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True) 1012 hamiltonian = normal_ordered(hamiltonian) 1013 hamiltonian.compress() 1014 1015 occupied_states = numpy.array( 1016 lowest_single_particle_energy_states(hamiltonian, n_particles_big)) 1017 self.hf_state_index3 = numpy.sum(2**occupied_states) 1018 1019 self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])), 1020 shape=(2**n_qubits, 1)) 1021 1022 self.orbital_occupations3 = [ 1023 digit == '1' for digit in bin(self.hf_state_index3)[2:] 1024 ][::-1] 1025 self.occupied_orbitals3 = [ 1026 index for index, occupied in enumerate(self.orbital_occupations3) 1027 if occupied 1028 ] 1029 1030 self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3) 1031 for i in range(len(self.reversed_occupied_orbitals3)): 1032 self.reversed_occupied_orbitals3[i] = -1 + int( 1033 numpy.log2(self.hf_state3.shape[0]) 1034 ) - self.reversed_occupied_orbitals3[i] 1035 1036 self.reversed_hf_state_index3 = sum( 1037 2**index for index in self.reversed_occupied_orbitals3) 1038 1039 operator = (FermionOperator('4^ 2^ 3^ 5 5 4', 2) + 1040 FermionOperator('7^ 6^ 7 4', -3.7j) + 1041 FermionOperator('3^ 7', 2.1)) 1042 operator = normal_ordered(operator) 1043 normal_ordered(fourier_transform(operator, self.grid3, spinless)) 1044 expected = -0.2625 - 0.4625j 1045 # Calculated with expectation(get_sparse_operator( 1046 # transformed_operator), self.hf_state3) 1047 actual = expectation_db_operator_with_pw_basis_state( 1048 operator, self.reversed_occupied_orbitals3, n_spatial_orbitals, 1049 self.grid3, spinless) 1050 1051 self.assertAlmostEqual(expected, actual) 1052 1053 def test_3d2_with_spin(self): 1054 dimension = 3 1055 grid_length = 2 1056 n_spatial_orbitals = grid_length**dimension 1057 wigner_seitz_radius = 9.3 1058 1059 spinless = False 1060 n_qubits = n_spatial_orbitals 1061 if not spinless: 1062 n_qubits *= 2 1063 n_particles_big = 9 1064 1065 length_scale = wigner_seitz_length_scale(wigner_seitz_radius, 1066 n_particles_big, dimension) 1067 1068 self.grid3 = Grid(dimension, grid_length, length_scale) 1069 # Get the occupied orbitals of the plane-wave basis Hartree-Fock state. 1070 hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True) 1071 hamiltonian = normal_ordered(hamiltonian) 1072 hamiltonian.compress() 1073 1074 occupied_states = numpy.array( 1075 lowest_single_particle_energy_states(hamiltonian, n_particles_big)) 1076 self.hf_state_index3 = numpy.sum(2**occupied_states) 1077 1078 self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])), 1079 shape=(2**n_qubits, 1)) 1080 1081 self.orbital_occupations3 = [ 1082 digit == '1' for digit in bin(self.hf_state_index3)[2:] 1083 ][::-1] 1084 self.occupied_orbitals3 = [ 1085 index for index, occupied in enumerate(self.orbital_occupations3) 1086 if occupied 1087 ] 1088 1089 self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3) 1090 for i in range(len(self.reversed_occupied_orbitals3)): 1091 self.reversed_occupied_orbitals3[i] = -1 + int( 1092 numpy.log2(self.hf_state3.shape[0]) 1093 ) - self.reversed_occupied_orbitals3[i] 1094 1095 self.reversed_hf_state_index3 = sum( 1096 2**index for index in self.reversed_occupied_orbitals3) 1097 1098 operator = (FermionOperator('4^ 2^ 3^ 5 5 4', 2) + 1099 FermionOperator('7^ 6^ 7 4', -3.7j) + 1100 FermionOperator('3^ 7', 2.1)) 1101 operator = normal_ordered(operator) 1102 normal_ordered(fourier_transform(operator, self.grid3, spinless)) 1103 1104 expected = -0.2625 - 0.578125j 1105 # Calculated from expected = expectation(get_sparse_operator( 1106 # transformed_operator), self.hf_state3) 1107 actual = expectation_db_operator_with_pw_basis_state( 1108 operator, self.reversed_occupied_orbitals3, n_spatial_orbitals, 1109 self.grid3, spinless) 1110 1111 self.assertAlmostEqual(expected, actual) 1112 1113 1114class GetGapTest(unittest.TestCase): 1115 1116 def test_get_gap(self): 1117 operator = QubitOperator('Y0 X1') + QubitOperator('Z0 Z1') 1118 self.assertAlmostEqual(get_gap(get_sparse_operator(operator)), 2.0) 1119 1120 def test_get_gap_nonhermitian_error(self): 1121 operator = (QubitOperator('X0 Y1', 1 + 1j) + 1122 QubitOperator('Z0 Z1', 1j) + QubitOperator((), 2 + 1j)) 1123 with self.assertRaises(ValueError): 1124 get_gap(get_sparse_operator(operator)) 1125 1126 1127class InnerProductTest(unittest.TestCase): 1128 1129 def test_inner_product(self): 1130 state_1 = numpy.array([1., 1.j]) 1131 state_2 = numpy.array([1., -1.j]) 1132 1133 self.assertAlmostEqual(inner_product(state_1, state_1), 2.) 1134 self.assertAlmostEqual(inner_product(state_1, state_2), 0.) 1135 1136 1137class BosonSparseTest(unittest.TestCase): 1138 1139 def setUp(self): 1140 self.hbar = 1. 1141 self.d = 5 1142 self.b = numpy.diag(numpy.sqrt(numpy.arange(1, self.d)), 1) 1143 self.bd = self.b.conj().T 1144 self.q = numpy.sqrt(self.hbar / 2) * (self.b + self.bd) 1145 self.p = -1j * numpy.sqrt(self.hbar / 2) * (self.b - self.bd) 1146 self.Id = numpy.identity(self.d) 1147 1148 def test_boson_ladder_noninteger_trunc(self): 1149 with self.assertRaises(ValueError): 1150 boson_ladder_sparse(1, 0, 0, 0.1) 1151 1152 with self.assertRaises(ValueError): 1153 boson_ladder_sparse(1, 0, 0, -1) 1154 1155 with self.assertRaises(ValueError): 1156 boson_ladder_sparse(1, 0, 0, 0) 1157 1158 def test_boson_ladder_destroy_one_mode(self): 1159 b = boson_ladder_sparse(1, 0, 0, self.d).toarray() 1160 self.assertTrue(numpy.allclose(b, self.b)) 1161 1162 def test_boson_ladder_create_one_mode(self): 1163 bd = boson_ladder_sparse(1, 0, 1, self.d).toarray() 1164 self.assertTrue(numpy.allclose(bd, self.bd)) 1165 1166 def test_boson_ladder_single_adjoint(self): 1167 b = boson_ladder_sparse(1, 0, 0, self.d).toarray() 1168 bd = boson_ladder_sparse(1, 0, 1, self.d).toarray() 1169 self.assertTrue(numpy.allclose(b.conj().T, bd)) 1170 1171 def test_boson_ladder_two_mode(self): 1172 res = boson_ladder_sparse(2, 0, 0, self.d).toarray() 1173 expected = numpy.kron(self.b, self.Id) 1174 self.assertTrue(numpy.allclose(res, expected)) 1175 1176 res = boson_ladder_sparse(2, 1, 0, self.d).toarray() 1177 expected = numpy.kron(self.Id, self.b) 1178 self.assertTrue(numpy.allclose(res, expected)) 1179 1180 def test_single_quad_noninteger_trunc(self): 1181 with self.assertRaises(ValueError): 1182 single_quad_op_sparse(1, 0, 'q', self.hbar, 0.1) 1183 1184 with self.assertRaises(ValueError): 1185 single_quad_op_sparse(1, 0, 'q', self.hbar, -1) 1186 1187 with self.assertRaises(ValueError): 1188 single_quad_op_sparse(1, 0, 'q', self.hbar, 0) 1189 1190 def test_single_quad_q_one_mode(self): 1191 res = single_quad_op_sparse(1, 0, 'q', self.hbar, self.d).toarray() 1192 self.assertTrue(numpy.allclose(res, self.q)) 1193 self.assertTrue(numpy.allclose(res, res.conj().T)) 1194 1195 def test_single_quad_p_one_mode(self): 1196 res = single_quad_op_sparse(1, 0, 'p', self.hbar, self.d).toarray() 1197 self.assertTrue(numpy.allclose(res, self.p)) 1198 self.assertTrue(numpy.allclose(res, res.conj().T)) 1199 1200 def test_single_quad_two_mode(self): 1201 res = single_quad_op_sparse(2, 0, 'q', self.hbar, self.d).toarray() 1202 expected = numpy.kron(self.q, self.Id) 1203 self.assertTrue(numpy.allclose(res, expected)) 1204 1205 res = single_quad_op_sparse(2, 1, 'p', self.hbar, self.d).toarray() 1206 expected = numpy.kron(self.Id, self.p) 1207 self.assertTrue(numpy.allclose(res, expected)) 1208 1209 def test_boson_operator_sparse_trunc(self): 1210 op = BosonOperator('0') 1211 with self.assertRaises(ValueError): 1212 boson_operator_sparse(op, 0.1) 1213 1214 with self.assertRaises(ValueError): 1215 boson_operator_sparse(op, -1) 1216 1217 with self.assertRaises(ValueError): 1218 boson_operator_sparse(op, 0) 1219 1220 def test_boson_operator_invalid_op(self): 1221 op = FermionOperator('0') 1222 with self.assertRaises(ValueError): 1223 boson_operator_sparse(op, self.d) 1224 1225 def test_boson_operator_sparse_empty(self): 1226 for op in (BosonOperator(), QuadOperator()): 1227 res = boson_operator_sparse(op, self.d) 1228 self.assertEqual(res, numpy.array([[0]])) 1229 1230 def test_boson_operator_sparse_identity(self): 1231 for op in (BosonOperator(''), QuadOperator('')): 1232 res = boson_operator_sparse(op, self.d) 1233 self.assertEqual(res, numpy.array([[1]])) 1234 1235 def test_boson_operator_sparse_single(self): 1236 op = BosonOperator('0') 1237 res = boson_operator_sparse(op, self.d).toarray() 1238 self.assertTrue(numpy.allclose(res, self.b)) 1239 1240 op = BosonOperator('0^') 1241 res = boson_operator_sparse(op, self.d).toarray() 1242 self.assertTrue(numpy.allclose(res, self.bd)) 1243 1244 op = QuadOperator('q0') 1245 res = boson_operator_sparse(op, self.d, self.hbar).toarray() 1246 self.assertTrue(numpy.allclose(res, self.q)) 1247 1248 op = QuadOperator('p0') 1249 res = boson_operator_sparse(op, self.d, self.hbar).toarray() 1250 self.assertTrue(numpy.allclose(res, self.p)) 1251 1252 def test_boson_operator_sparse_number(self): 1253 op = BosonOperator('0^ 0') 1254 res = boson_operator_sparse(op, self.d).toarray() 1255 self.assertTrue(numpy.allclose(res, numpy.dot(self.bd, self.b))) 1256 1257 def test_boson_operator_sparse_multi_mode(self): 1258 op = BosonOperator('0^ 1 1^ 2') 1259 res = boson_operator_sparse(op, self.d).toarray() 1260 1261 b0 = boson_ladder_sparse(3, 0, 0, self.d).toarray() 1262 b1 = boson_ladder_sparse(3, 1, 0, self.d).toarray() 1263 b2 = boson_ladder_sparse(3, 2, 0, self.d).toarray() 1264 1265 expected = multi_dot([b0.T, b1, b1.T, b2]) 1266 self.assertTrue(numpy.allclose(res, expected)) 1267 1268 op = QuadOperator('q0 p0 p1') 1269 res = boson_operator_sparse(op, self.d, self.hbar).toarray() 1270 1271 expected = numpy.identity(self.d**2) 1272 for term in op.terms: 1273 for i, j in term: 1274 expected = expected.dot( 1275 single_quad_op_sparse(2, i, j, self.hbar, self.d).toarray()) 1276 self.assertTrue(numpy.allclose(res, expected)) 1277 1278 def test_boson_operator_sparse_addition(self): 1279 op = BosonOperator('0^ 1') 1280 op += BosonOperator('0 0^') 1281 res = boson_operator_sparse(op, self.d).toarray() 1282 1283 b0 = boson_ladder_sparse(2, 0, 0, self.d).toarray() 1284 b1 = boson_ladder_sparse(2, 1, 0, self.d).toarray() 1285 1286 expected = numpy.dot(b0.T, b1) + numpy.dot(b0, b0.T) 1287 self.assertTrue(numpy.allclose(res, expected)) 1288 1289 1290class GetNumberPreservingSparseOperatorIntegrationTestLiH(unittest.TestCase): 1291 1292 def setUp(self): 1293 # Set up molecule. 1294 geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))] 1295 basis = 'sto-3g' 1296 multiplicity = 1 1297 filename = os.path.join(DATA_DIRECTORY, 'H1-Li1_sto-3g_singlet_1.45') 1298 self.molecule = MolecularData(geometry, 1299 basis, 1300 multiplicity, 1301 filename=filename) 1302 self.molecule.load() 1303 1304 # Get molecular Hamiltonian. 1305 self.molecular_hamiltonian = self.molecule.get_molecular_hamiltonian() 1306 1307 self.hubbard_hamiltonian = fermi_hubbard(2, 1308 2, 1309 1.0, 1310 4.0, 1311 chemical_potential=.2, 1312 magnetic_field=0.0, 1313 spinless=False) 1314 1315 def test_exceptions(self): 1316 op = FermionOperator('1') 1317 with self.assertRaises(ValueError): 1318 _ = get_number_preserving_sparse_operator(op, 2, 1) 1319 1320 def test_number_on_reference(self): 1321 sum_n_op = FermionOperator() 1322 sum_sparse_n_op = get_number_preserving_sparse_operator( 1323 sum_n_op, 1324 self.molecule.n_qubits, 1325 self.molecule.n_electrons, 1326 spin_preserving=False) 1327 1328 space_size = sum_sparse_n_op.shape[0] 1329 reference = numpy.zeros((space_size)) 1330 reference[0] = 1.0 1331 1332 for i in range(self.molecule.n_qubits): 1333 n_op = FermionOperator(((i, 1), (i, 0))) 1334 sum_n_op += n_op 1335 1336 sparse_n_op = get_number_preserving_sparse_operator( 1337 n_op, 1338 self.molecule.n_qubits, 1339 self.molecule.n_electrons, 1340 spin_preserving=False) 1341 1342 sum_sparse_n_op += sparse_n_op 1343 1344 expectation = reference.dot(sparse_n_op.dot(reference)) 1345 1346 if i < self.molecule.n_electrons: 1347 assert expectation == 1.0 1348 else: 1349 assert expectation == 0.0 1350 1351 convert_after_adding = get_number_preserving_sparse_operator( 1352 sum_n_op, 1353 self.molecule.n_qubits, 1354 self.molecule.n_electrons, 1355 spin_preserving=False) 1356 1357 assert scipy.sparse.linalg.norm(convert_after_adding - 1358 sum_sparse_n_op) < 1E-9 1359 1360 assert reference.dot(sum_sparse_n_op.dot(reference)) - \ 1361 self.molecule.n_electrons < 1E-9 1362 1363 def test_space_size_correct(self): 1364 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1365 1366 sparse_ham = get_number_preserving_sparse_operator( 1367 hamiltonian_fop, 1368 self.molecule.n_qubits, 1369 self.molecule.n_electrons, 1370 spin_preserving=True) 1371 1372 space_size = sparse_ham.shape[0] 1373 1374 # Naive Hilbert space size is 2**12, or 4096. 1375 assert space_size == 225 1376 1377 def test_hf_energy(self): 1378 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1379 1380 sparse_ham = get_number_preserving_sparse_operator( 1381 hamiltonian_fop, 1382 self.molecule.n_qubits, 1383 self.molecule.n_electrons, 1384 spin_preserving=True) 1385 1386 space_size = sparse_ham.shape[0] 1387 reference = numpy.zeros((space_size)) 1388 reference[0] = 1.0 1389 1390 sparse_hf_energy = reference.dot(sparse_ham.dot(reference)) 1391 1392 assert numpy.linalg.norm(sparse_hf_energy - 1393 self.molecule.hf_energy) < 1E-9 1394 1395 def test_one_body_hf_energy(self): 1396 one_body_part = self.molecular_hamiltonian 1397 one_body_part.two_body_tensor = numpy.zeros_like( 1398 one_body_part.two_body_tensor) 1399 1400 one_body_fop = get_fermion_operator(one_body_part) 1401 one_body_regular_sparse_op = get_sparse_operator(one_body_fop) 1402 1403 make_hf_fop = FermionOperator(((3, 1), (2, 1), (1, 1), (0, 1))) 1404 make_hf_sparse_op = get_sparse_operator(make_hf_fop, n_qubits=12) 1405 1406 hf_state = numpy.zeros((2**12)) 1407 hf_state[0] = 1.0 1408 hf_state = make_hf_sparse_op.dot(hf_state) 1409 1410 regular_sparse_hf_energy = \ 1411 (hf_state.dot(one_body_regular_sparse_op.dot(hf_state))).real 1412 1413 one_body_sparse_op = get_number_preserving_sparse_operator( 1414 one_body_fop, 1415 self.molecule.n_qubits, 1416 self.molecule.n_electrons, 1417 spin_preserving=True) 1418 1419 space_size = one_body_sparse_op.shape[0] 1420 reference = numpy.zeros((space_size)) 1421 reference[0] = 1.0 1422 1423 sparse_hf_energy = reference.dot(one_body_sparse_op.dot(reference)) 1424 1425 assert numpy.linalg.norm(sparse_hf_energy - 1426 regular_sparse_hf_energy) < 1E-9 1427 1428 def test_ground_state_energy(self): 1429 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1430 1431 sparse_ham = get_number_preserving_sparse_operator( 1432 hamiltonian_fop, 1433 self.molecule.n_qubits, 1434 self.molecule.n_electrons, 1435 spin_preserving=True) 1436 1437 eig_val, _ = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA') 1438 1439 assert numpy.abs(eig_val[0] - self.molecule.fci_energy) < 1E-9 1440 1441 def test_doubles_are_subset(self): 1442 reference_determinants = [[ 1443 True, True, True, True, False, False, False, False, False, False, 1444 False, False 1445 ], 1446 [ 1447 True, True, False, False, False, False, 1448 True, True, False, False, False, False 1449 ]] 1450 1451 for reference_determinant in reference_determinants: 1452 reference_determinant = numpy.asarray(reference_determinant) 1453 doubles_state_array = numpy.asarray( 1454 list( 1455 _iterate_basis_(reference_determinant, 1456 excitation_level=2, 1457 spin_preserving=True))) 1458 doubles_int_state_array = doubles_state_array.dot( 1459 1 << numpy.arange(doubles_state_array.shape[1])[::-1]) 1460 1461 all_state_array = numpy.asarray( 1462 list( 1463 _iterate_basis_(reference_determinant, 1464 excitation_level=4, 1465 spin_preserving=True))) 1466 all_int_state_array = all_state_array.dot( 1467 1 << numpy.arange(all_state_array.shape[1])[::-1]) 1468 1469 for item in doubles_int_state_array: 1470 assert item in all_int_state_array 1471 1472 for reference_determinant in reference_determinants: 1473 reference_determinant = numpy.asarray(reference_determinant) 1474 doubles_state_array = numpy.asarray( 1475 list( 1476 _iterate_basis_(reference_determinant, 1477 excitation_level=2, 1478 spin_preserving=True))) 1479 doubles_int_state_array = doubles_state_array.dot( 1480 1 << numpy.arange(doubles_state_array.shape[1])[::-1]) 1481 1482 all_state_array = numpy.asarray( 1483 list( 1484 _iterate_basis_(reference_determinant, 1485 excitation_level=4, 1486 spin_preserving=False))) 1487 all_int_state_array = all_state_array.dot( 1488 1 << numpy.arange(all_state_array.shape[1])[::-1]) 1489 1490 for item in doubles_int_state_array: 1491 assert item in all_int_state_array 1492 1493 for reference_determinant in reference_determinants: 1494 reference_determinant = numpy.asarray(reference_determinant) 1495 doubles_state_array = numpy.asarray( 1496 list( 1497 _iterate_basis_(reference_determinant, 1498 excitation_level=2, 1499 spin_preserving=False))) 1500 doubles_int_state_array = doubles_state_array.dot( 1501 1 << numpy.arange(doubles_state_array.shape[1])[::-1]) 1502 1503 all_state_array = numpy.asarray( 1504 list( 1505 _iterate_basis_(reference_determinant, 1506 excitation_level=4, 1507 spin_preserving=False))) 1508 all_int_state_array = all_state_array.dot( 1509 1 << numpy.arange(all_state_array.shape[1])[::-1]) 1510 1511 for item in doubles_int_state_array: 1512 assert item in all_int_state_array 1513 1514 def test_full_ham_hermitian(self): 1515 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1516 1517 sparse_ham = get_number_preserving_sparse_operator( 1518 hamiltonian_fop, 1519 self.molecule.n_qubits, 1520 self.molecule.n_electrons, 1521 spin_preserving=True) 1522 1523 assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9 1524 1525 def test_full_ham_hermitian_non_spin_preserving(self): 1526 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1527 1528 sparse_ham = get_number_preserving_sparse_operator( 1529 hamiltonian_fop, 1530 self.molecule.n_qubits, 1531 self.molecule.n_electrons, 1532 spin_preserving=False) 1533 1534 assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9 1535 1536 def test_singles_simple_one_body_term_hermitian(self): 1537 fop = FermionOperator(((3, 1), (1, 0))) 1538 fop_conj = FermionOperator(((1, 1), (3, 0))) 1539 1540 sparse_op = get_number_preserving_sparse_operator( 1541 fop, 1542 self.molecule.n_qubits, 1543 self.molecule.n_electrons, 1544 spin_preserving=True, 1545 excitation_level=1) 1546 1547 sparse_op_conj = get_number_preserving_sparse_operator( 1548 fop_conj, 1549 self.molecule.n_qubits, 1550 self.molecule.n_electrons, 1551 spin_preserving=True, 1552 excitation_level=1) 1553 1554 assert scipy.sparse.linalg.norm(sparse_op - 1555 sparse_op_conj.getH()) < 1E-9 1556 1557 def test_singles_simple_two_body_term_hermitian(self): 1558 fop = FermionOperator(((3, 1), (8, 1), (1, 0), (4, 0))) 1559 fop_conj = FermionOperator(((4, 1), (1, 1), (8, 0), (3, 0))) 1560 1561 sparse_op = get_number_preserving_sparse_operator( 1562 fop, 1563 self.molecule.n_qubits, 1564 self.molecule.n_electrons, 1565 spin_preserving=True, 1566 excitation_level=1) 1567 1568 sparse_op_conj = get_number_preserving_sparse_operator( 1569 fop_conj, 1570 self.molecule.n_qubits, 1571 self.molecule.n_electrons, 1572 spin_preserving=True, 1573 excitation_level=1) 1574 1575 assert scipy.sparse.linalg.norm(sparse_op - 1576 sparse_op_conj.getH()) < 1E-9 1577 1578 def test_singles_repeating_two_body_term_hermitian(self): 1579 fop = FermionOperator(((3, 1), (1, 1), (5, 0), (1, 0))) 1580 fop_conj = FermionOperator(((5, 1), (1, 1), (3, 0), (1, 0))) 1581 1582 sparse_op = get_number_preserving_sparse_operator( 1583 fop, 1584 self.molecule.n_qubits, 1585 self.molecule.n_electrons, 1586 spin_preserving=True, 1587 excitation_level=1) 1588 1589 sparse_op_conj = get_number_preserving_sparse_operator( 1590 fop_conj, 1591 self.molecule.n_qubits, 1592 self.molecule.n_electrons, 1593 spin_preserving=True, 1594 excitation_level=1) 1595 1596 assert scipy.sparse.linalg.norm(sparse_op - 1597 sparse_op_conj.getH()) < 1E-9 1598 1599 def test_singles_ham_hermitian(self): 1600 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1601 1602 sparse_ham = get_number_preserving_sparse_operator( 1603 hamiltonian_fop, 1604 self.molecule.n_qubits, 1605 self.molecule.n_electrons, 1606 spin_preserving=True, 1607 excitation_level=1) 1608 1609 assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9 1610 1611 def test_singles_ham_hermitian_non_spin_preserving(self): 1612 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1613 1614 sparse_ham = get_number_preserving_sparse_operator( 1615 hamiltonian_fop, 1616 self.molecule.n_qubits, 1617 self.molecule.n_electrons, 1618 spin_preserving=False, 1619 excitation_level=1) 1620 1621 assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9 1622 1623 def test_cisd_energy(self): 1624 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1625 1626 sparse_ham = get_number_preserving_sparse_operator( 1627 hamiltonian_fop, 1628 self.molecule.n_qubits, 1629 self.molecule.n_electrons, 1630 spin_preserving=True, 1631 excitation_level=2) 1632 1633 eig_val, _ = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA') 1634 1635 assert numpy.abs(eig_val[0] - self.molecule.cisd_energy) < 1E-9 1636 1637 def test_cisd_energy_non_spin_preserving(self): 1638 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1639 1640 sparse_ham = get_number_preserving_sparse_operator( 1641 hamiltonian_fop, 1642 self.molecule.n_qubits, 1643 self.molecule.n_electrons, 1644 spin_preserving=False, 1645 excitation_level=2) 1646 1647 eig_val, _ = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA') 1648 1649 assert numpy.abs(eig_val[0] - self.molecule.cisd_energy) < 1E-9 1650 1651 def test_cisd_matches_fci_energy_two_electron_hubbard(self): 1652 hamiltonian_fop = self.hubbard_hamiltonian 1653 1654 sparse_ham_cisd = get_number_preserving_sparse_operator( 1655 hamiltonian_fop, 8, 2, spin_preserving=True, excitation_level=2) 1656 1657 sparse_ham_fci = get_sparse_operator(hamiltonian_fop, n_qubits=8) 1658 1659 eig_val_cisd, _ = scipy.sparse.linalg.eigsh(sparse_ham_cisd, 1660 k=1, 1661 which='SA') 1662 eig_val_fci, _ = scipy.sparse.linalg.eigsh(sparse_ham_fci, 1663 k=1, 1664 which='SA') 1665 1666 assert numpy.abs(eig_val_cisd[0] - eig_val_fci[0]) < 1E-9 1667 1668 def test_weird_determinant_matches_fci_energy_two_electron_hubbard(self): 1669 hamiltonian_fop = self.hubbard_hamiltonian 1670 1671 sparse_ham_cisd = get_number_preserving_sparse_operator( 1672 hamiltonian_fop, 1673 8, 1674 2, 1675 spin_preserving=True, 1676 excitation_level=2, 1677 reference_determinant=numpy.asarray( 1678 [False, False, True, True, False, False, False, False])) 1679 1680 sparse_ham_fci = get_sparse_operator(hamiltonian_fop, n_qubits=8) 1681 1682 eig_val_cisd, _ = scipy.sparse.linalg.eigsh(sparse_ham_cisd, 1683 k=1, 1684 which='SA') 1685 eig_val_fci, _ = scipy.sparse.linalg.eigsh(sparse_ham_fci, 1686 k=1, 1687 which='SA') 1688 1689 assert numpy.abs(eig_val_cisd[0] - eig_val_fci[0]) < 1E-9 1690 1691 def test_number_restricted_spectra_match_molecule(self): 1692 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1693 1694 sparse_ham_number_preserving = get_number_preserving_sparse_operator( 1695 hamiltonian_fop, 1696 self.molecule.n_qubits, 1697 self.molecule.n_electrons, 1698 spin_preserving=False) 1699 1700 sparse_ham = get_sparse_operator(hamiltonian_fop, 1701 self.molecule.n_qubits) 1702 1703 sparse_ham_restricted_number_preserving = jw_number_restrict_operator( 1704 sparse_ham, 1705 n_electrons=self.molecule.n_electrons, 1706 n_qubits=self.molecule.n_qubits) 1707 1708 spectrum_from_new_sparse_method = sparse_eigenspectrum( 1709 sparse_ham_number_preserving) 1710 1711 spectrum_from_old_sparse_method = sparse_eigenspectrum( 1712 sparse_ham_restricted_number_preserving) 1713 1714 spectral_deviation = numpy.amax( 1715 numpy.absolute(spectrum_from_new_sparse_method - 1716 spectrum_from_old_sparse_method)) 1717 self.assertAlmostEqual(spectral_deviation, 0.) 1718 1719 def test_number_restricted_spectra_match_hubbard(self): 1720 hamiltonian_fop = self.hubbard_hamiltonian 1721 1722 sparse_ham_number_preserving = get_number_preserving_sparse_operator( 1723 hamiltonian_fop, 8, 4, spin_preserving=False) 1724 1725 sparse_ham = get_sparse_operator(hamiltonian_fop, 8) 1726 1727 sparse_ham_restricted_number_preserving = jw_number_restrict_operator( 1728 sparse_ham, n_electrons=4, n_qubits=8) 1729 1730 spectrum_from_new_sparse_method = sparse_eigenspectrum( 1731 sparse_ham_number_preserving) 1732 1733 spectrum_from_old_sparse_method = sparse_eigenspectrum( 1734 sparse_ham_restricted_number_preserving) 1735 1736 spectral_deviation = numpy.amax( 1737 numpy.absolute(spectrum_from_new_sparse_method - 1738 spectrum_from_old_sparse_method)) 1739 self.assertAlmostEqual(spectral_deviation, 0.) 1740 1741 def test_number_sz_restricted_spectra_match_molecule(self): 1742 hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian) 1743 1744 sparse_ham_number_sz_preserving = get_number_preserving_sparse_operator( 1745 hamiltonian_fop, 1746 self.molecule.n_qubits, 1747 self.molecule.n_electrons, 1748 spin_preserving=True) 1749 1750 sparse_ham = get_sparse_operator(hamiltonian_fop, 1751 self.molecule.n_qubits) 1752 1753 sparse_ham_restricted_number_sz_preserving = jw_sz_restrict_operator( 1754 sparse_ham, 1755 0, 1756 n_electrons=self.molecule.n_electrons, 1757 n_qubits=self.molecule.n_qubits) 1758 1759 spectrum_from_new_sparse_method = sparse_eigenspectrum( 1760 sparse_ham_number_sz_preserving) 1761 1762 spectrum_from_old_sparse_method = sparse_eigenspectrum( 1763 sparse_ham_restricted_number_sz_preserving) 1764 1765 spectral_deviation = numpy.amax( 1766 numpy.absolute(spectrum_from_new_sparse_method - 1767 spectrum_from_old_sparse_method)) 1768 self.assertAlmostEqual(spectral_deviation, 0.) 1769 1770 def test_number_sz_restricted_spectra_match_hubbard(self): 1771 hamiltonian_fop = self.hubbard_hamiltonian 1772 1773 sparse_ham_number_sz_preserving = get_number_preserving_sparse_operator( 1774 hamiltonian_fop, 8, 4, spin_preserving=True) 1775 1776 sparse_ham = get_sparse_operator(hamiltonian_fop, 8) 1777 1778 sparse_ham_restricted_number_sz_preserving = jw_sz_restrict_operator( 1779 sparse_ham, 0, n_electrons=4, n_qubits=8) 1780 1781 spectrum_from_new_sparse_method = sparse_eigenspectrum( 1782 sparse_ham_number_sz_preserving) 1783 1784 spectrum_from_old_sparse_method = sparse_eigenspectrum( 1785 sparse_ham_restricted_number_sz_preserving) 1786 1787 spectral_deviation = numpy.amax( 1788 numpy.absolute(spectrum_from_new_sparse_method - 1789 spectrum_from_old_sparse_method)) 1790 self.assertAlmostEqual(spectral_deviation, 0.) 1791 1792 1793class GetSparseOperatorQubitTest(unittest.TestCase): 1794 1795 def test_sparse_matrix_Y(self): 1796 term = QubitOperator(((0, 'Y'),)) 1797 sparse_operator = get_sparse_operator(term) 1798 self.assertEqual(list(sparse_operator.data), [1j, -1j]) 1799 self.assertEqual(list(sparse_operator.indices), [1, 0]) 1800 self.assertTrue(is_hermitian(sparse_operator)) 1801 1802 def test_sparse_matrix_ZX(self): 1803 coefficient = 2. 1804 operators = ((0, 'Z'), (1, 'X')) 1805 term = QubitOperator(operators, coefficient) 1806 sparse_operator = get_sparse_operator(term) 1807 self.assertEqual(list(sparse_operator.data), [2., 2., -2., -2.]) 1808 self.assertEqual(list(sparse_operator.indices), [1, 0, 3, 2]) 1809 self.assertTrue(is_hermitian(sparse_operator)) 1810 1811 def test_sparse_matrix_ZIZ(self): 1812 operators = ((0, 'Z'), (2, 'Z')) 1813 term = QubitOperator(operators) 1814 sparse_operator = get_sparse_operator(term) 1815 self.assertEqual(list(sparse_operator.data), 1816 [1, -1, 1, -1, -1, 1, -1, 1]) 1817 self.assertEqual(list(sparse_operator.indices), list(range(8))) 1818 self.assertTrue(is_hermitian(sparse_operator)) 1819 1820 def test_sparse_matrix_combo(self): 1821 qop = (QubitOperator(((0, 'Y'), (1, 'X')), -0.1j) + QubitOperator( 1822 ((0, 'X'), (1, 'Z')), 3. + 2.j)) 1823 sparse_operator = get_sparse_operator(qop) 1824 1825 self.assertEqual( 1826 list(sparse_operator.data), 1827 [3 + 2j, 0.1, 0.1, -3 - 2j, 3 + 2j, -0.1, -0.1, -3 - 2j]) 1828 self.assertEqual(list(sparse_operator.indices), 1829 [2, 3, 2, 3, 0, 1, 0, 1]) 1830 1831 def test_sparse_matrix_zero_1qubit(self): 1832 sparse_operator = get_sparse_operator(QubitOperator((), 0.0), 1) 1833 sparse_operator.eliminate_zeros() 1834 self.assertEqual(len(list(sparse_operator.data)), 0) 1835 self.assertEqual(sparse_operator.shape, (2, 2)) 1836 1837 def test_sparse_matrix_zero_5qubit(self): 1838 sparse_operator = get_sparse_operator(QubitOperator((), 0.0), 5) 1839 sparse_operator.eliminate_zeros() 1840 self.assertEqual(len(list(sparse_operator.data)), 0) 1841 self.assertEqual(sparse_operator.shape, (32, 32)) 1842 1843 def test_sparse_matrix_identity_1qubit(self): 1844 sparse_operator = get_sparse_operator(QubitOperator(()), 1) 1845 self.assertEqual(list(sparse_operator.data), [1] * 2) 1846 self.assertEqual(sparse_operator.shape, (2, 2)) 1847 1848 def test_sparse_matrix_identity_5qubit(self): 1849 sparse_operator = get_sparse_operator(QubitOperator(()), 5) 1850 self.assertEqual(list(sparse_operator.data), [1] * 32) 1851 self.assertEqual(sparse_operator.shape, (32, 32)) 1852 1853 def test_sparse_matrix_linearity(self): 1854 identity = QubitOperator(()) 1855 zzzz = QubitOperator(tuple((i, 'Z') for i in range(4)), 1.0) 1856 1857 sparse1 = get_sparse_operator(identity + zzzz) 1858 sparse2 = get_sparse_operator(identity, 4) + get_sparse_operator(zzzz) 1859 1860 self.assertEqual(list(sparse1.data), [2] * 8) 1861 self.assertEqual(list(sparse1.indices), [0, 3, 5, 6, 9, 10, 12, 15]) 1862 self.assertEqual(list(sparse2.data), [2] * 8) 1863 self.assertEqual(list(sparse2.indices), [0, 3, 5, 6, 9, 10, 12, 15]) 1864 1865 1866class GetSparseOperatorFermionTest(unittest.TestCase): 1867 1868 def test_sparse_matrix_zero_n_qubit(self): 1869 sparse_operator = get_sparse_operator(FermionOperator.zero(), 4) 1870 sparse_operator.eliminate_zeros() 1871 self.assertEqual(len(list(sparse_operator.data)), 0) 1872 self.assertEqual(sparse_operator.shape, (16, 16)) 1873 1874 1875class GetSparseOperatorBosonTest(unittest.TestCase): 1876 1877 def setUp(self): 1878 self.hbar = 1. 1879 self.d = 4 1880 self.b = numpy.diag(numpy.sqrt(numpy.arange(1, self.d)), 1) 1881 self.bd = self.b.conj().T 1882 self.q = numpy.sqrt(self.hbar / 2) * (self.b + self.bd) 1883 1884 def test_sparse_matrix_ladder(self): 1885 sparse_operator = get_sparse_operator(BosonOperator('0'), trunc=self.d) 1886 self.assertTrue(numpy.allclose(sparse_operator.toarray(), self.b)) 1887 self.assertEqual(sparse_operator.shape, (self.d, self.d)) 1888 1889 def test_sparse_matrix_quad(self): 1890 sparse_operator = get_sparse_operator(QuadOperator('q0'), trunc=self.d) 1891 self.assertTrue(numpy.allclose(sparse_operator.toarray(), self.q)) 1892 self.assertEqual(sparse_operator.shape, (self.d, self.d)) 1893 1894 def test_sparse_matrix_error(self): 1895 with self.assertRaises(TypeError): 1896 _ = get_sparse_operator(1) 1897 1898 1899class GetSparseOperatorDiagonalCoulombHamiltonianTest(unittest.TestCase): 1900 1901 def test_diagonal_coulomb_hamiltonian(self): 1902 n_qubits = 5 1903 one_body = random_hermitian_matrix(n_qubits, real=False) 1904 two_body = random_hermitian_matrix(n_qubits, real=True) 1905 constant = numpy.random.randn() 1906 op = DiagonalCoulombHamiltonian(one_body, two_body, constant) 1907 1908 op1 = get_sparse_operator(op) 1909 op2 = get_sparse_operator(jordan_wigner(get_fermion_operator(op))) 1910 diff = op1 - op2 1911 discrepancy = 0. 1912 if diff.nnz: 1913 discrepancy = max(abs(diff.data)) 1914 self.assertAlmostEqual(discrepancy, 0.) 1915