Automatically center box when converting to GSD.

src/lammpsio/snapshot.py

@@ -68,8 +68,9 @@ def from_hoomd_gsd(cls, frame):
         frame.validate()
 
         # process HOOMD box to LAMMPS box
-        L = frame.configuration.box[:3]
-        tilt = frame.configuration.box[3:]
+        box = numpy.array(frame.configuration.box, copy=True)
+        L = box[:3]
+        tilt = box[3:]
         if frame.configuration.dimensions == 3:
             tilt[0] *= L[1]
             tilt[1:] *= L[2]
@@ -78,7 +79,12 @@ def from_hoomd_gsd(cls, frame):
             # HOOMD boxes can have Lz = 0, but LAMMPS does not allow this.
             if L[2] == 0:
                 L[2] = 1.0
-        box = Box(low=-0.5 * L, high=0.5 * L, tilt=tilt)
+
+        matrix = numpy.array(
+            [[L[0], tilt[0], tilt[1]], [0, L[1], tilt[2]], [0, 0, L[2]]]
+        )
+        low = -0.5 * numpy.sum(matrix, axis=1)
+        box = Box.from_matrix(low=low, matrix=matrix)
 
         snap = Snapshot(
             N=frame.particles.N,
@@ -225,12 +231,11 @@ def to_hoomd_gsd(self, type_map=None):
         if self.step is not None:
             frame.configuration.step = int(self.step)
 
-        # we could shift the box later, but for now this is an error
-        if not numpy.allclose(-self.box.low, self.box.high):
-            raise ValueError("GSD boxes must be centered around 0")
         L = self.box.high - self.box.low
         if self.box.tilt is not None:
-            tilt = self.box.tilt
+            tilt = self.box.tilt.copy()
+            tilt[0] /= L[1]
+            tilt[1:] /= L[2]
         else:
             tilt = [0, 0, 0]
         frame.configuration.box = numpy.concatenate((L, tilt))
@@ -247,7 +252,16 @@ def to_hoomd_gsd(self, type_map=None):
 
         frame.particles.N = self.N
         if self.has_position():
-            frame.particles.position = self.position.copy()
+            # Center the positions using HOOMD tilt factors (computed above)
+            matrix = numpy.array(
+                [
+                    [L[0], tilt[0] * L[1], tilt[1] * L[2]],
+                    [0, L[1], tilt[2] * L[2]],
+                    [0, 0, L[2]],
+                ]
+            )
+            center = self.box.low + 0.5 * numpy.sum(matrix, axis=1)
+            frame.particles.position = self.position - center
         if self.has_velocity():
             frame.particles.velocity = self.velocity.copy()
         if self.has_image():

tests/test_snapshot.py

@@ -46,8 +46,8 @@ def test_gsd_conversion():
     # make Snapshot from GSD
     snap, type_map = lammpsio.Snapshot.from_hoomd_gsd(frame)
     assert snap.step == 3
-    assert numpy.allclose(snap.box.low, [-2, -2.5, -3])
-    assert numpy.allclose(snap.box.high, [2, 2.5, 3])
+    assert numpy.allclose(snap.box.low, [-2.85, -3.4, -3])
+    assert numpy.allclose(snap.box.high, [1.15, 1.6, 3])
     assert numpy.allclose(snap.box.tilt, [0.5, 1.2, 1.8])
     assert snap.N == 2
     assert numpy.allclose(snap.position, [[0.1, 0.2, 0.3], [-0.1, -0.2, -0.3]])
@@ -99,10 +99,22 @@ def test_gsd_conversion():
     snap2.to_hoomd_gsd()
     assert numpy.all(snap2.id == [2, 1])
 
-    # check for error out on bad box
-    snap2.box.low = [-10, -10, -10]
-    with pytest.raises(ValueError):
-        snap2.to_hoomd_gsd()
+    # check that orthorhombic box is set correctly when not originally centered
+    snap3 = lammpsio.Snapshot(N=2, box=lammpsio.Box([0, 0, 0], [15, 15, 15]))
+    snap3.position = [[0, 0, 0], [1, 1, 1]]
+    # center calculated using box vectors as describe in LAMMPS documentation
+    center = [7.5, 7.5, 7.5]
+    frame5 = snap3.to_hoomd_gsd()
+    assert numpy.allclose(frame5.configuration.box, [15, 15, 15, 0, 0, 0])
+    assert numpy.allclose(frame5.particles.position, snap3.position - center)
+
+    # check that triclinic box is set correctly when not originally centered
+    snap3.box.tilt = [1.5, 3.0, 4.5]
+    # center calculated using box vectors as describe in LAMMPS documentation
+    center = [9.75, 9.75, 7.5]
+    frame6 = snap3.to_hoomd_gsd()
+    assert numpy.allclose(frame6.configuration.box, [15, 15, 15, 0.1, 0.2, 0.3])
+    assert numpy.allclose(frame6.particles.position, snap3.position - center)
 
 
 @pytest.mark.skipif(not has_gsd, reason="gsd not installed")