From e0bf6ca83f43b55a097a0d04a030f8f7d16a6eb0 Mon Sep 17 00:00:00 2001
From: skatnagallu <85609781+skatnagallu@users.noreply.github.com>
Date: Thu, 28 Oct 2021 11:32:38 +0200
Subject: [PATCH 01/14] Update factory.py

A new function to get a high index surface based on the step length, step orientation and kink orientation.
---
 .../atomistics/structure/factory.py           | 35 +++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index eade8aa1a..51d768686 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -407,3 +407,38 @@ def from_pymatgen(pymatgen_obj):
     @wraps(ovito_to_pyiron)
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
+     
+    @staticmethod
+    def get_high_ind_surf_sk(pr,element='Ni',Xstr='fcc',lattice_constant=3.526,terraceOrientation=[1,1,1],\
+              stepOrientation=[1,1,0],\
+              kinkOrientation=[1,0,1],stepDownVector=[1,1,0],\
+              lengthStep=3,lengthTerrace=3, lengthKink=1):
+                
+                
+
+            basis=pr.create_structure(element=element, bravais_basis=Xstr, lattice_constant=lattice_constant)
+            sym = basis.get_symmetry()
+            eqvdirS = np.unique(np.matmul(sym.rotations[:] , (np.array(stepOrientation))),axis=0)
+            eqvdirK = np.unique(np.matmul(sym.rotations[:] , (np.array(kinkOrientation))),axis=0)
+            eqvdirS_ind = np.where(np.dot(np.squeeze(eqvdirS),terraceOrientation)==0)[0]
+            eqvdirK_ind = np.where(np.dot(np.squeeze(eqvdirK),terraceOrientation)==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(terraceOrientation)
+            fin_kinkOrientation = eqvdirK[eqvdirK_ind[0]]
+            fin_stepOrientation = eqvdirS[temp]
+            vec1 = (np.asanyarray(fin_stepOrientation).dot(lengthStep)) + (np.asanyarray(fin_kinkOrientation).dot(lengthKink))
+            vec2 = (np.asanyarray(fin_kinkOrientation).dot(lengthTerrace)) + stepDownVector
+            highIndexSurface = np.cross(np.asanyarray(vec1),np.asanyarray(vec2))
+            highIndexSurface = np.array(highIndexSurface/np.gcd.reduce(highIndexSurface),dtype=int)
+
+        return fin_kinkOrientation, fin_stepOrientation, highIndexSurface

From 96e3ecd3f8b6c412641b09de5cf0cbd643bb3aad Mon Sep 17 00:00:00 2001
From: skatnagallu <85609781+skatnagallu@users.noreply.github.com>
Date: Thu, 28 Oct 2021 11:46:39 +0200
Subject: [PATCH 02/14] Indentation

Updating the indentations
---
 pyiron_atomistics/atomistics/structure/factory.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index 51d768686..ee223be95 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -441,4 +441,4 @@ def get_high_ind_surf_sk(pr,element='Ni',Xstr='fcc',lattice_constant=3.526,terra
             highIndexSurface = np.cross(np.asanyarray(vec1),np.asanyarray(vec2))
             highIndexSurface = np.array(highIndexSurface/np.gcd.reduce(highIndexSurface),dtype=int)
 
-        return fin_kinkOrientation, fin_stepOrientation, highIndexSurface
+         return fin_kinkOrientation, fin_stepOrientation, highIndexSurface

From 66d35ed69c63c070e8f070d498b6796caf771f35 Mon Sep 17 00:00:00 2001
From: sudarsan1989 <surendralal@mpie.de>
Date: Thu, 28 Oct 2021 16:46:07 +0200
Subject: [PATCH 03/14] :bug: Fix formatting

---
 .../atomistics/structure/factory.py           | 61 +++++++++----------
 1 file changed, 28 insertions(+), 33 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index ee223be95..d421d47be 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -409,36 +409,31 @@ def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
      
     @staticmethod
-    def get_high_ind_surf_sk(pr,element='Ni',Xstr='fcc',lattice_constant=3.526,terraceOrientation=[1,1,1],\
-              stepOrientation=[1,1,0],\
-              kinkOrientation=[1,0,1],stepDownVector=[1,1,0],\
-              lengthStep=3,lengthTerrace=3, lengthKink=1):
-                
-                
-
-            basis=pr.create_structure(element=element, bravais_basis=Xstr, lattice_constant=lattice_constant)
-            sym = basis.get_symmetry()
-            eqvdirS = np.unique(np.matmul(sym.rotations[:] , (np.array(stepOrientation))),axis=0)
-            eqvdirK = np.unique(np.matmul(sym.rotations[:] , (np.array(kinkOrientation))),axis=0)
-            eqvdirS_ind = np.where(np.dot(np.squeeze(eqvdirS),terraceOrientation)==0)[0]
-            eqvdirK_ind = np.where(np.dot(np.squeeze(eqvdirK),terraceOrientation)==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(terraceOrientation)
-            fin_kinkOrientation = eqvdirK[eqvdirK_ind[0]]
-            fin_stepOrientation = eqvdirS[temp]
-            vec1 = (np.asanyarray(fin_stepOrientation).dot(lengthStep)) + (np.asanyarray(fin_kinkOrientation).dot(lengthKink))
-            vec2 = (np.asanyarray(fin_kinkOrientation).dot(lengthTerrace)) + stepDownVector
-            highIndexSurface = np.cross(np.asanyarray(vec1),np.asanyarray(vec2))
-            highIndexSurface = np.array(highIndexSurface/np.gcd.reduce(highIndexSurface),dtype=int)
-
-         return fin_kinkOrientation, fin_stepOrientation, highIndexSurface
+    def get_high_ind_surf_sk(pr, element='Ni', Xstr='fcc', lattice_constant=3.526, terraceOrientation=[1,1,1],
+                             stepOrientation=[1,1,0], kinkOrientation=[1,0,1],stepDownVector=[1,1,0], lengthStep=3,
+                             lengthTerrace=3, lengthKink=1):
+        basis=pr.create_structure(element=element, bravais_basis=Xstr, lattice_constant=lattice_constant)
+        sym = basis.get_symmetry()
+        eqvdirS = np.unique(np.matmul(sym.rotations[:] , (np.array(stepOrientation))),axis=0)
+        eqvdirK = np.unique(np.matmul(sym.rotations[:] , (np.array(kinkOrientation))),axis=0)
+        eqvdirS_ind = np.where(np.dot(np.squeeze(eqvdirS),terraceOrientation)==0)[0]
+        eqvdirK_ind = np.where(np.dot(np.squeeze(eqvdirK),terraceOrientation)==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(terraceOrientation)
+        fin_kinkOrientation = eqvdirK[eqvdirK_ind[0]]
+        fin_stepOrientation = eqvdirS[temp]
+        vec1 = (np.asanyarray(fin_stepOrientation).dot(lengthStep)) + (np.asanyarray(fin_kinkOrientation).dot(lengthKink))
+        vec2 = (np.asanyarray(fin_kinkOrientation).dot(lengthTerrace)) + stepDownVector
+        highIndexSurface = np.cross(np.asanyarray(vec1),np.asanyarray(vec2))
+        highIndexSurface = np.array(highIndexSurface/np.gcd.reduce(highIndexSurface), dtype=int)
+        return fin_kinkOrientation, fin_stepOrientation, highIndexSurface

From 1787dba56666e263408ea2dd367c3b1db668f1ae Mon Sep 17 00:00:00 2001
From: sudarsan1989 <surendralal@mpie.de>
Date: Thu, 28 Oct 2021 16:54:33 +0200
Subject: [PATCH 04/14] Use crystal function from the class rather than a
 project

---
 pyiron_atomistics/atomistics/structure/factory.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index d421d47be..01c16549b 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -407,12 +407,11 @@ def from_pymatgen(pymatgen_obj):
     @wraps(ovito_to_pyiron)
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
-     
-    @staticmethod
-    def get_high_ind_surf_sk(pr, element='Ni', Xstr='fcc', lattice_constant=3.526, terraceOrientation=[1,1,1],
+
+    def get_high_ind_surf_sk(self, element='Ni', Xstr='fcc', lattice_constant=3.526, terraceOrientation=[1,1,1],
                              stepOrientation=[1,1,0], kinkOrientation=[1,0,1],stepDownVector=[1,1,0], lengthStep=3,
                              lengthTerrace=3, lengthKink=1):
-        basis=pr.create_structure(element=element, bravais_basis=Xstr, lattice_constant=lattice_constant)
+        basis = self.crystal(element=element, bravais_basis=Xstr, lattice_constant=lattice_constant)
         sym = basis.get_symmetry()
         eqvdirS = np.unique(np.matmul(sym.rotations[:] , (np.array(stepOrientation))),axis=0)
         eqvdirK = np.unique(np.matmul(sym.rotations[:] , (np.array(kinkOrientation))),axis=0)

From 1a786df1edef4f37b1cf5e3c6641b1bfed36e24a Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Fri, 29 Oct 2021 13:17:01 +0200
Subject: [PATCH 05/14] Made changes according to the input from the reviewers.

---
 .../atomistics/structure/factory.py           | 134 ++++++++++++------
 1 file changed, 91 insertions(+), 43 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index 01c16549b..66f2ff433 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/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
@@ -30,7 +30,8 @@
 from pyiron_atomistics.atomistics.structure.factories.aimsgb import AimsgbFactory
 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 Settings, PyironFactory, deprecate
 import types
@@ -49,6 +50,7 @@
 
 s = Settings()
 
+
 class StructureFactory(PyironFactory):
     def __init__(self):
         self._ase = AseFactory()
@@ -71,24 +73,29 @@ def aimsgb(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(
@@ -132,7 +139,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.
@@ -220,7 +227,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)
 
@@ -369,6 +376,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")
@@ -391,6 +399,7 @@ def aimsgb_build(
             delete_layer=delete_layer,
             add_if_dist=add_if_dist
         )
+
     aimsgb_build.__doc__ = AimsgbFactory.build.__doc__
 
     @staticmethod
@@ -408,31 +417,70 @@ def from_pymatgen(pymatgen_obj):
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
 
-    def get_high_ind_surf_sk(self, element='Ni', Xstr='fcc', lattice_constant=3.526, terraceOrientation=[1,1,1],
-                             stepOrientation=[1,1,0], kinkOrientation=[1,0,1],stepDownVector=[1,1,0], lengthStep=3,
-                             lengthTerrace=3, lengthKink=1):
-        basis = self.crystal(element=element, bravais_basis=Xstr, lattice_constant=lattice_constant)
-        sym = basis.get_symmetry()
-        eqvdirS = np.unique(np.matmul(sym.rotations[:] , (np.array(stepOrientation))),axis=0)
-        eqvdirK = np.unique(np.matmul(sym.rotations[:] , (np.array(kinkOrientation))),axis=0)
-        eqvdirS_ind = np.where(np.dot(np.squeeze(eqvdirS),terraceOrientation)==0)[0]
-        eqvdirK_ind = np.where(np.dot(np.squeeze(eqvdirK),terraceOrientation)==0)[0]
 
-        if len(eqvdirS_ind) == 0:
+    def get_high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
+                               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 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].
+
+        Additionally also returns a bottom slab with the calculated high index surface
+
+        [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
+            layers (int): Number of layers of the high_index_surface eg., 60
+            vacuum (float): Thickness of vacuum on the top of the slab
+
+        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
+            slab: pyiron_atomistics.atomistics.structure.atoms.Atoms instance Required surface
+        """
+        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:
+        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)
+        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
+        slab.set_initial_magnetic_moments(np.repeat(1.0, len(slab)))
+        return slab, high_index_surface, fin_kink_orientation, fin_step_orientation
+
 
-        temp = (np.cross(np.squeeze(eqvdirK[eqvdirK_ind[0]]), np.squeeze(eqvdirS))).tolist().index(terraceOrientation)
-        fin_kinkOrientation = eqvdirK[eqvdirK_ind[0]]
-        fin_stepOrientation = eqvdirS[temp]
-        vec1 = (np.asanyarray(fin_stepOrientation).dot(lengthStep)) + (np.asanyarray(fin_kinkOrientation).dot(lengthKink))
-        vec2 = (np.asanyarray(fin_kinkOrientation).dot(lengthTerrace)) + stepDownVector
-        highIndexSurface = np.cross(np.asanyarray(vec1),np.asanyarray(vec2))
-        highIndexSurface = np.array(highIndexSurface/np.gcd.reduce(highIndexSurface), dtype=int)
-        return fin_kinkOrientation, fin_stepOrientation, highIndexSurface

From becccc3974bfb49074bc1da1b06ffdc392ae9e89 Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Fri, 29 Oct 2021 13:21:53 +0200
Subject: [PATCH 06/14] Test for get high index surface

---
 .../structure/test_high_index_surface.py      | 22 +++++++++++++++++++
 1 file changed, 22 insertions(+)
 create mode 100644 tests/atomistics/structure/test_high_index_surface.py

diff --git a/tests/atomistics/structure/test_high_index_surface.py b/tests/atomistics/structure/test_high_index_surface.py
new file mode 100644
index 000000000..aa9d0ea26
--- /dev/null
+++ b/tests/atomistics/structure/test_high_index_surface.py
@@ -0,0 +1,22 @@
+import unittest
+from pyiron_atomistics.atomistics.structure.factory import StructureFactory
+
+
+class TestHighIndexSurface(unittest.TestCase):
+    def test_high_index_surface(self):
+        k, s, h, slab = StructureFactory().get_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(h, [-13, -7, -9])
+        self.assertEqual(k, [-1, 0, 1])
+        self.assertEqual(s, [0, -1, 1])
+
+
+if __name__ == '__main__':
+    unittest.main()

From d09c41f05a9f986f185b5d3f3643104cfda721fd Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Fri, 29 Oct 2021 14:20:22 +0200
Subject: [PATCH 07/14] Changes to the order of instances returned

---
 pyiron_atomistics/atomistics/structure/factory.py | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index 66f2ff433..0c81c408b 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -417,15 +417,14 @@ def from_pymatgen(pymatgen_obj):
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
 
-
     def get_high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
-                               terrace_orientation=None,step_orientation=None, kink_orientation=None,
+                               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 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].
+        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].
 
         Additionally also returns a bottom slab with the calculated high index surface
 
@@ -482,5 +481,3 @@ def get_high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_
         slab.set_pbc = True
         slab.set_initial_magnetic_moments(np.repeat(1.0, len(slab)))
         return slab, high_index_surface, fin_kink_orientation, fin_step_orientation
-
-

From 2ad076445cd890522f88aefb168eaa4bb0f57eb2 Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Wed, 3 Nov 2021 12:49:57 +0100
Subject: [PATCH 08/14] Changes to the unittest

---
 tests/atomistics/structure/test_high_index_surface.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/tests/atomistics/structure/test_high_index_surface.py b/tests/atomistics/structure/test_high_index_surface.py
index aa9d0ea26..50a37d969 100644
--- a/tests/atomistics/structure/test_high_index_surface.py
+++ b/tests/atomistics/structure/test_high_index_surface.py
@@ -4,7 +4,7 @@
 
 class TestHighIndexSurface(unittest.TestCase):
     def test_high_index_surface(self):
-        k, s, h, slab = StructureFactory().get_high_index_surface(element='Ni', crystal_structure='fcc',
+        slab, h, s, k = StructureFactory().get_high_index_surface(element='Ni', crystal_structure='fcc',
                                                                   lattice_constant=3.526,
                                                                   terrace_orientation=[1, 1, 1],
                                                                   step_orientation=[1, 1, 0],
@@ -13,9 +13,10 @@ def test_high_index_surface(self):
                                                                   length_terrace=3,
                                                                   length_kink=1, layers=60,
                                                                   vacuum=10)
-        self.assertEqual(h, [-13, -7, -9])
-        self.assertEqual(k, [-1, 0, 1])
-        self.assertEqual(s, [0, -1, 1])
+        self.assertEqual(len(h), 3)
+        self.assertEqual(h[0], -9)
+        self.assertEqual(len(k), 3)
+        self.assertEqual(len(s), 3)
 
 
 if __name__ == '__main__':

From 66d9d8295b3d805f35d85dadbb8b5c00996ffeff Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Mon, 8 Nov 2021 06:05:12 +0100
Subject: [PATCH 09/14] Changes to the unittest and the name of the function.
 Unittest now has a test for length of the slab.

---
 .../atomistics/structure/factory.py           |  8 ++++----
 .../structure/test_high_index_surface.py      | 19 ++++++++++---------
 2 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index 0c81c408b..f93e9b0a1 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -417,10 +417,10 @@ def from_pymatgen(pymatgen_obj):
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
 
-    def get_high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
-                               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):
+    def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
+                           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 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
diff --git a/tests/atomistics/structure/test_high_index_surface.py b/tests/atomistics/structure/test_high_index_surface.py
index 50a37d969..41a12764d 100644
--- a/tests/atomistics/structure/test_high_index_surface.py
+++ b/tests/atomistics/structure/test_high_index_surface.py
@@ -4,19 +4,20 @@
 
 class TestHighIndexSurface(unittest.TestCase):
     def test_high_index_surface(self):
-        slab, h, s, k = StructureFactory().get_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)
+        slab, h, s, k = 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(h), 3)
         self.assertEqual(h[0], -9)
         self.assertEqual(len(k), 3)
         self.assertEqual(len(s), 3)
+        self.assertEqual(len(slab), 60)
 
 
 if __name__ == '__main__':

From 50b2549d85250528d2db90309d698ef4dbe30135 Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Tue, 9 Nov 2021 08:36:30 +0100
Subject: [PATCH 10/14] Changes to the unittest. Added valueerror assertions.
 Refactored the function into two functions high_index_surface which returns
 only structure and high_index_surface_info which returns the Miller indices
 of the surface place, step and kink directions.

---
 .../atomistics/structure/factory.py           | 59 +++++++++++++++----
 .../structure/test_high_index_surface.py      | 40 +++++++++----
 2 files changed, 78 insertions(+), 21 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index f93e9b0a1..de59b802d 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -417,17 +417,14 @@ def from_pymatgen(pymatgen_obj):
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
 
-    def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
-                           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):
+    def high_index_surface_info(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
+                                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].
 
-        Additionally also returns a bottom slab with the calculated high index surface
-
         [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.
 
@@ -442,15 +439,18 @@ def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_cons
             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
-            layers (int): Number of layers of the high_index_surface eg., 60
-            vacuum (float): Thickness of vacuum on the top of the slab
+
 
         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
-            slab: pyiron_atomistics.atomistics.structure.atoms.Atoms instance Required surface
         """
+        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)
@@ -473,11 +473,48 @@ def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_cons
         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='Ni', crystal_structure='fcc', lattice_constant=3.526,
+                           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, \
+            fin_kink_orientation, \
+            fin_step_orientation = self.high_index_surface_info(element='Ni',
+                                                                crystal_structure='fcc', lattice_constant=3.526,
+                                                                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
-        slab.set_initial_magnetic_moments(np.repeat(1.0, len(slab)))
-        return slab, high_index_surface, fin_kink_orientation, fin_step_orientation
+        return slab
diff --git a/tests/atomistics/structure/test_high_index_surface.py b/tests/atomistics/structure/test_high_index_surface.py
index 41a12764d..73728d996 100644
--- a/tests/atomistics/structure/test_high_index_surface.py
+++ b/tests/atomistics/structure/test_high_index_surface.py
@@ -4,20 +4,40 @@
 
 class TestHighIndexSurface(unittest.TestCase):
     def test_high_index_surface(self):
-        slab, h, s, k = 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)
+        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)
-        self.assertEqual(len(slab), 60)
+        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__':

From cb7fe90081020cc3a36ac89bfe3c91d8d7533748 Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Thu, 11 Nov 2021 11:51:40 +0100
Subject: [PATCH 11/14] Changes to the unittest. Added valueerror assertions.
 Refactored the function into two functions high_index_surface which returns
 only structure and high_index_surface_info which returns the Miller indices
 of the surface place, step and kink directions.

made changes to high_index_surface function to use the given element, crystal_structure and lattice_constant.
---
 pyiron_atomistics/atomistics/structure/factory.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index de59b802d..e35ff080f 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -502,8 +502,9 @@ def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_cons
         basis = self.crystal(element=element, bravais_basis=crystal_structure, lattice_constant=lattice_constant)
         high_index_surface, \
             fin_kink_orientation, \
-            fin_step_orientation = self.high_index_surface_info(element='Ni',
-                                                                crystal_structure='fcc', lattice_constant=3.526,
+            fin_step_orientation = 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,

From a3f8f7f5ed322cf2d1d143e462a94eae7bb3ae33 Mon Sep 17 00:00:00 2001
From: skatnagallu <s.katnagallu@mpie.de>
Date: Wed, 17 Nov 2021 15:41:06 +0100
Subject: [PATCH 12/14] Changes to the unittest. Added valueerror assertions.
 Refactored the function into two functions high_index_surface which returns
 only structure and high_index_surface_info which returns the Miller indices
 of the surface place, step and kink directions.

corrected it to give periodic boundary conditions in all three directions.
---
 pyiron_atomistics/atomistics/structure/factory.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index e35ff080f..d3d46606a 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -517,5 +517,5 @@ def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_cons
         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
+        slab.set_pbc(True)
         return slab

From 3cd2981740a8b1d3edde0997ccd1b4b925873384 Mon Sep 17 00:00:00 2001
From: Sudarsan Surendralal <surendralal@mpie.de>
Date: Tue, 18 Jan 2022 17:38:43 +0100
Subject: [PATCH 13/14] Apply suggestions from code review

---
 pyiron_atomistics/atomistics/structure/factory.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index d3d46606a..1c725eeb2 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -500,9 +500,7 @@ def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_cons
             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, \
-            fin_kink_orientation, \
-            fin_step_orientation = self.high_index_surface_info(element=element,
+        high_index_surface, _, _ = self.high_index_surface_info(element=element,
                                                                 crystal_structure=crystal_structure,
                                                                 lattice_constant=lattice_constant,
                                                                 terrace_orientation=terrace_orientation,

From af372f6648d580fbd7dbfbb1f41f926a1fe31f26 Mon Sep 17 00:00:00 2001
From: sudarsan1989 <surendralal@mpie.de>
Date: Tue, 18 Jan 2022 17:40:53 +0100
Subject: [PATCH 14/14] Remove default values

---
 pyiron_atomistics/atomistics/structure/factory.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/pyiron_atomistics/atomistics/structure/factory.py b/pyiron_atomistics/atomistics/structure/factory.py
index 9053bb619..ab501a7ea 100644
--- a/pyiron_atomistics/atomistics/structure/factory.py
+++ b/pyiron_atomistics/atomistics/structure/factory.py
@@ -50,7 +50,6 @@
 __date__ = "May 1, 2020"
 
 
-
 class StructureFactory(PyironFactory):
     def __init__(self):
         self._ase = AseFactory()
@@ -422,7 +421,7 @@ def from_pymatgen(pymatgen_obj):
     def from_ovito(ovito_obj):
         return ovito_to_pyiron(ovito_obj)
 
-    def high_index_surface_info(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
+    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):
         """
@@ -481,7 +480,7 @@ def high_index_surface_info(self, element='Ni', crystal_structure='fcc', lattice
 
         return high_index_surface, fin_kink_orientation, fin_step_orientation
 
-    def high_index_surface(self, element='Ni', crystal_structure='fcc', lattice_constant=3.526,
+    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):