diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py
index 15cc467..fdd42be 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian.py
@@ -17,7 +17,7 @@ from __future__ import absolute_import
 import numpy
 
 from openfermion.config import EQ_TOLERANCE
-from openfermion.ops import PolynomialTensor
+from openfermion.ops import FermionOperator, PolynomialTensor
 
 
 class QuadraticHamiltonianError(Exception):
@@ -67,15 +67,17 @@ class QuadraticHamiltonian(PolynomialTensor):
                     hermitian_part -
                     chemical_potential * numpy.eye(n_qubits))
 
-        # Initialize antisymmetric part
         if antisymmetric_part is None:
-            antisymmetric_part = numpy.zeros((n_qubits, n_qubits), complex)
+            super(QuadraticHamiltonian, self).__init__(
+                    {(): constant,
+                     (1, 0): combined_hermitian_part})
+        else:
+            super(QuadraticHamiltonian, self).__init__(
+                    {(): constant,
+                     (1, 0): combined_hermitian_part,
+                     (1, 1): .5 * antisymmetric_part,
+                     (0, 0): -.5 * antisymmetric_part.conj()})
 
-        super(QuadraticHamiltonian, self).__init__(
-                {(): constant,
-                 (1, 0): combined_hermitian_part,
-                 (1, 1): .5 * antisymmetric_part,
-                 (0, 0): -.5 * antisymmetric_part.conj()})
         self.chemical_potential = chemical_potential
 
     def combined_hermitian_part(self):
@@ -92,9 +94,100 @@ class QuadraticHamiltonian(PolynomialTensor):
 
     def antisymmetric_part(self):
         """Return the antisymmetric part."""
-        return 2. * self.n_body_tensors[1, 1].copy()
+        if (1, 1) in self.n_body_tensors:
+            return 2. * self.n_body_tensors[1, 1].copy()
+        else:
+            return numpy.zeros((self.n_qubits, self.n_qubits), complex)
 
     def conserves_particle_number(self):
         """Return whether this Hamiltonian conserves particle number."""
         discrepancy = numpy.max(numpy.abs(self.antisymmetric_part()))
         return discrepancy < EQ_TOLERANCE
+
+    def majorana_form(self):
+        """Return the Majorana represention of the Hamiltonian.
+
+        Any quadratic Hamiltonian can be written in the form
+
+            constant + i / 2 \sum_{j, k} A_{jk} s_j s_k.
+
+        where the s_i are normalized Majorana fermion operators:
+
+            s_j = i / sqrt(2) (a^\dagger_j - a_j)
+            s_{j + n_qubits} = 1 / sqrt(2) (a^\dagger_j + a_j)
+
+        and A is a (2 * n_qubits) x (2 * n_qubits) real antisymmetric matrix.
+        This function returns the matrix A and the constant.
+        """
+        hermitian_part = self.combined_hermitian_part()
+        antisymmetric_part = self.antisymmetric_part()
+
+        # Compute the Majorana matrix using block matrix manipulations
+        majorana_matrix = numpy.zeros((2 * self.n_qubits, 2 * self.n_qubits))
+        # Set upper left block
+        majorana_matrix[:self.n_qubits, :self.n_qubits] = numpy.real(-.5j * (
+                hermitian_part - hermitian_part.conj() +
+                antisymmetric_part - antisymmetric_part.conj()))
+        # Set upper right block
+        majorana_matrix[:self.n_qubits, self.n_qubits:] = numpy.real(.5 * (
+                hermitian_part + hermitian_part.conj() -
+                antisymmetric_part - antisymmetric_part.conj()))
+        # Set lower left block
+        majorana_matrix[self.n_qubits:, :self.n_qubits] = numpy.real(-.5 * (
+                hermitian_part + hermitian_part.conj() +
+                antisymmetric_part + antisymmetric_part.conj()))
+        # Set lower right block
+        majorana_matrix[self.n_qubits:, self.n_qubits:] = numpy.real(-.5j * (
+                hermitian_part - hermitian_part.conj() -
+                antisymmetric_part + antisymmetric_part.conj()))
+
+        # Compute the constant
+        majorana_constant = (.5 * numpy.real(numpy.trace(hermitian_part)) +
+                             self.n_body_tensors[()])
+
+        return majorana_matrix, majorana_constant
+
+
+def majorana_operator(term=None, coefficient=1.):
+    """Initialize a Majorana operator.
+
+    Args:
+        term(tuple): The first element of the tuple indicates the mode
+            on which the Majorana operator acts, starting from zero.
+            The second element of the tuple is an integer, either 1 or 0,
+            indicating which type of Majorana operator it is:
+                type 1: 1 / sqrt(2) (a^\dagger_j + a_j)
+                type 0: i / sqrt(2) (a^\dagger_j - a_j)
+            where the a^\dagger_j and a_j are the usual fermionic ladder
+            operators.
+            Default will result in the zero operator.
+        coefficient(complex or float, optional): The coefficient of the term.
+            Default value is 1.0.
+
+    Returns:
+        FermionOperator
+    """
+    if not isinstance(coefficient, (int, float, complex)):
+        raise ValueError('Coefficient must be scalar.')
+
+    if term is None:
+        # Return zero operator
+        return FermionOperator()
+    elif isinstance(term, tuple):
+        mode, operator_type = term
+        if operator_type == 1:
+            majorana_op = FermionOperator(
+                    ((mode, 1),), coefficient / numpy.sqrt(2.))
+            majorana_op += FermionOperator(
+                    ((mode, 0),), coefficient / numpy.sqrt(2.))
+        elif operator_type == 0:
+            majorana_op = FermionOperator(
+                    ((mode, 1),), 1.j * coefficient / numpy.sqrt(2.))
+            majorana_op -= FermionOperator(
+                    ((mode, 0),), 1.j * coefficient / numpy.sqrt(2.))
+        else:
+            raise ValueError('Operator specified incorrectly.')
+        return majorana_op
+    # Invalid input.
+    else:
+        raise ValueError('Operator specified incorrectly.')
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index 9d5bb20..d75656f 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -20,8 +20,9 @@ from ._operator_utils import (commutator, count_qubits,
                               get_file_path, inverse_fourier_transform,
                               is_identity, load_operator, save_operator)
 
-from ._slater_determinants import (givens_decomposition,
-                                   fermionic_gaussian_decomposition)
+from ._slater_determinants import (fermionic_gaussian_decomposition,
+                                   givens_decomposition,
+                                   ground_state_preparation_circuit)
 
 from ._sparse_tools import (expectation,
                             expectation_computational_basis_state,
diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py
index 5ebecd6..b3ae971 100644
--- a/src/openfermion/utils/_slater_determinants.py
+++ b/src/openfermion/utils/_slater_determinants.py
@@ -18,6 +18,63 @@ import numpy
 from scipy.linalg import schur
 
 from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import QuadraticHamiltonian
+
+
+def ground_state_preparation_circuit(quadratic_hamiltonian):
+    """Obtain a description of a circuit which prepares the ground state
+    of a quadratic Hamiltonian.
+
+    Args:
+        quadratic_hamiltonian(QuadraticHamiltonian):
+            The Hamiltonian whose ground state is desired.
+
+    Returns:
+        circuit_description(list[tuple]):
+            A list of operations describing the circuit. Each operation
+            is a tuple of objects describing elementary operations that
+            can be performed in parallel. Each elementary operation
+            is either the string 'pht', indicating a particle-hole
+            transformation on the last fermionic mode, or a tuple of
+            the form (i, j, theta, phi), indicating a Givens rotation
+            of qubits i and j by angles theta and phi.
+    """
+    if not isinstance(quadratic_hamiltonian, QuadraticHamiltonian):
+        raise ValueError('Input must be an instance of QuadraticHamiltonian.')
+
+    if quadratic_hamiltonian.conserves_particle_number():
+        # The Hamiltonian conserves particle number, so we don't need
+        # to use the most general procedure.
+        hermitian_matrix = quadratic_hamiltonian.combined_hermitian_part()
+        # Get the unitary rows which represent the ground state
+        # Slater determinant
+        energies, diagonalizing_unitary = numpy.linalg.eigh(
+                hermitian_matrix)
+        num_negative_energies = numpy.count_nonzero(
+                energies < -EQ_TOLERANCE)
+        slater_determinant_matrix = diagonalizing_unitary.T[
+                :num_negative_energies]
+        # Get the circuit description
+        left_unitary, decomposition, diagonal = givens_decomposition(
+                slater_determinant_matrix)
+        circuit_description = decomposition
+    else:
+        # The Hamiltonian does not conserve particle number, so we
+        # need to use the most general procedure.
+        majorana_matrix, majorana_constant = (
+                quadratic_hamiltonian.majorana_form())
+        # Get the unitary rows which represent the ground state
+        # Slater determinant
+        diagonalizing_unitary = diagonalizing_fermionic_unitary(
+                majorana_matrix)
+        slater_determinant_matrix = diagonalizing_unitary[
+                quadratic_hamiltonian.n_qubits:]
+        # Get the circuit description
+        left_unitary, decomposition, diagonal = (
+                fermionic_gaussian_decomposition(
+                    slater_determinant_matrix))
+        circuit_description = decomposition
+    return circuit_description
 
 
 def fermionic_gaussian_decomposition(unitary_rows):
@@ -49,7 +106,7 @@ def fermionic_gaussian_decomposition(unitary_rows):
     something like [('pht', ), (G_1, ), ('pht', G_2), ... ].
     The objects within a tuple are either the string 'pht', which indicates
     a particle-hole transformation on the last fermionic mode, or a tuple
-    of the form (i, j, theta, phi), which indicates a Givens roation
+    of the form (i, j, theta, phi), which indicates a Givens rotation
     of rows i and j by angles theta and phi.
 
     Args:
@@ -148,109 +205,6 @@ def fermionic_gaussian_decomposition(unitary_rows):
     return left_unitary, decomposition, diagonal
 
 
-def diagonalizing_fermionic_unitary(antisymmetric_matrix):
-    """Compute the unitary that diagonalizes a quadratic Hamiltonian.
-
-    The input matrix represents a quadratic Hamiltonian in the Majorana basis.
-    The output matrix is a unitary that represents a transformation (mixing)
-    of the fermionic ladder operators. We use the convention that the
-    creation operators are listed before the annihilation operators.
-    The returned unitary has additional structure which ensures
-    that the transformed ladder operators also satisfy the fermionic
-    anticommutation relations.
-
-    Args:
-        antisymmetric_matrix(ndarray): A (2 * n_qubits) x (2 * n_qubits)
-            antisymmetric matrix representing a quadratic Hamiltonian in the
-            Majorana basis.
-    Returns:
-        diagonalizing_unitary(ndarray): A (2 * n_qubits) x (2 * n_qubits)
-            unitary matrix representing a transformation of the fermionic
-            ladder operators.
-    """
-    m, n = antisymmetric_matrix.shape
-    n_qubits = n // 2
-
-    # Get the orthogonal transformation that puts antisymmetric_matrix
-    # into canonical form
-    canonical, orthogonal = antisymmetric_canonical_form(antisymmetric_matrix)
-
-    # Create the matrix that converts between fermionic ladder and
-    # Majorana bases
-    normalized_identity = numpy.eye(n_qubits, dtype=complex) / numpy.sqrt(2.)
-    majorana_basis_change = numpy.eye(
-            2 * n_qubits, dtype=complex) / numpy.sqrt(2.)
-    majorana_basis_change[n_qubits:, n_qubits:] *= -1.j
-    majorana_basis_change[:n_qubits, n_qubits:] = normalized_identity
-    majorana_basis_change[n_qubits:, :n_qubits] = 1.j * normalized_identity
-
-    # Compute the unitary and return
-    diagonalizing_unitary = majorana_basis_change.T.conj().dot(
-            orthogonal.dot(majorana_basis_change))
-
-    return diagonalizing_unitary
-
-
-def antisymmetric_canonical_form(antisymmetric_matrix):
-    """Compute the canonical form of an antisymmetric matrix.
-
-    The input is a real, antisymmetric n x n matrix A, where n is even.
-    Its canonical form is::
-
-        A = R^T C R
-
-    where R is a real, orthogonal matrix and C has the form::
-
-        [  0     D ]
-        [ -D     0 ]
-
-    where D is a diagonal matrix with nonnegative entries.
-
-    Args:
-        antisymmetric_matrix(ndarray): An antisymmetric matrix with even
-            dimension.
-
-    Returns:
-        canonical(ndarray): The canonical form C of antisymmetric_matrix
-        orthogonal(ndarray): The orthogonal transformation R.
-    """
-    m, n = antisymmetric_matrix.shape
-
-    if m != n or n % 2 != 0:
-        raise ValueError('The input matrix must be square with even '
-                         'dimension.')
-
-    # Check that input matrix is antisymmetric
-    matrix_plus_transpose = antisymmetric_matrix + antisymmetric_matrix.T
-    maxval = numpy.max(numpy.abs(matrix_plus_transpose))
-    if maxval > EQ_TOLERANCE:
-        raise ValueError('The input matrix must be antisymmetric.')
-
-    # Compute Schur decomposition
-    canonical, orthogonal = schur(antisymmetric_matrix, output='real')
-
-    # The returned form is block diagonal; we need to permute rows and columns
-    # to put it into the form we want
-    num_blocks = n // 2
-    for i in range(1, num_blocks, 2):
-        swap_rows(canonical, i, num_blocks + i - 1)
-        swap_columns(canonical, i, num_blocks + i - 1)
-        swap_columns(orthogonal, i, num_blocks + i - 1)
-        if num_blocks % 2 != 0:
-            swap_rows(canonical, num_blocks - 1, num_blocks + i)
-            swap_columns(canonical, num_blocks - 1, num_blocks + i)
-            swap_columns(orthogonal, num_blocks - 1, num_blocks + i)
-
-    # Now we permute so that the upper right block is non-negative
-    for i in range(num_blocks):
-        if canonical[i, num_blocks + i] < -EQ_TOLERANCE:
-            swap_rows(canonical, i, num_blocks + i)
-            swap_columns(canonical, i, num_blocks + i)
-            swap_columns(orthogonal, i, num_blocks + i)
-
-    return canonical, orthogonal.T
-
-
 def givens_decomposition(unitary_rows):
     """Decompose a matrix into a sequence of Givens rotations.
 
@@ -379,6 +333,109 @@ def givens_decomposition(unitary_rows):
     return left_unitary, givens_rotations, diagonal
 
 
+def diagonalizing_fermionic_unitary(antisymmetric_matrix):
+    """Compute the unitary that diagonalizes a quadratic Hamiltonian.
+
+    The input matrix represents a quadratic Hamiltonian in the Majorana basis.
+    The output matrix is a unitary that represents a transformation (mixing)
+    of the fermionic ladder operators. We use the convention that the
+    creation operators are listed before the annihilation operators.
+    The returned unitary has additional structure which ensures
+    that the transformed ladder operators also satisfy the fermionic
+    anticommutation relations.
+
+    Args:
+        antisymmetric_matrix(ndarray): A (2 * n_qubits) x (2 * n_qubits)
+            antisymmetric matrix representing a quadratic Hamiltonian in the
+            Majorana basis.
+    Returns:
+        diagonalizing_unitary(ndarray): A (2 * n_qubits) x (2 * n_qubits)
+            unitary matrix representing a transformation of the fermionic
+            ladder operators.
+    """
+    m, n = antisymmetric_matrix.shape
+    n_qubits = n // 2
+
+    # Get the orthogonal transformation that puts antisymmetric_matrix
+    # into canonical form
+    canonical, orthogonal = antisymmetric_canonical_form(antisymmetric_matrix)
+
+    # Create the matrix that converts between fermionic ladder and
+    # Majorana bases
+    normalized_identity = numpy.eye(n_qubits, dtype=complex) / numpy.sqrt(2.)
+    majorana_basis_change = numpy.eye(
+            2 * n_qubits, dtype=complex) / numpy.sqrt(2.)
+    majorana_basis_change[n_qubits:, n_qubits:] *= -1.j
+    majorana_basis_change[:n_qubits, n_qubits:] = normalized_identity
+    majorana_basis_change[n_qubits:, :n_qubits] = 1.j * normalized_identity
+
+    # Compute the unitary and return
+    diagonalizing_unitary = majorana_basis_change.T.conj().dot(
+            orthogonal.dot(majorana_basis_change))
+
+    return diagonalizing_unitary
+
+
+def antisymmetric_canonical_form(antisymmetric_matrix):
+    """Compute the canonical form of an antisymmetric matrix.
+
+    The input is a real, antisymmetric n x n matrix A, where n is even.
+    Its canonical form is::
+
+        A = R^T C R
+
+    where R is a real, orthogonal matrix and C has the form::
+
+        [  0     D ]
+        [ -D     0 ]
+
+    where D is a diagonal matrix with nonnegative entries.
+
+    Args:
+        antisymmetric_matrix(ndarray): An antisymmetric matrix with even
+            dimension.
+
+    Returns:
+        canonical(ndarray): The canonical form C of antisymmetric_matrix
+        orthogonal(ndarray): The orthogonal transformation R.
+    """
+    m, n = antisymmetric_matrix.shape
+
+    if m != n or n % 2 != 0:
+        raise ValueError('The input matrix must be square with even '
+                         'dimension.')
+
+    # Check that input matrix is antisymmetric
+    matrix_plus_transpose = antisymmetric_matrix + antisymmetric_matrix.T
+    maxval = numpy.max(numpy.abs(matrix_plus_transpose))
+    if maxval > EQ_TOLERANCE:
+        raise ValueError('The input matrix must be antisymmetric.')
+
+    # Compute Schur decomposition
+    canonical, orthogonal = schur(antisymmetric_matrix, output='real')
+
+    # The returned form is block diagonal; we need to permute rows and columns
+    # to put it into the form we want
+    num_blocks = n // 2
+    for i in range(1, num_blocks, 2):
+        swap_rows(canonical, i, num_blocks + i - 1)
+        swap_columns(canonical, i, num_blocks + i - 1)
+        swap_columns(orthogonal, i, num_blocks + i - 1)
+        if num_blocks % 2 != 0:
+            swap_rows(canonical, num_blocks - 1, num_blocks + i)
+            swap_columns(canonical, num_blocks - 1, num_blocks + i)
+            swap_columns(orthogonal, num_blocks - 1, num_blocks + i)
+
+    # Now we permute so that the upper right block is non-negative
+    for i in range(num_blocks):
+        if canonical[i, num_blocks + i] < -EQ_TOLERANCE:
+            swap_rows(canonical, i, num_blocks + i)
+            swap_columns(canonical, i, num_blocks + i)
+            swap_columns(orthogonal, i, num_blocks + i)
+
+    return canonical, orthogonal.T
+
+
 def givens_matrix_elements(a, b):
     """Compute the matrix elements of the Givens rotation that zeroes out one
     of two row entries.

Circular dependency when importing from utils in ops
I'm working on the file `ops/_quadratic_hamiltonian.py`. In there, I add the line

```from openfermion.utils._slater_determinants import swap_columns, swap_rows```

This causes the tests in the `hamiltonians` folder to fail; i.e., when I run `pytest hamiltonians/`, I get

```
________________ ERROR collecting src/openfermion/hamiltonians/_chemical_series_test.py ________________
ImportError while importing test module '/home/kjs/Projects/OpenFermion/src/openfermion/hamiltonians/_chemical_series_test.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
hamiltonians/__init__.py:13: in <module>
    from ._chemical_series import (make_atomic_ring,
hamiltonians/_chemical_series.py:18: in <module>
    from openfermion.hamiltonians._molecular_data import (MolecularData,
hamiltonians/_molecular_data.py:21: in <module>
    from openfermion.ops import InteractionOperator, InteractionRDM
ops/__init__.py:23: in <module>
    from ._quadratic_hamiltonian import QuadraticHamiltonian
ops/_quadratic_hamiltonian.py:22: in <module>
    from openfermion.utils._slater_determinants import swap_columns, swap_rows
utils/__init__.py:39: in <module>
    from ._trotter_error import error_bound, error_operator
utils/_trotter_error.py:20: in <module>
    from openfermion.hamiltonians import MolecularData
E   ImportError: cannot import name 'MolecularData'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Interrupted: 1 errors during collection !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```

It seems that there is a circular dependency
MolecularData -> ops -> utils -> MolecularData
 and any file in the `ops` folder which imports anything from `utils` will cause this error. What is the best way to deal with it?
diff --git a/src/openfermion/ops/_quadratic_hamiltonian_test.py b/src/openfermion/ops/_quadratic_hamiltonian_test.py
index e852031..321f925 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian_test.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian_test.py
@@ -16,20 +16,29 @@ from __future__ import absolute_import
 import numpy
 import unittest
 
-from openfermion.ops import QuadraticHamiltonian
+from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import (FermionOperator,
+                             QuadraticHamiltonian,
+                             normal_ordered)
+from openfermion.ops._quadratic_hamiltonian import majorana_operator
+from openfermion.transforms import get_fermion_operator
 
 
 class QuadraticHamiltoniansTest(unittest.TestCase):
 
     def setUp(self):
         self.n_qubits = 5
-        self.constant = 1.
+        self.constant = 1.7
         self.chemical_potential = 2.
 
         # Obtain random Hermitian and antisymmetric matrices
-        rand_mat = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
         self.hermitian_mat = rand_mat + rand_mat.T.conj()
-        rand_mat = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
         self.antisymmetric_mat = rand_mat - rand_mat.T
 
         self.combined_hermitian = (
@@ -82,3 +91,43 @@ class QuadraticHamiltoniansTest(unittest.TestCase):
         """Test checking whether Hamiltonian conserves particle number."""
         self.assertTrue(self.quad_ham_pc.conserves_particle_number())
         self.assertFalse(self.quad_ham_npc.conserves_particle_number())
+
+    def test_majorana_form(self):
+        """Test getting the Majorana form."""
+        majorana_matrix, majorana_constant = self.quad_ham_npc.majorana_form()
+        # Convert the Majorana form to a FermionOperator
+        majorana_op = FermionOperator((), majorana_constant)
+        for i in range(2 * self.n_qubits):
+            if i < self.n_qubits:
+                left_op = majorana_operator((i, 1))
+            else:
+                left_op = majorana_operator((i - self.n_qubits, 0))
+            for j in range(2 * self.n_qubits):
+                if j < self.n_qubits:
+                    right_op = majorana_operator((j, 1), majorana_matrix[i, j])
+                else:
+                    right_op = majorana_operator((j - self.n_qubits, 0),
+                                                 majorana_matrix[i, j])
+                majorana_op += .5j * left_op * right_op
+        # Get FermionOperator for original Hamiltonian
+        fermion_operator = normal_ordered(
+                get_fermion_operator(self.quad_ham_npc))
+        self.assertTrue(
+                normal_ordered(majorana_op).isclose(fermion_operator))
+
+
+class MajoranaOperatorTest(unittest.TestCase):
+
+    def test_none_term(self):
+        majorana_op = majorana_operator()
+        self.assertTrue(majorana_operator().isclose(FermionOperator()))
+
+    def test_bad_coefficient(self):
+        with self.assertRaises(ValueError):
+            majorana_op = majorana_operator((1, 1), 'a')
+
+    def test_bad_term(self):
+        with self.assertRaises(ValueError):
+            majorana_op = majorana_operator((2, 2))
+        with self.assertRaises(ValueError):
+            majorana_op = majorana_operator('a')
diff --git a/src/openfermion/utils/_slater_determinants_test.py b/src/openfermion/utils/_slater_determinants_test.py
index d0c838c..da94ad1 100644
--- a/src/openfermion/utils/_slater_determinants_test.py
+++ b/src/openfermion/utils/_slater_determinants_test.py
@@ -18,8 +18,10 @@ import unittest
 from scipy.linalg import qr
 
 from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import QuadraticHamiltonian
 from openfermion.utils import (fermionic_gaussian_decomposition,
-                               givens_decomposition)
+                               givens_decomposition,
+                               ground_state_preparation_circuit)
 from openfermion.utils._slater_determinants import (
         antisymmetric_canonical_form,
         diagonalizing_fermionic_unitary,
@@ -28,6 +30,46 @@ from openfermion.utils._slater_determinants import (
         swap_rows)
 
 
+class GroundStatePreparationCircuitTest(unittest.TestCase):
+
+    def setUp(self):
+        self.n_qubits = 5
+        self.constant = 1.7
+        self.chemical_potential = 2.
+
+        # Obtain random Hermitian and antisymmetric matrices
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        self.hermitian_mat = rand_mat + rand_mat.T.conj()
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        self.antisymmetric_mat = rand_mat - rand_mat.T
+
+        self.combined_hermitian = (
+                self.hermitian_mat -
+                self.chemical_potential * numpy.eye(self.n_qubits))
+
+        # Initialize a particle-number-conserving Hamiltonian
+        self.quad_ham_pc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat)
+
+        # Initialize a non-particle-number-conserving Hamiltonian
+        self.quad_ham_npc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat, self.antisymmetric_mat,
+                self.chemical_potential)
+
+    def test_no_error(self):
+        """Test that the procedure runs without error."""
+        # Test a particle-number-conserving Hamiltonian
+        circuit_description = (
+                ground_state_preparation_circuit(self.quad_ham_pc))
+        # Test a non-particle-number-conserving Hamiltonian
+        circuit_description = (
+                ground_state_preparation_circuit(self.quad_ham_npc))
+
+
 class GivensDecompositionTest(unittest.TestCase):
 
     def test_bad_dimensions(self):

pip install -e .[dev]
parse_log_pytest
requirements.txt
pytest pytest-cov coveralls sphinx sphinx_rtd_theme
apt-get update apt-get install -y gcc g++
3.9
requirements.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.16
anyio==4.9.0
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
asttokens==3.0.0
async-lru==2.0.5
attrs==25.3.0
babel==2.17.0
beautifulsoup4==4.13.3
bleach==6.2.0
certifi==2025.1.31
cffi==1.17.1
charset-normalizer==3.4.1
comm==0.2.2
contourpy==1.3.0
coverage==7.8.0
coveralls==4.0.1
cvxopt==1.3.2
cycler==0.12.1
debugpy==1.8.13
decorator==5.2.1
defusedxml==0.7.1
docopt==0.6.2
docutils==0.21.2
exceptiongroup==1.2.2
executing==2.2.0
fastjsonschema==2.21.1
fonttools==4.56.0
fqdn==1.5.1
future==1.0.0
h11==0.14.0
h5py==3.13.0
httpcore==1.0.7
httpx==0.28.1
idna==3.10
imagesize==1.4.1
importlib_metadata==8.6.1
importlib_resources==6.5.2
iniconfig==2.1.0
ipykernel==6.29.5
ipython==8.18.1
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.10.0
jsonpointer==3.0.0
jsonschema==4.23.0
jsonschema-specifications==2024.10.1
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.12.0
jupyter-lsp==2.2.5
jupyter_client==8.6.3
jupyter_core==5.7.2
jupyter_server==2.15.0
jupyter_server_terminals==0.5.3
jupyterlab==4.3.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.27.3
jupyterlab_widgets==3.0.13
kiwisolver==1.4.7
MarkupSafe==3.0.2
matplotlib==3.9.4
matplotlib-inline==0.1.7
mistune==3.1.3
nbclient==0.10.2
nbconvert==7.16.6
nbformat==5.10.4
nest-asyncio==1.6.0
networkx==3.2.1
notebook==7.3.3
notebook_shim==0.2.4
numpy==2.0.2
-e git+https://github.com/quantumlib/OpenFermion.git@e51fe9c381fe29c04a83d63c6d708d4af2bb7a63#egg=openfermion
overrides==7.7.0
packaging==24.2
pandocfilters==1.5.1
parso==0.8.4
pexpect==4.9.0
pillow==11.1.0
platformdirs==4.3.7
pluggy==1.5.0
prometheus_client==0.21.1
prompt_toolkit==3.0.50
psutil==7.0.0
ptyprocess==0.7.0
pure_eval==0.2.3
pycparser==2.22
Pygments==2.19.1
pyparsing==3.2.3
pytest==8.3.5
pytest-cov==6.0.0
python-dateutil==2.9.0.post0
python-json-logger==3.3.0
PyYAML==6.0.2
pyzmq==26.3.0
referencing==0.36.2
requests==2.32.3
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.24.0
scipy==1.13.1
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
snowballstemmer==2.2.0
soupsieve==2.6
Sphinx==7.4.7
sphinx-rtd-theme==3.0.2
sphinxcontrib-applehelp==2.0.0
sphinxcontrib-devhelp==2.0.0
sphinxcontrib-htmlhelp==2.1.0
sphinxcontrib-jquery==4.1
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==2.0.0
sphinxcontrib-serializinghtml==2.0.0
stack-data==0.6.3
terminado==0.18.1
tinycss2==1.4.0
tomli==2.2.1
tornado==6.4.2
traitlets==5.14.3
types-python-dateutil==2.9.0.20241206
typing_extensions==4.13.0
uri-template==1.3.0
urllib3==2.3.0
wcwidth==0.2.13
webcolors==24.11.1
webencodings==0.5.1
websocket-client==1.8.0
widgetsnbextension==4.0.13
zipp==3.21.0

name: OpenFermion
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=3.0.16=h5eee18b_0
  - pip=25.0=py39h06a4308_0
  - python=3.9.21=he870216_1
  - readline=8.2=h5eee18b_0
  - setuptools=75.8.0=py39h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - tzdata=2025a=h04d1e81_0
  - wheel=0.45.1=py39h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.16
      - anyio==4.9.0
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - arrow==1.3.0
      - asttokens==3.0.0
      - async-lru==2.0.5
      - attrs==25.3.0
      - babel==2.17.0
      - beautifulsoup4==4.13.3
      - bleach==6.2.0
      - certifi==2025.1.31
      - cffi==1.17.1
      - charset-normalizer==3.4.1
      - comm==0.2.2
      - contourpy==1.3.0
      - coverage==7.8.0
      - coveralls==4.0.1
      - cvxopt==1.3.2
      - cycler==0.12.1
      - debugpy==1.8.13
      - decorator==5.2.1
      - defusedxml==0.7.1
      - docopt==0.6.2
      - docutils==0.21.2
      - exceptiongroup==1.2.2
      - executing==2.2.0
      - fastjsonschema==2.21.1
      - fonttools==4.56.0
      - fqdn==1.5.1
      - future==1.0.0
      - h11==0.14.0
      - h5py==3.13.0
      - httpcore==1.0.7
      - httpx==0.28.1
      - idna==3.10
      - imagesize==1.4.1
      - importlib-metadata==8.6.1
      - importlib-resources==6.5.2
      - iniconfig==2.1.0
      - ipykernel==6.29.5
      - ipython==8.18.1
      - ipywidgets==8.1.5
      - isoduration==20.11.0
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.10.0
      - jsonpointer==3.0.0
      - jsonschema==4.23.0
      - jsonschema-specifications==2024.10.1
      - jupyter==1.1.1
      - jupyter-client==8.6.3
      - jupyter-console==6.6.3
      - jupyter-core==5.7.2
      - jupyter-events==0.12.0
      - jupyter-lsp==2.2.5
      - jupyter-server==2.15.0
      - jupyter-server-terminals==0.5.3
      - jupyterlab==4.3.6
      - jupyterlab-pygments==0.3.0
      - jupyterlab-server==2.27.3
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.7
      - markupsafe==3.0.2
      - matplotlib==3.9.4
      - matplotlib-inline==0.1.7
      - mistune==3.1.3
      - nbclient==0.10.2
      - nbconvert==7.16.6
      - nbformat==5.10.4
      - nest-asyncio==1.6.0
      - networkx==3.2.1
      - notebook==7.3.3
      - notebook-shim==0.2.4
      - numpy==2.0.2
      - overrides==7.7.0
      - packaging==24.2
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pexpect==4.9.0
      - pillow==11.1.0
      - platformdirs==4.3.7
      - pluggy==1.5.0
      - prometheus-client==0.21.1
      - prompt-toolkit==3.0.50
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pure-eval==0.2.3
      - pycparser==2.22
      - pygments==2.19.1
      - pyparsing==3.2.3
      - pytest==8.3.5
      - pytest-cov==6.0.0
      - python-dateutil==2.9.0.post0
      - python-json-logger==3.3.0
      - pyyaml==6.0.2
      - pyzmq==26.3.0
      - referencing==0.36.2
      - requests==2.32.3
      - rfc3339-validator==0.1.4
      - rfc3986-validator==0.1.1
      - rpds-py==0.24.0
      - scipy==1.13.1
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - snowballstemmer==2.2.0
      - soupsieve==2.6
      - sphinx==7.4.7
      - sphinx-rtd-theme==3.0.2
      - sphinxcontrib-applehelp==2.0.0
      - sphinxcontrib-devhelp==2.0.0
      - sphinxcontrib-htmlhelp==2.1.0
      - sphinxcontrib-jquery==4.1
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==2.0.0
      - sphinxcontrib-serializinghtml==2.0.0
      - stack-data==0.6.3
      - terminado==0.18.1
      - tinycss2==1.4.0
      - tomli==2.2.1
      - tornado==6.4.2
      - traitlets==5.14.3
      - types-python-dateutil==2.9.0.20241206
      - typing-extensions==4.13.0
      - uri-template==1.3.0
      - urllib3==2.3.0
      - wcwidth==0.2.13
      - webcolors==24.11.1
      - webencodings==0.5.1
      - websocket-client==1.8.0
      - widgetsnbextension==4.0.13
      - zipp==3.21.0
prefix: /opt/conda/envs/OpenFermion

src/openfermion/ops/_quadratic_hamiltonian_test.py::QuadraticHamiltoniansTest::test_antisymmetric_part src/openfermion/ops/_quadratic_hamiltonian_test.py::QuadraticHamiltoniansTest::test_combined_hermitian_part src/openfermion/ops/_quadratic_hamiltonian_test.py::QuadraticHamiltoniansTest::test_conserves_particle_number src/openfermion/ops/_quadratic_hamiltonian_test.py::QuadraticHamiltoniansTest::test_hermitian_part src/openfermion/ops/_quadratic_hamiltonian_test.py::QuadraticHamiltoniansTest::test_majorana_form src/openfermion/ops/_quadratic_hamiltonian_test.py::MajoranaOperatorTest::test_bad_coefficient src/openfermion/ops/_quadratic_hamiltonian_test.py::MajoranaOperatorTest::test_bad_term src/openfermion/ops/_quadratic_hamiltonian_test.py::MajoranaOperatorTest::test_none_term src/openfermion/utils/_slater_determinants_test.py::GroundStatePreparationCircuitTest::test_no_error src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_3 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_4 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_5 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_6 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_7 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_8 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_3_by_9 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_4_by_5 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_4_by_9 src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_antidiagonal src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_bad_dimensions src/openfermion/utils/_slater_determinants_test.py::GivensDecompositionTest::test_identity src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_bad_constraints src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_bad_dimensions src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_3 src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_4 src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_5 src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_6 src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_7 src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_8 src/openfermion/utils/_slater_determinants_test.py::FermionicGaussianDecompositionTest::test_n_equals_9 src/openfermion/utils/_slater_determinants_test.py::DiagonalizingFermionicUnitaryTest::test_bad_dimensions src/openfermion/utils/_slater_determinants_test.py::DiagonalizingFermionicUnitaryTest::test_n_equals_3 src/openfermion/utils/_slater_determinants_test.py::DiagonalizingFermionicUnitaryTest::test_not_antisymmetric src/openfermion/utils/_slater_determinants_test.py::AntisymmetricCanonicalFormTest::test_canonical src/openfermion/utils/_slater_determinants_test.py::AntisymmetricCanonicalFormTest::test_equality
src/openfermion/ops/_quadratic_hamiltonian.py src/openfermion/utils/__init__.py src/openfermion/utils/_slater_determinants.py
src/openfermion/ops/_quadratic_hamiltonian.py src/openfermion/utils/__init__.py src/openfermion/utils/_slater_determinants.py