Merge pull request #400 from pyiron/step_surface

High index surfaces

pyiron_atomistics/atomistics/structure/factory.py

@@ -3,26 +3,26 @@
 # Distributed under the terms of "New BSD License", see the LICENSE file.
 
 from ase.build import (
-        add_adsorbate,
-        add_vacuum,
-        bcc100,
-        bcc110,
-        bcc111,
-        diamond100,
-        diamond111,
-        fcc100,
-        fcc110,
-        fcc111,
-        fcc211,
-        hcp0001,
-        hcp10m10,
-        mx2,
-        hcp0001_root,
-        fcc111_root,
-        bcc111_root,
-        root_surface,
-        root_surface_analysis,
-        surface as ase_surf
+    add_adsorbate,
+    add_vacuum,
+    bcc100,
+    bcc110,
+    bcc111,
+    diamond100,
+    diamond111,
+    fcc100,
+    fcc110,
+    fcc111,
+    fcc211,
+    hcp0001,
+    hcp10m10,
+    mx2,
+    hcp0001_root,
+    fcc111_root,
+    bcc111_root,
+    root_surface,
+    root_surface_analysis,
+    surface as ase_surf
 )
 import numpy as np
 from pyiron_atomistics.atomistics.structure.factories.ase import AseFactory
@@ -31,7 +31,8 @@
 from pyiron_atomistics.atomistics.structure.factories.compound import CompoundFactory
 from pyiron_atomistics.atomistics.structure.pyironase import publication as publication_ase
 from pyiron_atomistics.atomistics.structure.atoms import CrystalStructure, Atoms, \
-    ase_to_pyiron, pymatgen_to_pyiron, ovito_to_pyiron
+    ase_to_pyiron, pymatgen_to_pyiron, ovito_to_pyiron, pyiron_to_pymatgen
+from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 from pyiron_atomistics.atomistics.structure.periodic_table import PeriodicTable
 from pyiron_base import state, PyironFactory, deprecate
 import types
@@ -76,24 +77,29 @@ def compound(self):
 
     def cut(self, *args, **kwargs):
         return self.ase.cut(*args, **kwargs)
+
     cut.__doc__ = AseFactory.cut.__doc__
 
     def stack(self, *args, **kwargs):
         return self.ase.stack(*args, **kwargs)
+
     stack.__doc__ = AseFactory.stack.__doc__
 
     def read(self, *args, **kwargs):
         return self.ase.read(*args, **kwargs)
+
     read.__doc__ = AseFactory.read.__doc__
 
     @deprecate(message="Please use .read or .ase.read", version="0.2.2")
     def ase_read(self, *args, **kwargs):
         return self.ase.read(*args, **kwargs)
+
     ase_read.__doc__ = AseFactory.read.__doc__
 
     @deprecate(message="Please use .bulk or .ase.bulk", version="0.2.2")
     def ase_bulk(self, *args, **kwargs):
         return self.ase.bulk(*args, **kwargs)
+
     ase_bulk.__doc__ = AseFactory.bulk.__doc__
 
     def bulk(
@@ -137,7 +143,7 @@ def bulk(
 
     @staticmethod
     def surface(
-        element, surface_type, size=(1, 1, 1), vacuum=1.0, center=False, pbc=True, **kwargs
+            element, surface_type, size=(1, 1, 1), vacuum=1.0, center=False, pbc=True, **kwargs
     ):
         """
         Generate a surface based on the ase.build.surface module.
@@ -225,7 +231,7 @@ def surface_hkl(lattice, hkl, layers, vacuum=1.0, center=False, pbc=True):
         z_max = np.max(surface.positions[:, 2])
         surface.cell[2, 2] = z_max + vacuum
         if center:
-            surface.positions += 0.5 * surface.cell[2] - [0, 0, z_max/2]
+            surface.positions += 0.5 * surface.cell[2] - [0, 0, z_max / 2]
         surface.pbc = pbc
         return ase_to_pyiron(surface)
 
@@ -374,6 +380,7 @@ def element(parent_element, new_element_name=None, spin=None, potential_file=Non
     @deprecate(message="Use .aimsgb.info", version="0.2.2")
     def aimsgb_info(self, axis, max_sigma):
         return self.aimsgb.info(axis=axis, max_sigma=max_sigma)
+
     aimsgb_info.__doc__ = AimsgbFactory.info.__doc__
 
     @deprecate(message="Use .aimsgb.build", version="0.2.2")
@@ -396,6 +403,7 @@ def aimsgb_build(
             delete_layer=delete_layer,
             add_if_dist=add_if_dist
         )
+
     aimsgb_build.__doc__ = AimsgbFactory.build.__doc__
 
     @staticmethod
@@ -412,3 +420,104 @@ def from_pymatgen(pymatgen_obj):
     @wraps(ovito_to_pyiron)
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
+
+    def high_index_surface_info(self, element, crystal_structure, lattice_constant,
+                                terrace_orientation=None, step_orientation=None, kink_orientation=None,
+                                step_down_vector=None, length_step=3, length_terrace=3, length_kink=1):
+        """
+        Gives the miller indices of high index surface required to create a stepped and kink surface, based
+        on the general orientation and length of terrace, step and kinks respectively. The microfacet notation used is
+        based on the work of Van Hove et al.,[1].
+
+        [1] Van Hove, M. A., and G. A. Somorjai. "A new microfacet notation for high-Miller-index surfaces of cubic
+        materials with terrace, step and kink structures." Surface Science 92.2-3 (1980): 489-518.
+
+        Args:
+            element (str): The parent element eq. "N", "O", "Mg" etc.
+            crystal_structure (str): The crystal structure of the lattice
+            lattice_constant (float): The lattice constant
+            terrace_orientation (list): The miller index of the terrace eg., [1,1,1]
+            step_orientation (list): The miller index of the step eg., [1,1,0]
+            kink_orientation (list): The miller index of the kink eg., [1,1,1]
+            step_down_vector (list): The direction for stepping down from the step to next terrace eg., [1,1,0]
+            length_terrace (int): The length of the terrace along the kink direction in atoms eg., 3
+            length_step (int): The length of the step along the step direction in atoms eg., 3
+            length_kink (int): The length of the kink along the kink direction in atoms eg., 1
+
+
+        Returns:
+            high_index_surface: The high miller index surface which can be used to create slabs
+            fin_kink_orientation: The kink orientation lying in the terrace
+            fin_step_orientation: The step orientation lying in the terrace
+        """
+        terrace_orientation = terrace_orientation if terrace_orientation is not None else [1, 1, 1]
+        step_orientation = step_orientation if step_orientation is not None else [1, 1, 0]
+        kink_orientation = kink_orientation if kink_orientation is not None else [1, 1, 1]
+        step_down_vector = step_down_vector if step_down_vector is not None else [1, 1, 0]
+
+        basis = self.crystal(element=element, bravais_basis=crystal_structure, lattice_constant=lattice_constant)
+        sym = basis.get_symmetry()
+        eqvdirs = np.unique(np.matmul(sym.rotations[:], (np.array(step_orientation))), axis=0)
+        eqvdirk = np.unique(np.matmul(sym.rotations[:], (np.array(kink_orientation))), axis=0)
+        eqvdirs_ind = np.where(np.dot(np.squeeze(eqvdirs), terrace_orientation) == 0)[0]
+        eqvdirk_ind = np.where(np.dot(np.squeeze(eqvdirk), terrace_orientation) == 0)[0]
+        if len(eqvdirs_ind) == 0:
+            raise ValueError('Step orientation vector should lie in terrace.\
+            For the given choice I could not find any symmetrically equivalent vector that lies in the terrace.\
+            please change the stepOrientation and try again')
+        if len(eqvdirk_ind) == 0:
+            raise ValueError('Kink orientation vector should lie in terrace.\
+            For the given choice I could not find any symmetrically equivalent vector that lies in the terrace.\
+            please change the kinkOrientation and try again')
+        temp = (np.cross(np.squeeze(eqvdirk[eqvdirk_ind[0]]), np.squeeze(eqvdirs))).tolist().index(terrace_orientation)
+        fin_kink_orientation = eqvdirk[eqvdirk_ind[0]]
+        fin_step_orientation = eqvdirs[temp]
+        vec1 = (np.asanyarray(fin_step_orientation).dot(length_step)) + \
+               (np.asanyarray(fin_kink_orientation).dot(length_kink))
+        vec2 = (np.asanyarray(fin_kink_orientation).dot(length_terrace)) + step_down_vector
+        high_index_surface = np.cross(np.asanyarray(vec1), np.asanyarray(vec2))
+        high_index_surface = np.array(high_index_surface / np.gcd.reduce(high_index_surface), dtype=int)
+
+        return high_index_surface, fin_kink_orientation, fin_step_orientation
+
+    def high_index_surface(self, element, crystal_structure, lattice_constant,
+                           terrace_orientation=None, step_orientation=None, kink_orientation=None,
+                           step_down_vector=None, length_step=3, length_terrace=3, length_kink=1, layers=60,
+                           vacuum=10):
+        """
+        Gives a slab positioned at the bottom with the high index surface computed by high_index_surface_info().
+        Args:
+            element (str): The parent element eq. "N", "O", "Mg" etc.
+            crystal_structure (str): The crystal structure of the lattice
+            lattice_constant (float): The lattice constant
+            terrace_orientation (list): The miller index of the terrace. default: [1,1,1]
+            step_orientation (list): The miller index of the step. default: [1,1,0]
+            kink_orientation (list): The miller index of the kink. default: [1,1,1]
+            step_down_vector (list): The direction for stepping down from the step to next terrace. default: [1,1,0]
+            length_terrace (int): The length of the terrace along the kink direction in atoms. default: 3
+            length_step (int): The length of the step along the step direction in atoms. default: 3
+            length_kink (int): The length of the kink along the kink direction in atoms. default: 1
+            layers (int): Number of layers of the high_index_surface. default: 60
+            vacuum (float): Thickness of vacuum on the top of the slab. default:10
+
+        Returns:
+            slab: pyiron_atomistics.atomistics.structure.atoms.Atoms instance Required surface
+        """
+        basis = self.crystal(element=element, bravais_basis=crystal_structure, lattice_constant=lattice_constant)
+        high_index_surface, _, _ = self.high_index_surface_info(element=element,
+                                                                crystal_structure=crystal_structure,
+                                                                lattice_constant=lattice_constant,
+                                                                terrace_orientation=terrace_orientation,
+                                                                step_orientation=step_orientation,
+                                                                kink_orientation=kink_orientation,
+                                                                step_down_vector=step_down_vector,
+                                                                length_step=length_step,
+                                                                length_terrace=length_terrace,
+                                                                length_kink=length_kink)
+        surf = ase_surf(basis, high_index_surface, layers, vacuum)
+        sga = SpacegroupAnalyzer(pyiron_to_pymatgen(ase_to_pyiron(surf)))
+        pmg_refined = sga.get_refined_structure()
+        slab = pymatgen_to_pyiron(pmg_refined)
+        slab.positions[:, 2] = slab.positions[:, 2] - np.min(slab.positions[:, 2])
+        slab.set_pbc(True)
+        return slab

tests/atomistics/structure/test_high_index_surface.py

@@ -0,0 +1,44 @@
+import unittest
+from pyiron_atomistics.atomistics.structure.factory import StructureFactory
+
+
+class TestHighIndexSurface(unittest.TestCase):
+    def test_high_index_surface(self):
+        slab = StructureFactory().high_index_surface(element='Ni', crystal_structure='fcc',
+                                                     lattice_constant=3.526,
+                                                     terrace_orientation=[1, 1, 1],
+                                                     step_orientation=[1, 1, 0],
+                                                     kink_orientation=[1, 0, 1],
+                                                     step_down_vector=[1, 1, 0], length_step=2,
+                                                     length_terrace=3,
+                                                     length_kink=1, layers=60,
+                                                     vacuum=10)
+
+        self.assertEqual(len(slab), 60)
+
+    def test_high_index_surface_info(self):
+        h, s, k = StructureFactory().high_index_surface_info(element='Ni', crystal_structure='fcc',
+                                                             lattice_constant=3.526,
+                                                             terrace_orientation=[1, 1, 1],
+                                                             step_orientation=[1, 1, 0],
+                                                             kink_orientation=[1, 0, 1],
+                                                             step_down_vector=[1, 1, 0], length_step=2,
+                                                             length_terrace=3,
+                                                             length_kink=1)
+        self.assertEqual(len(h), 3)
+        self.assertEqual(h[0], -9)
+        self.assertEqual(len(k), 3)
+        self.assertEqual(len(s), 3)
+        with self.assertRaises(ValueError):
+            StructureFactory().high_index_surface_info(element='Ni', crystal_structure='fcc',
+                                                       lattice_constant=3.526,
+                                                       terrace_orientation=[1, 1, 1],
+                                                       step_orientation=[1, 0, 0],
+                                                       kink_orientation=[1, 0, 0],
+                                                       step_down_vector=[1, 1, 0], length_step=2,
+                                                       length_terrace=3,
+                                                       length_kink=1)
+
+
+if __name__ == '__main__':
+    unittest.main()