{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Create Surface Slabs in FCC Crystals\n", "\n", "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.\n", "\n", "1. Load relevant modules\n", "2. Define necessary parameters\n", " + Miller indices\n", " + Lattice (e.g. Aluminum)\n", " + Cut-off radius (r_cut) and area (a_cut)\n", "3. Compute basis-vectos of the 2D lattice of the $(h \\, k \\, l)$ plane\n", "4. Compute the basis-vectors of sublattice with r_cut and a_cut constraints\n", " + Determine appropriate $\\Sigma$ using a_cut\n", " + Generate unique set of sub-lattices (HNF matrices and LLL reduction) " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import byxtal.lattice as gbl\n", "import byxtal.integer_manipulations as iman\n", "import byxtal.bp_basis as bpb\n", "import numpy as np\n", "import numpy.linalg as nla\n", "import math as mt\n", "import gbpy.util_funcs_create_byxtal as uf\n", "import gbpy.generate_hkl_indices as ghi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. `(h k l)` : Miller indices defined in primitive lattice\n", "2. `l1` : Lattice object\n", "3. `lat_par` : Lattice parameter\n", "4. `r_cut` : Cutoff radius\n", "5. `a_cut` : Cutoff Area\n", "6. `zCut` : The lenght of the surface-slab is equal to 2(zCut)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "h, k, l = -2, -2 , -2\n", "hkl_p = np.array([h,k,l])\n", "\n", "l1 = gbl.Lattice('Al')\n", "\n", "lat_par = l1.lat_params['a']\n", "r_cut = lat_par*4\n", "a_cut = r_cut**2\n", "zCut = 25*lat_par" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. `l_p_po` : $\\Lambda_{p}^{po}$\n", "2. `l_po_p` : $\\Lambda_{po}^{p}$\n", "3. `l_bpb_p` : $\\Lambda_{\\mathrm{2D}}^{p}$\n", "4. `l_bpb_po` : $\\Lambda_{\\mathrm{2D}}^{po}$" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "# import numpy as np\n", "# from math import gcd\n", "# import os\n", "\n", "# def bp_basis(miller_ind):\n", "# \"\"\"\n", "# The function computes the primitve basis of the plane if the\n", "# boundary plane indices are specified\n", "\n", "# Parameters\n", "# ---------------\n", "# miller_ind: numpy array\n", "# Miller indices of the plane (h k l)\n", "\n", "# Returns\n", "# -----------\n", "# l_pl_g1: numpy array\n", "# The primitive basis of the plane in 'g1' reference frame\n", "# \"\"\"\n", "# # If *miller_inds* are not integers or if the gcd != 1\n", "# # miller_ind = int_man.int_finder(miller_ind)\n", "# if (np.ndim(miller_ind) == 2):\n", "# Sz = np.shape(miller_ind)\n", "# if ((Sz[0] == 1) or (Sz[1] == 1)):\n", "# miller_ind = miller_ind.flatten()\n", "# else:\n", "# raise Exception(\"Wrong Input Type.\")\n", "# h = miller_ind[0]\n", "# k = miller_ind[1]\n", "# l = miller_ind[2]\n", "\n", "# if h == 0 and k == 0 and l == 0:\n", "# raise Exception('hkl indices cannot all be zero')\n", "# else:\n", "# if h != 0 and k != 0 and l != 0:\n", "# gc_f1_p = gcd(k, l)\n", "# bv1_g1 = np.array([[0], [-l / gc_f1_p], [k / gc_f1_p]])\n", "# bv2_g1 = bpb.compute_basis_vec([h, k, l])\n", "# bv2_g1 = bv2_g1.reshape(np.shape(bv2_g1)[0],1)\n", "# else:\n", "# if h == 0:\n", "# if k == 0:\n", "# bv1_g1 = np.array([[1], [0], [0]])\n", "# bv2_g1 = np.array([[0], [1], [0]])\n", "# elif l == 0:\n", "# bv1_g1 = np.array([[0], [0], [1]])\n", "# bv2_g1 = np.array([[1], [0], [0]])\n", "# else:\n", "# gc_f1_p = gcd(k, l)\n", "# bv1_g1 = np.array([[0], [-l / gc_f1_p],\n", "# [k / gc_f1_p]])\n", "# bv2_g1 = np.array([[1], [-l / gc_f1_p],\n", "# [k / gc_f1_p]])\n", "# else:\n", "# if k == 0:\n", "# if l == 0:\n", "# bv1_g1 = np.array([[0], [1], [0]])\n", "# bv2_g1 = np.array([[0], [0], [1]])\n", "# else:\n", "# gc_f1_p = gcd(h, l)\n", "# bv1_g1 = np.array([[-l / gc_f1_p], [0],\n", "# [h / gc_f1_p]])\n", "# bv2_g1 = np.array([[-l / gc_f1_p], [1],\n", "# [h / gc_f1_p]])\n", "# else:\n", "# if l == 0:\n", "# gc_f1_p = gcd(h, k)\n", "# bv1_g1 = np.array([[-k / gc_f1_p],\n", "# [h / gc_f1_p], [0]])\n", "# bv2_g1 = np.array([[-k / gc_f1_p],\n", "# [h / gc_f1_p], [1]])\n", "\n", "# # The reduced basis vectors for the plane\n", "# Tmat = np.array(np.column_stack([bv1_g1, bv2_g1]), dtype='int64')\n", "# l_pl_g1 = lll_reduction(Tmat)\n", "# return l_pl_g1\n", "\n", "# def lll_reduction(int_mat):\n", "# \"\"\"\n", "# \"\"\"\n", "\n", "# exec_str = '/compute_lll.py'\n", "# inp_args = {}\n", "# inp_args['mat'] = int_mat\n", "# lll_int_mat = call_sage_math(exec_str, inp_args)\n", "# return lll_int_mat\n", "\n", "# def call_sage_math(exec_str, inp_args):\n", "# print(exec_str)\n", "# print(inp_args)\n", "# byxtal_dir = os.path.dirname((inspect.getfile(byxtal)))\n", "# exec_str1 = byxtal_dir+exec_str\n", "# run_lst = []\n", "# run_lst.append(exec_str1)\n", "\n", "# A = inp_args['mat']\n", "# Sz = np.shape(A)\n", "# run_lst.append(str(Sz[0]))\n", "# run_lst.append(str(Sz[1]))\n", "\n", "# if len(inp_args.keys()) == 2:\n", "# sig_num = inp_args['sig_num']\n", "# str1 = str(sig_num)\n", "# run_lst.append(str1)\n", "\n", "# for i1 in range(Sz[0]):\n", "# for j1 in range(Sz[1]):\n", "# str1 = str(A[i1, j1])\n", "# run_lst.append(str1)\n", "\n", "# result = subprocess.run(run_lst, stdout=subprocess.PIPE)\n", "# print('result', result)\n", "# print('result.stdout', result.stdout)\n", "# str_out = (result.stdout).split()\n", "# print('str_out', str_out)\n", "\n", "# sz1 = len(str_out)\n", "# M_out = np.zeros((Sz[0],Sz[1]), dtype='int64')\n", "# print(M_out)\n", "# ct1 = 0\n", "# for i1 in range(Sz[0]):\n", "# for j1 in range(Sz[1]):\n", "# print('i1', i1)\n", "# print('j1', j1)\n", "# print('int', str_out)\n", "# M_out[i1, j1] = int(str_out[ct1])\n", "# ct1 = ct1 + 1\n", "\n", "# return M_out" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/compute_lll.py\n", "{'mat': array([[ 0, 1],\n", " [ 1, 0],\n", " [-1, -1]])}\n", "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'')\n", "result.stdout b''\n", "str_out []\n", "[[0 0]\n", " [0 0]\n", " [0 0]]\n", "i1 0\n", "j1 0\n", "int []\n" ] }, { "ename": "IndexError", "evalue": "list index out of range", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0ml_bpb_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbp_basis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhkl_p\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m\u001b[0m in \u001b[0;36mbp_basis\u001b[0;34m(miller_ind)\u001b[0m\n\u001b[1;32m 73\u001b[0m \u001b[0;31m# The reduced basis vectors for the plane\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 74\u001b[0m \u001b[0mTmat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumn_stack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbv1_g1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbv2_g1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'int64'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 75\u001b[0;31m \u001b[0ml_pl_g1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlll_reduction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTmat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 76\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0ml_pl_g1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 77\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36mlll_reduction\u001b[0;34m(int_mat)\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[0minp_args\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 84\u001b[0m \u001b[0minp_args\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mat'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint_mat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 85\u001b[0;31m \u001b[0mlll_int_mat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcall_sage_math\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexec_str\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minp_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 86\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mlll_int_mat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m\u001b[0m in \u001b[0;36mcall_sage_math\u001b[0;34m(exec_str, inp_args)\u001b[0m\n\u001b[1;32m 124\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'j1'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 125\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'int'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr_out\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 126\u001b[0;31m \u001b[0mM_out\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstr_out\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mct1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 127\u001b[0m \u001b[0mct1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mct1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 128\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mIndexError\u001b[0m: list index out of range" ] } ], "source": [ "# import inspect\n", "# import byxtal\n", "# import subprocess\n", "# l_p_po = l1.l_p_po\n", "# l_po_p = nla.inv(l_p_po)\n", "\n", "# ## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)\n", "# l_bpb_p = bp_basis(hkl_p)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "ename": "IndexError", "evalue": "list index out of range", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;31m## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0ml_bpb_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbpb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbp_basis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhkl_p\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0ml_bpb_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ml_bpb_p\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m## l_bpb_p: Primitive Basis vectors of the boundary-plane (in po reference frame)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/Leila_sndhard/codes/byxtal/byxtal/byxtal/bp_basis.py\u001b[0m in \u001b[0;36mbp_basis\u001b[0;34m(miller_ind)\u001b[0m\n\u001b[1;32m 205\u001b[0m \u001b[0;31m# The reduced basis vectors for the plane\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 206\u001b[0m \u001b[0mTmat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumn_stack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbv1_g1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbv2_g1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'int64'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 207\u001b[0;31m \u001b[0ml_pl_g1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlll_reduction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mTmat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 208\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0ml_pl_g1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 209\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/Leila_sndhard/codes/byxtal/byxtal/byxtal/lll_tools.py\u001b[0m in \u001b[0;36mlll_reduction\u001b[0;34m(int_mat)\u001b[0m\n\u001b[1;32m 107\u001b[0m \u001b[0minp_args\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 108\u001b[0m \u001b[0minp_args\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mat'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint_mat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m \u001b[0mlll_int_mat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrpl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcall_sage_math\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexec_str\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minp_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 110\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mlll_int_mat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 111\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/Leila_sndhard/codes/byxtal/byxtal/byxtal/reduce_po_lat.py\u001b[0m in \u001b[0;36mcall_sage_math\u001b[0;34m(exec_str, inp_args)\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi1\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mSz\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mj1\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mSz\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m \u001b[0mM_out\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstr_out\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mct1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 39\u001b[0m \u001b[0mct1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mct1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 40\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mIndexError\u001b[0m: list index out of range" ] } ], "source": [ "l_p_po = l1.l_p_po\n", "l_po_p = nla.inv(l_p_po)\n", "\n", "## l_bpb_p: Primitive Basis vectors of the boundary-plane (in p reference frame)\n", "l_bpb_p = bpb.bp_basis(hkl_p)\n", "l_bpb_p = l_bpb_p.astype(int)\n", "## l_bpb_p: Primitive Basis vectors of the boundary-plane (in po reference frame)\n", "l_bpb_po = l_p_po.dot(l_bpb_p)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. `sig_num`: $\\Sigma$ number\n", "2. `l_bpbSig_p` : $\\Lambda_{\\Sigma-\\mathrm{2D}}^{p}$\n", "3. `l_bpbSig_po` : $\\Lambda_{\\Sigma-\\mathrm{2D}}^{po}$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'l_bpb_po' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m## area_bpl: Area of the 2D-primitive-unit-cell\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0marea_bpl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnla\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcross\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ml_bpb_po\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ml_bpb_po\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0msig_num\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mceil\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma_cut\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0marea_bpl\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mind2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'int64'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'l_bpb_po' is not defined" ] } ], "source": [ "## area_bpl: Area of the 2D-primitive-unit-cell\n", "area_bpl = nla.norm((np.cross(l_bpb_po[:,0], l_bpb_po[:,1])))\n", "sig_num = np.ceil(a_cut/area_bpl)\n", "\n", "ind2 = np.array([], dtype='int64');\n", "while (np.size(ind2) == 0):\n", " # Generate 2D hermite normal forms for sig_num (hnf_mats)\n", " hnf_mats = ghi.sig_hnf_mats(sig_num)\n", " # Compute the properties of the sub-lattices\n", " l_sig_p_mats, l_sig_po_mats = ghi.compute_hnf_props(hnf_mats, l_bpb_p, l_p_po, 1e-2)\n", " # Get the index for the sub-lattice that has the minimum cost\n", " ind2 = ghi.ind_min_cost(l_sig_po_mats, r_cut)\n", " sig_num = sig_num + 1\n", "\n", "## l_bpbSig_p: Basis vectors of the sub-lattice of the boundary-plane (in p reference frame)\n", "l_bpbSig_p = l_sig_p_mats[ind2];\n", "l_bpbSig_p = l_bpbSig_p.astype(int)\n", "l_bpbSig_po = l_p_po.dot(l_bpbSig_p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l_sig_p_mats" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. `l_po_go`: $\\Lambda_{PO}^{GO}$ - Orientation of crystal with respect to global reference frame $\\mathrm{GO}$\n", " + `go`: Orthogonal global reference frame (for LAMMPS)\n", " + `po`: Orhtogonal reference frame of the lattice\n", " + 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.\n", "\n", "2. `l2D_bpbSig_go`: $\\Lambda_{\\Sigma-\\mathrm{2D}}^{go}$\n", "3. `l2D_bpb_go`: $\\Lambda_{\\Sigma-\\mathrm{2D}}^{go}$\n", "\n", "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.\n", "\n", "4. `twoD_mat`: The $x$, $y$ components of the 2D basis vectors of $\\Lambda_{\\mathrm{2D}}^{go}$\n", "5. `twoDSig_mat`: The $x$, $y$ components of the 2D basis vectors of $\\Lambda_{\\Sigma-\\mathrm{2D}}^{go}$\n", "6. `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.\n", "7. `twoD_pts`: The lattice points in the 2D periodic box for the 3D slab." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l_po_go = uf.compute_orientation(l_bpbSig_po)\n", "l2D_bpbSig_go = l_po_go.dot(l_bpbSig_po)\n", "l2D_bpb_go = l_po_go.dot(l_bpb_po)\n", "twoD_mat = l2D_bpb_go[:2,:]\n", "twoDSig_mat = l2D_bpbSig_go[:2,:]\n", "r_cut1 = uf.compute_rCut(l_bpbSig_po)\n", "\n", "## Create 2D periodic box\n", "twoD_pts = uf.replicate_pts(twoD_mat, r_cut1)\n", "twoD_pts = uf.remove_periodic_overlaps(twoD_pts, twoDSig_mat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. `avec`: Basis vector $\\vec{a}$ of the 3D lattice defined using $\\mathcal{B}_{P2}$.\n", "2. `bvec`: Basis vector $\\vec{b}$ of the 3D lattice defined using $\\mathcal{B}_{P2}$.\n", "3. `l_p2_p1`: Equivalent basis for the 3D lattice, $\\Lambda_{P2}^{P1}$. \n", " + The function `find_int_solns` computes the third basis-vector, $\\vec{c}$, such that $\\mathcal{B}_{P1} \\equiv \\mathcal{B}_{P2}$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "avec = l_bpb_p[:,0]\n", "bvec = l_bpb_p[:,1]\n", "l_p2_p1 = uf.find_int_solns(avec, bvec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l2D_bpbSig_go" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. `l_p2_po1`: $\\Lambda_{p2}^{po1}$\n", "2. `l_p2_go`: $\\Lambda_{p2}^{go}$\n", "3. `tz_vec`: The components of $\\vec{c}$ of $\\mathcal{B}_{P2}$ in $go$ reference frame.\n", "4. `threeD_pts`: The 3D slab of atoms for simulating free surfaces.\n", " + This is obtained by replicating `twoD_pts` along the `tz_vec` in both +$z$ and -$z$ directions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l_p2_po1 = l_p_po.dot(l_p2_p1)\n", "l_p2_go = l_po_go.dot(l_p2_po1)\n", "tz_vec = np.array(l_p2_go[:,2], dtype='double')\n", "tz_vec = np.reshape(tz_vec, (3,))\n", "\n", "################################################################################\n", "## Translate 2D points in the Z-direction with zCut\n", "num_rep = np.abs(int(np.ceil(zCut/tz_vec[2])))\n", "num_2d = np.shape(twoD_pts)[0]\n", "num_3d_pts = int((2*num_rep+1)*num_2d)\n", "threeD_pts = np.zeros((num_3d_pts,3));\n", "\n", "twoD_pts1 = np.hstack((twoD_pts, np.zeros((num_2d,1))));\n", "\n", "for ct1 in np.arange(-num_rep, num_rep+1):\n", " ct2 = ct1 + num_rep\n", " ind_st = (ct2)*num_2d\n", " ind_stop = ind_st + num_2d\n", " trans_vec = tz_vec*ct1\n", " threeD_pts[ind_st:ind_stop, :] = twoD_pts1 + np.tile(trans_vec, (num_2d,1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The simulation box vectors are defined using Ovito's convenction.\n", "1. `sim_avec`: The components of $\\vec{a}$ in the global reference, $go$, of LAMMPS. $\\vec{a}$ is along $\\hat{e}_x$.\n", "2. `sim_bvec`: The components of $\\vec{b}$ in the global reference, $go$, of LAMMPS. $\\vec{b}$ lies in the $x-y$ plane.\n", "3. `sim_cvec`: The components of $\\vec{c}$ in the global reference, $go$, of LAMMPS.\n", "4. `sim_orig`: The origin of the simulation box." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "### Simulation Cell Box\n", "### Following Ovito's convention\n", "l_bpbSig_po_arr = np.array(l2D_bpbSig_go, dtype='double')\n", "sim_cell = np.zeros((3,4))\n", "sim_avec = l_bpbSig_po_arr[:,0]\n", "sim_bvec = l_bpbSig_po_arr[:,1]\n", "\n", "### Change this with inter-planar spacing\n", "sim_cvec = np.array([0,0,2*zCut]); # sim_cvec = np.array([0,0,zCut]);\n", "sim_orig = np.array([0,0,-zCut]); # sim_orig = np.array([0,0,0]);\n", "\n", "sim_cell[:,0] = sim_avec\n", "sim_cell[:,1] = sim_bvec\n", "sim_cell[:,2] = sim_cvec\n", "sim_cell[:,3] = sim_orig\n", "\n", "box_vecs = sim_cell[:,0:3]\n", "threeD_pts1 = uf.wrap_cc(sim_cell, threeD_pts)\n", "# tpts1 = np.dot(np.linalg.inv(box_vecs), threeD_pts.transpose()).transpose();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot of 3D lattice points along with the simulation box." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import plotting_routines as plr\n", "%matplotlib inline\n", "\n", "fig1 = plt.figure()\n", "plr.plot_3d_pts_box(fig1, threeD_pts1, sim_cell[:,0:3], sim_orig)\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.9" } }, "nbformat": 4, "nbformat_minor": 4 }