Create Surface Slabs in FCC Crystals¶

In this code, a slab for simulating free-surfaces are created. For given Miller indices $(h \, k \, l)$ of the plane, input configuration containing the simulation-cell and atoms is created.

Load relevant modules
Define necessary parameters
- Miller indices
- Lattice (e.g. Aluminum)
- Cut-off radius (r_cut) and area (a_cut)
Compute basis-vectos of the 2D lattice of the $(h \, k \, l)$ plane
Compute the basis-vectors of sublattice with r_cut and a_cut constraints
- Determine appropriate $\Sigma$ using a_cut
- Generate unique set of sub-lattices (HNF matrices and LLL reduction)

[1]:

import byxtal.lattice as gbl
import byxtal.integer_manipulations as iman
import byxtal.bp_basis as bpb
import numpy as np
import numpy.linalg as nla
import math as mt
import gbpy.util_funcs_create_byxtal as uf
import gbpy.generate_hkl_indices as ghi

      
    

(h k l) : Miller indices defined in primitive lattice
l1 : Lattice object
lat_par : Lattice parameter
r_cut : Cutoff radius
a_cut : Cutoff Area
zCut : The lenght of the surface-slab is equal to 2(zCut)

[2]:

h, k, l = -2, -2 , -2
hkl_p = np.array([h,k,l])

l1 = gbl.Lattice('Al')

lat_par = l1.lat_params['a']
r_cut = lat_par*4
a_cut = r_cut**2
zCut = 25*lat_par

l_p_po : $\Lambda_{p}^{po}$
l_po_p : $\Lambda_{po}^{p}$
l_bpb_p : $\Lambda_{\mathrm{2D}}^{p}$
l_bpb_po : $\Lambda_{\mathrm{2D}}^{po}$

[65]:

# import numpy as np
# from math import gcd
# import os

# def bp_basis(miller_ind):
#     """
#     The function computes the primitve basis of the plane if the
#     boundary plane indices are specified

#     Parameters
#     ---------------
#     miller_ind: numpy array
#         Miller indices of the plane (h k l)

#     Returns
#     -----------
#     l_pl_g1: numpy array
#         The primitive basis of the plane in 'g1' reference frame
#     """
#     # If *miller_inds* are not integers or if the gcd != 1
#     # miller_ind = int_man.int_finder(miller_ind)
#     if (np.ndim(miller_ind) == 2):
#         Sz = np.shape(miller_ind)
#         if ((Sz[0] == 1) or (Sz[1] == 1)):
#             miller_ind = miller_ind.flatten()
#         else:
#             raise Exception("Wrong Input Type.")
#     h = miller_ind[0]
#     k = miller_ind[1]
#     l = miller_ind[2]

#     if h == 0 and k == 0 and l == 0:
#         raise Exception('hkl indices cannot all be zero')
#     else:
#         if h != 0 and k != 0 and l != 0:
#             gc_f1_p = gcd(k, l)
#             bv1_g1 = np.array([[0], [-l / gc_f1_p], [k / gc_f1_p]])
#             bv2_g1 = bpb.compute_basis_vec([h, k, l])
#             bv2_g1 = bv2_g1.reshape(np.shape(bv2_g1)[0],1)
#         else:
#                 if h == 0:
#                     if k == 0:
#                         bv1_g1 = np.array([[1], [0], [0]])
#                         bv2_g1 = np.array([[0], [1], [0]])
#                     elif l == 0:
#                         bv1_g1 = np.array([[0], [0], [1]])
#                         bv2_g1 = np.array([[1], [0], [0]])
#                     else:
#                         gc_f1_p = gcd(k, l)
#                         bv1_g1 = np.array([[0], [-l / gc_f1_p],
#                                            [k / gc_f1_p]])
#                         bv2_g1 = np.array([[1], [-l / gc_f1_p],
#                                            [k / gc_f1_p]])
#                 else:
#                     if k == 0:
#                         if l == 0:
#                             bv1_g1 = np.array([[0], [1], [0]])
#                             bv2_g1 = np.array([[0], [0], [1]])
#                         else:
#                             gc_f1_p = gcd(h, l)
#                             bv1_g1 = np.array([[-l / gc_f1_p], [0],
#                                                [h / gc_f1_p]])
#                             bv2_g1 = np.array([[-l / gc_f1_p], [1],
#                                                [h / gc_f1_p]])
#                     else:
#                         if l == 0:
#                             gc_f1_p = gcd(h, k)
#                             bv1_g1 = np.array([[-k / gc_f1_p],
#                                                [h / gc_f1_p], [0]])
#                             bv2_g1 = np.array([[-k / gc_f1_p],
#                                                [h / gc_f1_p], [1]])

#     #  The reduced basis vectors for the plane
#     Tmat = np.array(np.column_stack([bv1_g1, bv2_g1]), dtype='int64')
#     l_pl_g1 = lll_reduction(Tmat)
#     return l_pl_g1

# def lll_reduction(int_mat):
#     """
#     """

#     exec_str = '/compute_lll.py'
#     inp_args = {}
#     inp_args['mat'] = int_mat
#     lll_int_mat = call_sage_math(exec_str, inp_args)
#     return lll_int_mat

# def call_sage_math(exec_str, inp_args):
#     print(exec_str)
#     print(inp_args)
#     byxtal_dir = os.path.dirname((inspect.getfile(byxtal)))
#     exec_str1 = byxtal_dir+exec_str
#     run_lst = []
#     run_lst.append(exec_str1)

#     A = inp_args['mat']
#     Sz = np.shape(A)
#     run_lst.append(str(Sz[0]))
#     run_lst.append(str(Sz[1]))

#     if len(inp_args.keys()) == 2:
#         sig_num = inp_args['sig_num']
#         str1 = str(sig_num)
#         run_lst.append(str1)

#     for i1 in range(Sz[0]):
#         for j1 in range(Sz[1]):
#             str1 = str(A[i1, j1])
#             run_lst.append(str1)

#     result = subprocess.run(run_lst, stdout=subprocess.PIPE)
#     print('result', result)
#     print('result.stdout', result.stdout)
#     str_out = (result.stdout).split()
#     print('str_out', str_out)

#     sz1 = len(str_out)
#     M_out = np.zeros((Sz[0],Sz[1]), dtype='int64')
#     print(M_out)
#     ct1 = 0
#     for i1 in range(Sz[0]):
#         for j1 in range(Sz[1]):
#             print('i1', i1)
#             print('j1', j1)
#             print('int', str_out)
#             M_out[i1, j1] = int(str_out[ct1])
#             ct1 = ct1 + 1

#     return M_out

      
    

[66]:

# import inspect
# import byxtal
# import subprocess
# l_p_po = l1.l_p_po
# l_po_p = nla.inv(l_p_po)

# ## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)
# l_bpb_p = bp_basis(hkl_p)

      
    

/compute_lll.py
{'mat': array([[ 0,  1],
       [ 1,  0],
       [-1, -1]])}
result CompletedProcess(args=['/home/leila/Leila_sndhard/codes/byxtal/byxtal/byxtal/compute_lll.py', '3', '2', '0', '1', '1', '0', '-1', '-1'], returncode=126, stdout=b'')
result.stdout b''
str_out []
[[0 0]
 [0 0]
 [0 0]]
i1 0
j1 0
int []

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-66-4bb23afa5176> in <module>
      6
      7 ## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)
----> 8 l_bpb_p = bp_basis(hkl_p)

<ipython-input-65-0ac95a520c52> in bp_basis(miller_ind)
     73     #  The reduced basis vectors for the plane
     74     Tmat = np.array(np.column_stack([bv1_g1, bv2_g1]), dtype='int64')
---> 75     l_pl_g1 = lll_reduction(Tmat)
     76     return l_pl_g1
     77

<ipython-input-65-0ac95a520c52> in lll_reduction(int_mat)
     83     inp_args = {}
     84     inp_args['mat'] = int_mat
---> 85     lll_int_mat = call_sage_math(exec_str, inp_args)
     86     return lll_int_mat
     87

<ipython-input-65-0ac95a520c52> in call_sage_math(exec_str, inp_args)
    124             print('j1', j1)
    125             print('int', str_out)
--> 126             M_out[i1, j1] = int(str_out[ct1])
    127             ct1 = ct1 + 1
    128

IndexError: list index out of range

[16]:

l_p_po = l1.l_p_po
l_po_p = nla.inv(l_p_po)

## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)
l_bpb_p = bpb.bp_basis(hkl_p)
l_bpb_p = l_bpb_p.astype(int)
## l_bpb_p: Primitive Basis vectors of the boundary-plane (in po reference frame)
l_bpb_po = l_p_po.dot(l_bpb_p)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-16-4338a158805e> in <module>
      3
      4 ## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)
----> 5 l_bpb_p = bpb.bp_basis(hkl_p)
      6 l_bpb_p = l_bpb_p.astype(int)
      7 ## l_bpb_p: Primitive Basis vectors of the boundary-plane (in po reference frame)

~/Leila_sndhard/codes/byxtal/byxtal/byxtal/bp_basis.py in bp_basis(miller_ind)
    205     #  The reduced basis vectors for the plane
    206     Tmat = np.array(np.column_stack([bv1_g1, bv2_g1]), dtype='int64')
--> 207     l_pl_g1 = lt.lll_reduction(Tmat)
    208     return l_pl_g1
    209

~/Leila_sndhard/codes/byxtal/byxtal/byxtal/lll_tools.py in lll_reduction(int_mat)
    107     inp_args = {}
    108     inp_args['mat'] = int_mat
--> 109     lll_int_mat = rpl.call_sage_math(exec_str, inp_args)
    110     return lll_int_mat
    111

~/Leila_sndhard/codes/byxtal/byxtal/byxtal/reduce_po_lat.py in call_sage_math(exec_str, inp_args)
     36     for i1 in range(Sz[0]):
     37         for j1 in range(Sz[1]):
---> 38             M_out[i1, j1] = int(str_out[ct1])
     39             ct1 = ct1 + 1
     40

IndexError: list index out of range

sig_num: $\Sigma$ number
l_bpbSig_p : $\Lambda_{\Sigma-\mathrm{2D}}^{p}$
l_bpbSig_po : $\Lambda_{\Sigma-\mathrm{2D}}^{po}$

[12]:

## area_bpl: Area of the 2D-primitive-unit-cell
area_bpl = nla.norm((np.cross(l_bpb_po[:,0], l_bpb_po[:,1])))
sig_num = np.ceil(a_cut/area_bpl)

ind2 = np.array([], dtype='int64');
while (np.size(ind2) == 0):
    # Generate 2D hermite normal forms for sig_num (hnf_mats)
    hnf_mats = ghi.sig_hnf_mats(sig_num)
    # Compute the properties of the sub-lattices
    l_sig_p_mats, l_sig_po_mats = ghi.compute_hnf_props(hnf_mats, l_bpb_p, l_p_po, 1e-2)
    # Get the index for the sub-lattice that has the minimum cost
    ind2 = ghi.ind_min_cost(l_sig_po_mats, r_cut)
    sig_num = sig_num + 1

## l_bpbSig_p: Basis vectors of the sub-lattice of the boundary-plane (in p reference frame)
l_bpbSig_p = l_sig_p_mats[ind2];
l_bpbSig_p = l_bpbSig_p.astype(int)
l_bpbSig_po = l_p_po.dot(l_bpbSig_p)

      
    

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-12-94b79fcf1741> in <module>
      1 ## area_bpl: Area of the 2D-primitive-unit-cell
----> 2 area_bpl = nla.norm((np.cross(l_bpb_po[:,0], l_bpb_po[:,1])))
      3 sig_num = np.ceil(a_cut/area_bpl)
      4
      5 ind2 = np.array([], dtype='int64');

NameError: name 'l_bpb_po' is not defined

[ ]:

l_sig_p_mats

l_po_go: $\Lambda_{PO}^{GO}$ - Orientation of crystal with respect to global reference frame $\mathrm{GO}$
- go: Orthogonal global reference frame (for LAMMPS)
- po: Orhtogonal reference frame of the lattice
- function compute_orientation: computes the orientation of the crystal such that the basis vectors in $\Lambda_{\Sigma-\mathrm{2D}}^{po}$ are along the x-axis and in the xy-plane.
l2D_bpbSig_go: $\Lambda_{\Sigma-\mathrm{2D}}^{go}$
l2D_bpb_go: $\Lambda_{\Sigma-\mathrm{2D}}^{go}$

If computed accurately, the $z$ -components of the basis-vectors in the $\mathrm{GO}$ frame, $\Lambda_{\mathrm{2D}}^{go}$ and $\Lambda_{\Sigma-\mathrm{2D}}^{go}$ , are zero.

twoD_mat: The $x$ , $y$ components of the 2D basis vectors of $\Lambda_{\mathrm{2D}}^{go}$
twoDSig_mat: The $x$ , $y$ components of the 2D basis vectors of $\Lambda_{\Sigma-\mathrm{2D}}^{go}$
r_cut1: The cut-off radius for replicating the lattice basis. Given two vectors that form the basis of the interface plane, $\Lambda_{\Sigma-\mathrm{2D}}^{go}$ , the function compute_rCut determines the maximum of the norms of the two vectors.
twoD_pts: The lattice points in the 2D periodic box for the 3D slab.

[ ]:

l_po_go = uf.compute_orientation(l_bpbSig_po)
l2D_bpbSig_go = l_po_go.dot(l_bpbSig_po)
l2D_bpb_go = l_po_go.dot(l_bpb_po)
twoD_mat = l2D_bpb_go[:2,:]
twoDSig_mat = l2D_bpbSig_go[:2,:]
r_cut1 = uf.compute_rCut(l_bpbSig_po)

## Create 2D periodic box
twoD_pts = uf.replicate_pts(twoD_mat, r_cut1)
twoD_pts = uf.remove_periodic_overlaps(twoD_pts, twoDSig_mat)

      
    

avec: Basis vector $\vec{a}$ of the 3D lattice defined using $\mathcal{B}_{P2}$ .
bvec: Basis vector $\vec{b}$ of the 3D lattice defined using $\mathcal{B}_{P2}$ .
l_p2_p1: Equivalent basis for the 3D lattice, $\Lambda_{P2}^{P1}$ .
- The function find_int_solns computes the third basis-vector, $\vec{c}$ , such that $\mathcal{B}_{P1} \equiv \mathcal{B}_{P2}$

[ ]:

avec = l_bpb_p[:,0]
bvec = l_bpb_p[:,1]
l_p2_p1 = uf.find_int_solns(avec, bvec)

[ ]:

l2D_bpbSig_go

l_p2_po1: $\Lambda_{p2}^{po1}$
l_p2_go: $\Lambda_{p2}^{go}$
tz_vec: The components of $\vec{c}$ of $\mathcal{B}_{P2}$ in $go$ reference frame.
threeD_pts: The 3D slab of atoms for simulating free surfaces.
- This is obtained by replicating twoD_pts along the tz_vec in both + $z$ and - $z$ directions.

[ ]:

l_p2_po1 = l_p_po.dot(l_p2_p1)
l_p2_go = l_po_go.dot(l_p2_po1)
tz_vec = np.array(l_p2_go[:,2], dtype='double')
tz_vec = np.reshape(tz_vec, (3,))

################################################################################
## Translate 2D points in the Z-direction with zCut
num_rep = np.abs(int(np.ceil(zCut/tz_vec[2])))
num_2d = np.shape(twoD_pts)[0]
num_3d_pts = int((2*num_rep+1)*num_2d)
threeD_pts = np.zeros((num_3d_pts,3));

twoD_pts1 = np.hstack((twoD_pts, np.zeros((num_2d,1))));

for ct1 in np.arange(-num_rep, num_rep+1):
    ct2 = ct1 + num_rep
    ind_st = (ct2)*num_2d
    ind_stop = ind_st + num_2d
    trans_vec = tz_vec*ct1
    threeD_pts[ind_st:ind_stop, :] = twoD_pts1 + np.tile(trans_vec, (num_2d,1))

      
    

The simulation box vectors are defined using Ovito’s convenction. 1. sim_avec: The components of $\vec{a}$ in the global reference, $go$ , of LAMMPS. $\vec{a}$ is along $\hat{e}_x$ . 2. sim_bvec: The components of $\vec{b}$ in the global reference, $go$ , of LAMMPS. $\vec{b}$ lies in the $x-y$ plane. 3. sim_cvec: The components of $\vec{c}$ in the global reference, $go$ , of LAMMPS. 4. sim_orig: The origin of the simulation box.

[ ]:

### Simulation Cell Box
### Following Ovito's convention
l_bpbSig_po_arr = np.array(l2D_bpbSig_go, dtype='double')
sim_cell = np.zeros((3,4))
sim_avec = l_bpbSig_po_arr[:,0]
sim_bvec = l_bpbSig_po_arr[:,1]

### Change this with inter-planar spacing
sim_cvec = np.array([0,0,2*zCut]); # sim_cvec = np.array([0,0,zCut]);
sim_orig = np.array([0,0,-zCut]); # sim_orig = np.array([0,0,0]);

sim_cell[:,0] = sim_avec
sim_cell[:,1] = sim_bvec
sim_cell[:,2] = sim_cvec
sim_cell[:,3] = sim_orig

box_vecs = sim_cell[:,0:3]
threeD_pts1 = uf.wrap_cc(sim_cell, threeD_pts)
# tpts1 = np.dot(np.linalg.inv(box_vecs), threeD_pts.transpose()).transpose();

      
    

Plot of 3D lattice points along with the simulation box.

[ ]:

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import plotting_routines as plr
%matplotlib inline

fig1 = plt.figure()
plr.plot_3d_pts_box(fig1, threeD_pts1, sim_cell[:,0:3], sim_orig)
plt.show()

      
    