From aefe218fbb0779ccb88b6aa22acad14c570d3418 Mon Sep 17 00:00:00 2001 From: aoowweenn Date: Mon, 10 Oct 2022 11:39:01 +0800 Subject: [PATCH] Fix convention of PolynomialTensor basis change In general, we rotate a matrix M to M' is to apply a rotation matrix R: M' = R @ M @ R.T The corresponding Einstein notation is M'^{p_1p_2} = R^{p_1}_{a_1} R^{p_2}_{b_1} M^{a_1a_2} --- .../ops/representations/polynomial_tensor.py | 7 ++- .../representations/polynomial_tensor_test.py | 45 ++++++++++++++++++- 2 files changed, 47 insertions(+), 5 deletions(-) diff --git a/src/openfermion/ops/representations/polynomial_tensor.py b/src/openfermion/ops/representations/polynomial_tensor.py index 09110341c..ebed64d1e 100644 --- a/src/openfermion/ops/representations/polynomial_tensor.py +++ b/src/openfermion/ops/representations/polynomial_tensor.py @@ -31,8 +31,7 @@ def general_basis_change(general_tensor, rotation_matrix, key): r"""Change the basis of an general interaction tensor. M'^{p_1p_2...p_n} = R^{p_1}_{a_1} R^{p_2}_{a_2} ... - R^{p_n}_{a_n} M^{a_1a_2...a_n} R^{p_n}_{a_n}^T ... - R^{p_2}_{a_2}^T R_{p_1}_{a_1}^T + R^{p_n}_{a_n} M^{a_1a_2...a_n} where R is the rotation matrix, M is the general tensor, M' is the transformed general tensor, and a_k and p_k are indices. The formula uses @@ -67,7 +66,7 @@ def general_basis_change(general_tensor, rotation_matrix, key): # Do the basis change through a single call of numpy.einsum. For example, # for the (1, 1, 0, 0) tensor, the call is: - # numpy.einsum('abcd,aA,bB,cC,dD', + # numpy.einsum('abcd,Aa,Bb,Cc,Dd', # general_tensor, # rotation_matrix.conj(), # rotation_matrix.conj(), @@ -79,7 +78,7 @@ def general_basis_change(general_tensor, rotation_matrix, key): # The 'Aa,Bb,Cc,Dd' part of the subscripts subscripts_rest = ','.join( - chr(ord('a') + i) + chr(ord('A') + i) for i in range(order)) + chr(ord('A') + i) + chr(ord('a') + i) for i in range(order)) subscripts = subscripts_first + ',' + subscripts_rest diff --git a/src/openfermion/ops/representations/polynomial_tensor_test.py b/src/openfermion/ops/representations/polynomial_tensor_test.py index c608309fd..3852e72fb 100644 --- a/src/openfermion/ops/representations/polynomial_tensor_test.py +++ b/src/openfermion/ops/representations/polynomial_tensor_test.py @@ -526,6 +526,47 @@ def test_rotate_basis_reverse(self): polynomial_tensor.rotate_basis(rotation_matrix_reverse) self.assertEqual(polynomial_tensor, want_polynomial_tensor) + def test_rotate_basis_90deg(self): + rotation_matrix_90deg = numpy.zeros((self.n_qubits, self.n_qubits)) + rotation_matrix_90deg[0, 1] = -1 + rotation_matrix_90deg[1, 0] = 1 + + one_body = numpy.zeros((self.n_qubits, self.n_qubits)) + two_body = numpy.zeros( + (self.n_qubits, self.n_qubits, self.n_qubits, self.n_qubits)) + one_body_90deg = numpy.zeros((self.n_qubits, self.n_qubits)) + two_body_90deg = numpy.zeros( + (self.n_qubits, self.n_qubits, self.n_qubits, self.n_qubits)) + i = 0 + j = 0 + i_90deg = pow(self.n_qubits, 2) - 1 + j_90deg = pow(self.n_qubits, 4) - 1 + for p in range(self.n_qubits): + for q in range(self.n_qubits): + one_body[p, q] = i + i = i + 1 + one_body_90deg[p, q] = (-1)**([p, q].count(0))*i_90deg + i_90deg = i_90deg - 1 + for r in range(self.n_qubits): + for s in range(self.n_qubits): + two_body[p, q, r, s] = j + j = j + 1 + two_body_90deg[p, q, r, s] = \ + (-1)**([p, q, r, s].count(0))*j_90deg + j_90deg = j_90deg - 1 + polynomial_tensor = PolynomialTensor({ + (): self.constant, + (1, 0): one_body, + (1, 1, 0, 0): two_body + }) + want_polynomial_tensor = PolynomialTensor({ + (): self.constant, + (1, 0): one_body_90deg, + (1, 1, 0, 0): two_body_90deg + }) + polynomial_tensor.rotate_basis(rotation_matrix_90deg) + self.assertEqual(polynomial_tensor, want_polynomial_tensor) + def test_rotate_basis_quadratic_hamiltonian_real(self): self.do_rotate_basis_quadratic_hamiltonian(True) @@ -545,7 +586,7 @@ def do_rotate_basis_quadratic_hamiltonian(self, real): # Rotate a basis where the Hamiltonian is diagonal _, diagonalizing_unitary, _ = ( quad_ham.diagonalizing_bogoliubov_transform()) - quad_ham.rotate_basis(diagonalizing_unitary.T) + quad_ham.rotate_basis(diagonalizing_unitary) # Check that the rotated Hamiltonian is diagonal with the correct # orbital energies @@ -567,6 +608,8 @@ def test_rotate_basis_max_order(self): self.assertEqual(tensor, want_tensor) # I originally wanted to test 25 and 26, but it turns out that # numpy.einsum complains "too many subscripts in einsum" before 26. + # Currently, the sum of the number of output labels and combined labels + # can't not exceed NPY_MAXDIMS(=32). for order in [27, 28]: with self.assertRaises(ValueError):