structure sanitation

bindings/rascal/models/asemd.py

@@ -39,13 +39,14 @@ def calculate(
     ):
         Calculator.calculate(self, atoms, properties, system_changes)
 
+        self.atoms.wrap(eps=1e-11)
+
         if self.manager is None:
-            #  happens at the begining of the MD run
+            # happens at the begining of the MD run
             at = self.atoms.copy()
-            at.wrap(eps=1e-11)
             self.manager = [at]
         elif isinstance(self.manager, AtomsList):
-            structure = unpack_ase(self.atoms, wrap_pos=True)
+            structure = unpack_ase(self.atoms)
             structure.pop("center_atoms_mask")
             self.manager[0].update(**structure)
 

bindings/rascal/models/genericmd.py

@@ -106,15 +106,15 @@ def calculate(self, positions, cell_matrix):
         # re-wrapping of the atoms that needs to take place)
         self.atoms.set_cell(cell_matrix)
         self.atoms.set_positions(positions)
-
+        self.atoms.wrap(eps=1e-11)
+
         # Convert from ASE to librascal
         if self.manager is None:
             #  happens at the begining of the MD run
             at = self.atoms.copy()
-            at.wrap(eps=1e-11)
             self.manager = [at]
         elif isinstance(self.manager, AtomsList):
-            structure = unpack_ase(self.atoms, wrap_pos=True)
+            structure = unpack_ase(self.atoms)
             structure.pop("center_atoms_mask")
             self.manager[0].update(**structure)
 

bindings/rascal/neighbourlist/structure_manager.py

@@ -1,6 +1,7 @@
 from collections.abc import Iterable
 from functools import reduce
 from operator import and_
+from copy import deepcopy
 
 import numpy as np
 
@@ -248,50 +249,111 @@ def convert_to_structure_list(frames):
     structure_list = neighbour_list.AtomicStructureList()
     for frame in frames:
         if is_valid_structure(frame):
-            structure = frame
+            structure = sanitize_non_periodic_structure(frame)
+            structure = wrap_structure(structure)
         else:
             if is_ase_Atoms(frame):
+                frame.wrap(eps=1e-11)
+                frame = sanitize_non_periodic_ase(frame)
                 structure = unpack_ase(frame)
             else:
                 raise RuntimeError(
                     "Cannot convert structure of type {}".format(type(frame))
                 )
 
-        structure = sanitize_non_periodic_structure(structure)
         structure_list.append(**structure)
     return structure_list
 
 
+def convert_structure_to_ase(structure):
+    """Convert a structure as per is_valid_structure into an ASE Atoms object"""
+    from ase import Atoms
+
+    return Atoms(
+        positions=structure["positions"].T,
+        cell=structure["cell"].T,
+        numbers=structure["atom_types"].flatten(),
+        pbc=np.array(structure["pbc"].flatten(), dtype=bool),
+    )
+
+
+def sanitize_non_periodic_ase(frame):
+    """
+    Rascal expects a unit cell that contains all the atoms even if the
+    structure is not periodic.
+    In this case a unit cell is set to be the smallest rectangle that contains
+    all of the atoms and that starts at the origin. The atoms are moved inside
+    this box.
+
+    Parameters
+    ----------
+    frame : ase.Atoms
+
+
+    Returns
+    -------
+    frame : ase.Atoms
+        cell and positions have been modified if structure is not periodic
+    """
+
+    if not np.all(frame.get_pbc()):
+        frame = deepcopy(frame)
+        pos = frame.get_positions()
+        bounds = np.array([pos.min(axis=0), pos.max(axis=0)])
+        bounding_box_lengths = (bounds[1] - bounds[0]) * 1.1
+        new_cell = np.diag(bounding_box_lengths)
+        frame.set_cell(new_cell)
+        frame.center()
+    return frame
+
+
 def sanitize_non_periodic_structure(structure):
     """
     Rascal expects a unit cell that contains all the atoms even if the
     structure is not periodic.
-    If the cell is set to 0 and the structure is not periodic then it
-    is adapted to contain the atoms and the atoms are shifted inside the unit
-    cell.
+    In this case a unit cell is set to be the smallest rectangle that contains
+    all of the atoms and that starts at the origin. The atoms are moved inside
+    this box.
 
     Parameters
     ----------
-    structure : a valid structure as per is_valid_structure
+    structure : dict
+         a valid structure as per is_valid_structure
 
 
     Returns
     -------
-    a valid structure as per is_valid_structure
+    a valid structure as per is_valid_structure: dict
         cell and positions have been modified if structure is not periodic
     """
 
     if np.all(structure["pbc"] == 0):
-        cell = structure["cell"]
-        if np.allclose(cell, np.zeros((3, 3))):
-            pos = structure["positions"]
-            bounds = np.array([pos.min(axis=1), pos.max(axis=1)])
-            bounding_box_lengths = (bounds[1] - bounds[0]) * 1.05
-            new_cell = np.diag(bounding_box_lengths)
-            CoM = pos.mean(axis=1)
-            disp = 0.5 * bounding_box_lengths - CoM
-            new_pos = pos + disp[:, None]
-            structure["positions"] = new_pos
+        atoms = convert_structure_to_ase(structure)
+        atoms = sanitize_non_periodic_ase(atoms)
+        structure = unpack_ase(atoms)
+    return structure
+
+
+def wrap_structure(structure, tol=1e-11):
+    """
+    Wrap the atoms inside the unit cell for periodic structure. This a wrapper
+    around `ase.geometry.wrap_positions`.
+
+    Parameters
+    ----------
+    structure : a valid structure as per is_valid_structure
+
+
+    Returns
+    -------
+    a valid structure as per is_valid_structure
+        cell and positions have been modified if structure is not periodic
+    """
+    positions = structure["positions"].T
+    cell = structure["cell"].T
+    pbc = np.array(structure["pbc"], dtype=bool).flatten()
+    positions = wrap_positions(positions, cell, pbc, eps=tol)
+    structure["positions"] = np.array(positions.T, order="F")
     return structure
 
 
@@ -308,7 +370,7 @@ def is_ase_Atoms(frame):
     return is_ase
 
 
-def unpack_ase(frame, wrap_pos=False):
+def unpack_ase(frame):
     """
     Convert ASE Atoms object to rascal's equivalent
 
@@ -331,9 +393,6 @@ def unpack_ase(frame, wrap_pos=False):
     numbers = frame.get_atomic_numbers()
     pbc = frame.get_pbc().astype(int)
 
-    if wrap_pos:
-        positions = wrap_positions(positions, cell, frame.get_pbc(), eps=1e-11)
-
     if "center_atoms_mask" in frame.arrays.keys():
         center_atoms_mask = frame.get_array("center_atoms_mask")
     else:

reference_data/inputs/1-Propanol.xyz

@@ -1,5 +1,5 @@
 12
-Lattice="10.0 0.0 0.0 0.0 10.0 0.0 0.0 0.0 10.0" Properties=species:S:1:pos:R:3 309=T #=T cellangstrom=T 14.46900=T Bead:=T positionsangstrom=T CELLabcABC:=T 0=T 99.79000=T 90.00000=T 6.54200=T 13.41500=T Step:=T pbc="T T T"
+Lattice="10.0 0.0 0.0 0.0 10.0 0.0 0.0 0.0 10.0" Properties=species:S:1:pos:R:3 pbc="F F F"
 C       0.67522167       0.84488083       0.07696250
 C       0.10989167      -0.55559917      -0.16835750
 C      -1.41576833      -0.54222917      -0.17510750

reference_data/inputs/2-Propanol.xyz

@@ -1,5 +1,5 @@
 12
-Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 309=T #=T cellangstrom=T 14.46900=T Bead:=T positionsangstrom=T CELLabcABC:=T 0=T 99.79000=T 90.00000=T 6.54200=T 13.41500=T Step:=T pbc="T T T"
+Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 pbc="F F F"
 O       0.34376458       0.23988583       1.61715000
 H       1.16430458      -0.26240417       1.89935000
 C      -0.10618542      -0.29654417       0.33695000

reference_data/inputs/H2O.xyz

@@ -1,5 +1,5 @@
 3
-Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 #=T 19.28600=T cellangstrom=T Bead:=T positionsangstrom=T CELLabcABC:=T 0=T 238=T 90.00000=T 19.04200=T 118.72000=T Step:=T 3.56990=T pbc="T T T"
+Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 pbc="F F F"
 O      37.48522999      37.54021008      37.45053162
 H      37.93953999      36.70175408      37.34403162
 H      37.29488999      37.66003053      38.44113162

reference_data/inputs/Methoxyethane.xyz

@@ -1,5 +1,5 @@
 12
-Lattice="10.0 0.0 0.0 0.0 10.0 0.0 0.0 0.0 10.0" Properties=species:S:1:pos:R:3 309=T #=T cellangstrom=T 14.46900=T Bead:=T positionsangstrom=T CELLabcABC:=T 0=T 99.79000=T 90.00000=T 6.54200=T 13.41500=T Step:=T pbc="T T T"
+Lattice="10.0 0.0 0.0 0.0 10.0 0.0 0.0 0.0 10.0" Properties=species:S:1:pos:R:3 309=T pbc="F F F"
 O      -0.04901667       0.86314167       0.00512500
 C      -0.02961667      -0.56575833      -0.00277500
 C      -1.46501667      -1.09525833       0.00812500

reference_data/inputs/N3.xyz

@@ -1,5 +1,5 @@
 3
-Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 #=T cellangstrom=T 7.76940=T 0=T 13.86620=T Bead:=T positionsangstrom=T CELLabcABC:=T 8.24590=T 99.72100=T 90.00000=T 184=T Step:=T pbc="T T T"
+Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 pbc="F F F"
 N      38.11170167      38.22351667      36.78380000
 N      37.50339267      37.49964667      37.49280000
 N      36.88490567      36.77683667      38.22340000

reference_data/inputs/benzene.xyz

@@ -1,14 +1,14 @@
 12
-Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 #=T 90.86800=T 13.26960=T 108.09300=T 0=T 91.74300=T Bead:=T positionsangstrom=T CELLabcABC:=T 121=T 6.74540=T cellangstrom=T Step:=T 6.48750=T pbc="T T T"
-C      37.44163263      36.44837987      36.57577912 
-H      37.38995263      35.63112987      35.85770912 
-C      37.41823263      36.17213987      37.94639912 
-H      37.32973263      35.14457987      38.29609912 
-C      37.48350263      37.22622987      38.86999912 
-H      37.48745263      37.00926987      39.93869912 
-C      37.55605263      38.55031987      38.42649912 
-H      37.60366263      39.36741987      39.14629912 
-C      37.57583263      38.82755987      37.05138912 
-H      37.65543263      39.85797987      36.70560912 
-C      37.52267263      37.77536987      36.12947912 
-H      37.55848263      37.98962987      35.06100912 
+Lattice="75.0 0.0 0.0 0.0 75.0 0.0 0.0 0.0 75.0" Properties=species:S:1:pos:R:3 pbc="F F F"
+C      37.44163263      36.44837987      36.57577912
+H      37.38995263      35.63112987      35.85770912
+C      37.41823263      36.17213987      37.94639912
+H      37.32973263      35.14457987      38.29609912
+C      37.48350263      37.22622987      38.86999912
+H      37.48745263      37.00926987      39.93869912
+C      37.55605263      38.55031987      38.42649912
+H      37.60366263      39.36741987      39.14629912
+C      37.57583263      38.82755987      37.05138912
+H      37.65543263      39.85797987      36.70560912
+C      37.52267263      37.77536987      36.12947912
+H      37.55848263      37.98962987      35.06100912

reference_data/inputs/methane.xyz

@@ -1,7 +1,7 @@
-           5
- Lattice="100 0 0 0 100 0 0 0 100" pbc="F F F" Properties=species:S:1:pos:R:3
- C   0.0 0.0 0.0
- H   0.629118 0.629118 0.629118
- H   -0.629118 -0.629118 0.629118
- H   0.629118 -0.629118 -0.629118 
- H   -0.629118 0.629118 -0.629118
+5
+Lattice="100 0 0 0 100 0 0 0 100" Properties=species:S:1:pos:R:3 pbc="F F F"
+C   0.0 0.0 0.0
+H   0.629118 0.629118 0.629118
+H   -0.629118 -0.629118 0.629118
+H   0.629118 -0.629118 -0.629118
+H   -0.629118 0.629118 -0.629118

reference_data/inputs/water_rotations.xyz

@@ -1,10 +1,10 @@
 3
-Properties=species:S:1:pos:R:3 Lattice="100 0 0 0 100 0 0 0 100" 
-O      2.09821000       0.92801900       1.89574000 
-H      1.30969000       1.64896000       1.48342000 
-H      1.84956000       0.10002900       1.20268000 
+Properties=species:S:1:pos:R:3 Lattice="100 0 0 0 100 0 0 0 100" pbc="F F F"
+O      2.09821000       0.92801900       1.89574000
+H      1.30969000       1.64896000       1.48342000
+H      1.84956000       0.10002900       1.20268000
 3
-Properties=species:S:1:pos:R:3 Lattice="100 0 0 0 100 0 0 0 100" 
+Properties=species:S:1:pos:R:3 Lattice="100 0 0 0 100 0 0 0 100" pbc="F F F"
 O      2.17346042  1.10841022  1.70798335
 H      1.15970281  1.52033811  1.37013956
 H      1.92429677  0.04825967  1.50371709

tests/python/python_structure_manager_test.py

@@ -1,4 +1,4 @@
-import ase.io
+from ase.io import read
 from ase.build import molecule
 from rascal.neighbourlist import get_neighbourlist, AtomsList
 from rascal.neighbourlist.structure_manager import (
@@ -164,6 +164,32 @@ def test_manager_iteration_2(self):
                     dist = np.linalg.norm(neigh.position - center.position)
 
 
+class TestNLsanitation(unittest.TestCase):
+    def setUp(self):
+        """
+        test that improper structures are sanitized by the
+        """
+
+        self.frames = []
+        # methane.xyz is not inside the unit cell and is not periodic
+        # dummy_structure.json is periodic but not inside the unitcell
+        fns = ["methane.xyz", "dummy_structure.json"]
+        for fn in fns:
+            self.frames += [read(os.path.join(inputs_path, fn))]
+
+        self.cutoff = 3.0
+        self.nl_options = [
+            dict(name="centers", args=dict()),
+            dict(name="neighbourlist", args=dict(cutoff=self.cutoff)),
+            dict(name="strict", args=dict(cutoff=self.cutoff)),
+        ]
+
+    def test_sanitation(self):
+        for frame in self.frames:
+            # it should not raise errors
+            _ = AtomsList(frame, self.nl_options)
+
+
 class TestNLStrict(unittest.TestCase):
     def setUp(self):
         """
@@ -318,7 +344,7 @@ class CenterSelectTest(unittest.TestCase):
 
     def setUp(self):
         filename = "reference_data/inputs/small_molecule.json"
-        self.frame = ase.io.read(filename)
+        self.frame = read(filename)
         self.natoms = len(self.frame)
 
     def get_mask(self):