diff --git a/bindings/rascal/models/asemd.py b/bindings/rascal/models/asemd.py index ce34e9395..f7de7c65c 100644 --- a/bindings/rascal/models/asemd.py +++ b/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) diff --git a/bindings/rascal/models/genericmd.py b/bindings/rascal/models/genericmd.py index f61c8d205..7cf315599 100644 --- a/bindings/rascal/models/genericmd.py +++ b/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) diff --git a/bindings/rascal/neighbourlist/structure_manager.py b/bindings/rascal/neighbourlist/structure_manager.py index 97b5b0d89..3d415d3b4 100644 --- a/bindings/rascal/neighbourlist/structure_manager.py +++ b/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: diff --git a/reference_data/inputs/1-Propanol.xyz b/reference_data/inputs/1-Propanol.xyz index a506e90c9..61e85d0c9 100644 --- a/reference_data/inputs/1-Propanol.xyz +++ b/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 diff --git a/reference_data/inputs/2-Propanol.xyz b/reference_data/inputs/2-Propanol.xyz index 580410f43..f13926550 100644 --- a/reference_data/inputs/2-Propanol.xyz +++ b/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 diff --git a/reference_data/inputs/H2O.xyz b/reference_data/inputs/H2O.xyz index 869f14986..7edfad10f 100644 --- a/reference_data/inputs/H2O.xyz +++ b/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 diff --git a/reference_data/inputs/Methoxyethane.xyz b/reference_data/inputs/Methoxyethane.xyz index 949cc7850..ae2472a28 100644 --- a/reference_data/inputs/Methoxyethane.xyz +++ b/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 diff --git a/reference_data/inputs/N3.xyz b/reference_data/inputs/N3.xyz index cb293cccc..0c10f4e90 100644 --- a/reference_data/inputs/N3.xyz +++ b/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 diff --git a/reference_data/inputs/benzene.xyz b/reference_data/inputs/benzene.xyz index 433a6fb1c..42790fd90 100644 --- a/reference_data/inputs/benzene.xyz +++ b/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 diff --git a/reference_data/inputs/methane.xyz b/reference_data/inputs/methane.xyz index 5c4f4ae5e..1d4984a22 100644 --- a/reference_data/inputs/methane.xyz +++ b/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 diff --git a/reference_data/inputs/water_rotations.xyz b/reference_data/inputs/water_rotations.xyz index e46417691..6015eeba5 100644 --- a/reference_data/inputs/water_rotations.xyz +++ b/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 diff --git a/tests/python/python_structure_manager_test.py b/tests/python/python_structure_manager_test.py index 0e3d788af..de7c42766 100644 --- a/tests/python/python_structure_manager_test.py +++ b/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):