Added code to compute length scales based on Wigner-Seitz radius (#87)

* Added code to compute length scales based on Wigner-Seitz radius

* Added test and docs for exception that is raised for invalid dimensions

src/fermilib/utils/_plane_wave_hamiltonian.py

@@ -26,6 +26,36 @@
 from projectq.ops import QubitOperator
 
 
+def wigner_seitz_length_scale(wigner_seitz_radius, n_particles, dimension):
+    """Function to give length_scale associated with Wigner-Seitz radius.
+
+    Args:
+        wigner_seitz_radius (float): The radius per particle in atomic units.
+        n_particles (int): The number of particles in the simulation cell.
+        dimension (int): The dimension of the system.
+
+    Returns:
+        length_scale (float): The length scale for the simulation.
+
+    Raises:
+        ValueError: System dimension must be 1, 2 or 3.
+    """
+    if dimension == 1:
+        volume_per_particle = 2. * numpy.pi * wigner_seitz_radius
+        length_scale = volume_per_particle * float(n_particles)
+    elif dimension == 2:
+        volume_per_particle = numpy.pi * wigner_seitz_radius ** 2.
+        volume = volume_per_particle * float(n_particles)
+        length_scale = volume ** (1. / 2.)
+    elif dimension == 3:
+        volume_per_particle = (4. * numpy.pi / 3.) * wigner_seitz_radius ** 3.
+        volume = volume_per_particle * float(n_particles)
+        length_scale = volume ** (1. / 3.)
+    else:
+        raise ValueError('System dimension must be 1, 2 or 3.')
+    return length_scale
+
+
 def dual_basis_u_operator(grid, geometry, spinless):
     """Return the external potential operator in plane wave dual basis.
 

src/fermilib/utils/_plane_wave_hamiltonian_test.py

@@ -27,11 +27,32 @@
     plane_wave_u_operator,
     dual_basis_u_operator,
     jordan_wigner_dual_basis_hamiltonian,
+    wigner_seitz_length_scale,
 )
 
 
 class PlaneWaveHamiltonianTest(unittest.TestCase):
 
+    def test_wigner_seitz_radius(self):
+        wigner_seitz_radius = 3.17
+        n_particles = 20
+        one_d_test = wigner_seitz_length_scale(
+            wigner_seitz_radius, n_particles, 1)
+        self.assertAlmostEqual(
+            one_d_test, n_particles * 2. * numpy.pi * wigner_seitz_radius)
+        two_d_test = wigner_seitz_length_scale(
+            wigner_seitz_radius, n_particles, 2) ** 2.
+        self.assertAlmostEqual(
+            two_d_test, n_particles * numpy.pi * wigner_seitz_radius ** 2.)
+        three_d_test = wigner_seitz_length_scale(
+            wigner_seitz_radius, n_particles, 3) ** 3.
+        self.assertAlmostEqual(
+            three_d_test, n_particles * (4. * numpy.pi / 3. *
+                                         wigner_seitz_radius ** 3.))
+        with self.assertRaises(ValueError):
+            wigner_seitz_length_scale(
+                wigner_seitz_radius, n_particles, 4)
+
     def test_fourier_transform(self):
         grid = Grid(dimensions=1, scale=1.5, length=3)
         spinless_set = [True, False]