Generalized wigner_seitz_length_scale and corrected error in code and…

… test. (#88)

* Generalized wigner_seitz_length_scale and corrected error in code and test.

* Fixed truediv

src/fermilib/utils/_plane_wave_hamiltonian.py

@@ -38,21 +38,25 @@ def wigner_seitz_length_scale(wigner_seitz_radius, n_particles, dimension):
         length_scale (float): The length scale for the simulation.
 
     Raises:
-        ValueError: System dimension must be 1, 2 or 3.
+        ValueError: System dimension must be a positive integer.
     """
-    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.)
+    if not isinstance(dimension, int) or dimension < 1:
+        raise ValueError('System dimension must be a positive integer.')
+
+    half_dimension = dimension // 2
+    if dimension % 2:
+        volume_per_particle = (2 * numpy.math.factorial(half_dimension) *
+                               (4 * numpy.pi) ** half_dimension /
+                               numpy.math.factorial(dimension) *
+                               wigner_seitz_radius ** dimension)
     else:
-        raise ValueError('System dimension must be 1, 2 or 3.')
+        volume_per_particle = (numpy.pi ** half_dimension /
+                               numpy.math.factorial(half_dimension) *
+                               wigner_seitz_radius ** dimension)
+
+    volume = volume_per_particle * n_particles
+    length_scale = volume ** (1. / dimension)
+
     return length_scale
 
 

src/fermilib/utils/_plane_wave_hamiltonian_test.py

@@ -33,25 +33,47 @@
 
 class PlaneWaveHamiltonianTest(unittest.TestCase):
 
-    def test_wigner_seitz_radius(self):
+    def test_wigner_seitz_radius_1d(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)
+            one_d_test, n_particles * 2. * wigner_seitz_radius)
+
+    def test_wigner_seitz_radius_2d(self):
+        wigner_seitz_radius = 0.5
+        n_particles = 3
         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.)
+
+    def test_wigner_seitz_radius_3d(self):
+        wigner_seitz_radius = 4.6
+        n_particles = 37
         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.))
+
+    def test_wigner_seitz_radius_6d(self):
+        wigner_seitz_radius = 5.
+        n_particles = 42
+        six_d_test = wigner_seitz_length_scale(
+            wigner_seitz_radius, n_particles, 6) ** 6
+        self.assertAlmostEqual(
+            six_d_test, n_particles * (numpy.pi ** 3 / 6 *
+                                       wigner_seitz_radius ** 6))
+
+    def test_wigner_seitz_radius_bad_dimension_not_integer(self):
+        with self.assertRaises(ValueError):
+            wigner_seitz_length_scale(3, 2, dimension=4.2)
+
+    def test_wigner_seitz_radius_bad_dimension_not_positive(self):
         with self.assertRaises(ValueError):
-            wigner_seitz_length_scale(
-                wigner_seitz_radius, n_particles, 4)
+            wigner_seitz_length_scale(3, 2, dimension=0)     
 
     def test_fourier_transform(self):
         grid = Grid(dimensions=1, scale=1.5, length=3)