Merge pull request #12 from pyiron/ase

ASE wrapper

.binder/environment.yml

@@ -6,3 +6,4 @@ dependencies:
 - numpy =1.26.0
 - black =24.8.0
 - h5py =3.12.1
+- ase =3.23.0

.ci_support/environment.yml

@@ -6,3 +6,4 @@ dependencies:
 - numpy =1.26.0
 - black =24.8.0
 - h5py =3.12.1
+- ase =3.23.0

docs/environment.yml

@@ -12,3 +12,4 @@ dependencies:
 - numpy =1.26.0
 - black =24.8.0
 - h5py =3.12.1
+- ase =3.23.0

stinx/ase.py

@@ -0,0 +1,70 @@
+import numpy as np
+import scipy.constants as sc
+from stinx.input import sphinx
+from ase.io.vasp import _handle_ase_constraints
+
+
+def get_constraints(atoms):
+    """
+    Get the constraints of the atoms object. The constraints are returned as a
+    boolean array. True means the atom is movable, False means it is fixed.
+    """
+    if atoms.constraints:
+        return ~_handle_ase_constraints(atoms)
+    else:
+        return np.full(shape=atoms.positions.shape, fill_value=True)
+
+
+def _to_angstrom(cell, positions):
+    bohr_to_angstrom = sc.physical_constants["Bohr radius"][0] / sc.angstrom
+    cell = np.array(cell) / bohr_to_angstrom
+    positions = np.array(positions) / bohr_to_angstrom
+    return cell, positions
+
+
+def get_structure_group(structure, use_symmetry=True):
+    """
+    create a SPHInX Group object based on structure
+
+    Args:
+        structure (Atoms): ASE structure object
+        use_symmetry (bool): Whether or not consider internal symmetry
+
+    Returns:
+        (Group): structure group
+    """
+    cell, positions = _to_angstrom(structure.cell, structure.positions)
+    movable = get_constraints(structure)
+    labels = structure.get_initial_magnetic_moments()
+    elements = np.array(structure.get_chemical_symbols())
+    species = []
+    for elm_species in np.unique(elements):
+        elm_list = elements == elm_species
+        atom_list = []
+        for elm_pos, elm_magmom, selective in zip(
+            positions[elm_list],
+            labels[elm_list],
+            movable[elm_list],
+        ):
+            atom_group = {
+                "coords": np.array(elm_pos),
+                "label": f'"spin_{elm_magmom}"',
+            }
+            if all(selective):
+                atom_group["movable"] = True
+            elif any(selective):
+                for xx in np.array(["X", "Y", "Z"])[selective]:
+                    atom_group["movable" + xx] = True
+            atom_list.append(sphinx.structure.species.atom.create(**atom_group))
+        species.append(
+            sphinx.structure.species.create(element=f'"{elm_species}"', atom=atom_list)
+        )
+    symmetry = None
+    if not use_symmetry:
+        symmetry = sphinx.structure.symmetry.create(
+            operator=sphinx.structure.symmetry.operator.create(S=np.eye(3).tolist())
+        )
+    structure_group = sphinx.structure.create(
+        cell=np.array(cell), species=species, symmetry=symmetry
+    )
+    return structure_group

tests/unit/test_ase.py

@@ -0,0 +1,54 @@
+import unittest
+from ase.build import bulk
+from stinx.ase import get_structure_group
+from stinx.toolkit import to_sphinx
+import re
+from ase.constraints import FixedPlane
+
+
+class TestStinx(unittest.TestCase):
+    def test_Ni_Al_bulk(self):
+        structure = bulk("Al", cubic=True)
+        structure[0].symbol = "Ni"
+        input_dict = get_structure_group(structure, use_symmetry=False)
+        for key in ["cell", "species", "species___0", "symmetry"]:
+            self.assertTrue(key in input_dict)
+        text = to_sphinx(input_dict)
+        self.assertEqual(
+            len(re.findall(r"\bspecies\b", text, re.IGNORECASE)),
+            2,
+            msg="There must be exactly two species (Al and Ni) in the structure",
+        )
+        self.assertEqual(
+            len(re.findall(r"\batom\b", text, re.IGNORECASE)),
+            4,
+            msg="There must be exactly four atoms in cubic fcc",
+        )
+
+    def test_constraint_bulk(self):
+        structure = bulk("Al", cubic=True)
+        c = FixedPlane([0], [1, 0, 0])
+        structure.set_constraint(c)
+        struct_group = get_structure_group(structure)
+        self.assertFalse(
+            "movableX" in struct_group["species"]["atom"],
+            msg="Not allowed to move along X",
+        )
+        for term in ["movableY", "movableZ"]:
+            self.assertTrue(
+                term in struct_group["species"]["atom"],
+                msg="Must be allowed to move along Y and Z",
+            )
+        for term in ["movableX", "movableY", "movableZ"]:
+            self.assertFalse(
+                term in struct_group["species"]["atom___0"],
+                msg="Must be using the global movable",
+            )
+        self.assertTrue(
+            "movable" in struct_group["species"]["atom___0"],
+            msg="Must be globally movable",
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()