MD2NEB-clamp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surface restructuring event analysis of LAMMPS MD trajectory:\n",
    "## Preparation of LAMMPS NEB calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lmpdump as lmpdump    # Use lmpdump.py parser\n",
    "import os\n",
    "import math\n",
    "import numpy as np\n",
    "import time\n",
    "from tqdm import tqdm\n",
    "from sitator.util import PBCCalculator\n",
    "\n",
    "pwd = os.getcwd()\n",
    "\n",
    "print('Please cite DOI: 10.1021/acs.jpcc.9b04863.')\n",
    "print('This script detects interlayer surface restructuring events from LAMMPS MD trajectory and generates images for subsequent LAMMPS NEB calculation of each event.')\n",
    "print('This script is interactive and requires user input. Please read carefully before proceeding:\\n')\n",
    "\n",
    "print('Note 1: Intended only for interlayer surface restructuring events in FCC metals.')\n",
    "print('Note 2: Intended only for periodic slab models with vacuum along the z-direction.')\n",
    "print('Note 3: Intended only for LAMMPS NVT simulations.')\n",
    "print('Note 4: Intended only for estimation of the event statistics, up to bimetallic systems. Please confirm with visual inspection as needed.\\n')\n",
    "print('Note 5: Must have used DFT unit cell if planning on subsequent DFT optimization (via NEB2DFT.ipynb).')\n",
    "\n",
    "print('Note 5: Requires a LAMMPS DATA file containing the equilibrated box dimensions.')\n",
    "print('Note 6: Requires two trajectories, ordered by atom ID & unscaled: (1) Clamped (via lmpclamp.py), dumped every ps, xy-unwrapped; (2) Raw, dumped every 0.1 ps, wrapped.')\n",
    "print('Note 7: Requires lmpdump.py to load the LAMMPS trajectories.\\n')\n",
    "print('Note 8: Requires PBC calculator from sitator package.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LAMMPS trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LAMMPS custom-style dump file - Clamped trajectory\n",
    "name = input('Enter the name of the clamped LAMMPS trajectory file (ID, type, xu, yu, z, c_cn) dumped every ps : \\nWarning: Must be clamped (via lmpclamp.py), ordered by atom ID, unscaled & xy-unwrapped! Please quit now if not so.\\n')\n",
    "print('Loading trajectory-1...')\n",
    "file = pwd + '/' + name\n",
    "xyz1 = lmpdump.lmpdump(file, loadmode='all')\n",
    "print('\\n')\n",
    "\n",
    "print('Loading time steps...')\n",
    "step = []\n",
    "for s in xyz1.finaldict.keys():    # keys = time steps\n",
    "    step.append(s)\n",
    "print('Done.\\n')\n",
    "\n",
    "size = np.array(step).size    # Number of time steps\n",
    "N_dump = xyz1.finaldict[0][1].shape[0]    # Number of atoms dumped\n",
    "xlo_ext = xyz1.finaldict[0][0][0]    # x_ext = x + xy\n",
    "xhi_ext = xyz1.finaldict[0][0][1]\n",
    "\n",
    "# LAMMPS atom-style dump file - Raw trajectory\n",
    "name = input('\\n Enter the name of the raw LAMMPS trajectory file (ID, type, x, y, z) dumped every 0.1 ps : \\nWarning: Must be raw, ordered by atom ID, unscaled & wrapped! Please quit now if not so.\\n')\n",
    "print('Loading trajectory-2...')\n",
    "file = pwd + '/' + name\n",
    "xyz2 = lmpdump.lmpdump(file, loadmode='all')\n",
    "print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract unit cell information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load LAMMPS data file after NPT equilibration\n",
    "name = input('Enter the name of the LAMMPS DATA file containing the equilibrated box dimensions : ')\n",
    "file_eq = pwd + '/' + name\n",
    "log = open(file_eq,'r')\n",
    "log_lines = log.readlines()\n",
    "log_lines = [line.split() for line in log_lines]\n",
    "\n",
    "for line_index, line in enumerate(log_lines):\n",
    "    if line:\n",
    "        for l in line:\n",
    "            \n",
    "            if l == 'atoms':\n",
    "                N_atom = int(line[0])\n",
    "                \n",
    "            elif l == 'types':\n",
    "                N_type = int(line[0])\n",
    "\n",
    "            elif l == 'xlo':\n",
    "                xlo = float(line[0])\n",
    "                xhi = float(line[1])\n",
    "            \n",
    "            elif l == 'ylo':\n",
    "                ylo = float(line[0])\n",
    "                yhi = float(line[1])\n",
    "            \n",
    "            elif l == 'zlo':\n",
    "                zlo = float(line[0])\n",
    "                zhi = float(line[1])\n",
    "            \n",
    "            elif l == 'xy':\n",
    "                xy = float(line[0])\n",
    "                xz = float(line[1])\n",
    "                yz = float(line[2])\n",
    "            \n",
    "            elif l == 'Atoms':\n",
    "                head = line_index + 1\n",
    "                start = line_index + 2\n",
    "            \n",
    "            elif l == 'Velocities':\n",
    "                end = line_index - 1\n",
    "\n",
    "Lx = xhi - xlo    # Box dimensions\n",
    "Ly = yhi - ylo\n",
    "Lz = zhi - zlo\n",
    "tan = Ly/xy\n",
    "\n",
    "xyz = []    # N*3 matrix\n",
    "for line in log_lines[start:end]:\n",
    "\n",
    "    line[:2] = np.array(line[:2]).astype(int)    # Atom ID, atom type\n",
    "    line[2:5] = np.array(line[2:5]).astype(float)    # x, y, z\n",
    "    line[5:8] = np.array(line[5:8]).astype(int)    # ix, iy, iz (inverse sign)\n",
    "\n",
    "    x = line[2] + line[5]*Lx + line[6]*xy    # Wrapped coordinates\n",
    "    y = line[3] + line[6]*Ly\n",
    "    z = line[4] + line[7]*Lz\n",
    "    \n",
    "    xyz.append([x,y,z])\n",
    "xyz = np.array(xyz)\n",
    "\n",
    "ID2 = int(input('Enter the ID of a nearest-neighbor atom that lies in the same layer as atom #1 : '))\n",
    "ID3 = int(input('Enter the ID of another nearest-neighbor atom that lies in the same layer as atom #1 : '))\n",
    "ID4 = int(input('Enter the ID of an atom that lies one layer above atom #1 : '))\n",
    "\n",
    "for atom, coord in enumerate(xyz):\n",
    "    \n",
    "    if atom+1 == 1:\n",
    "        atom1 = coord    # atom1 = Reference atom\n",
    "    \n",
    "    elif atom+1 == ID2:\n",
    "        atom2 = coord    # atom2 = 1st coplanar atom\n",
    "\n",
    "    elif atom+1 == ID3:\n",
    "        atom3 = coord    # atom3 = 2nd coplanar atom\n",
    "        \n",
    "    elif atom+1 == ID4:\n",
    "        atom4 = coord    # atom4 = Atom one layer above\n",
    "\n",
    "r12 = np.linalg.norm(np.subtract(atom1,atom2))\n",
    "r13 = np.linalg.norm(np.subtract(atom1,atom3))\n",
    "dr_avg = np.mean([r12, r13])    # In-plane NN distance\n",
    "\n",
    "a1 = np.subtract(atom2,atom1)    # 1st intralayer vector\n",
    "a1 = a1 / np.linalg.norm(a1)\n",
    "\n",
    "a2 = np.subtract(atom3,atom1)    # 2nd intralayer vector\n",
    "a2 = a2 / np.linalg.norm(a2)\n",
    "\n",
    "a3 = np.cross(a1,a2)    # Normal vector\n",
    "if a3[2] < 0:\n",
    "    a3 = -a3    # Ensure normal vector along +z\n",
    "a3 = a3 / np.linalg.norm(a3)\n",
    "\n",
    "dn_avg = np.dot(atom4,a3) - np.dot(atom1,a3)    # Interlayer distance; n = normal component\n",
    "\n",
    "N_slab = int(input('Enter the number of layers in the slab, excluding the deposit/adsorbate layer(s) : '))\n",
    "if xyz1.finaldict[0][1]['id'][0] != 1:    # Is bottommost layer muted?\n",
    "    N_slab -= 1\n",
    "\n",
    "z_layer = []    # Clamped layer heights\n",
    "for atom in range(1,N_dump):\n",
    "    z = xyz1.finaldict[step[-1]][1]['z'][atom]    # Use last frame\n",
    "    if z not in z_layer:\n",
    "        z_layer.append(z)\n",
    "N_layer = len(z_layer)\n",
    "\n",
    "log.close()\n",
    "\n",
    "print('Unit cell information & normal direction extracted.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Parsing trajectory for restructuring events...')\n",
    "t_begin = time.perf_counter()\n",
    "\n",
    "file_ev = pwd + '/event_flag.txt'    # List of time steps flagged as event points\n",
    "out_ev = open(file_ev,'w')\n",
    "\n",
    "header1 = ' tf = final time (flagged step) \\n dt = event duration \\n ni = initial normal component \\n nf = final normal \\n CNi = initial coordination number \\n CNf = final coordination number \\n \\n'\n",
    "header2 = 'tf (ps) \\t dt (ps) \\t ID \\t type \\t ni \\t nf \\t CNi \\t CNf \\n'\n",
    "out_ev.write(header1)\n",
    "out_ev.write(header2)\n",
    "\n",
    "skip_index = [0, 1, 2, size-2, size-1]    # Skip time steps at the edges of the trajectory\n",
    "\n",
    "for s_index, s in enumerate(tqdm(step)):\n",
    "    \n",
    "    if s_index in skip_index:\n",
    "        continue\n",
    "\n",
    "    ev_atom = []    # List of active atoms detected\n",
    "    for atom in range(0,N_dump):\n",
    "            \n",
    "        # f0 = current time step (final 0)\n",
    "        # f1 = 1 step after (final 1)\n",
    "        # f2 = 2 steps after (final 2)\n",
    "        # i1 = 1 step before (initial 1)\n",
    "        # i2 = 2 steps before (initial 2)\n",
    "        # i3 = 3 steps before (initial 3)\n",
    "        \n",
    "        ID = xyz1.finaldict[s][1]['id'][atom]\n",
    "        Type = xyz1.finaldict[s][1]['type'][atom]\n",
    "        \n",
    "        CNf0 = xyz1.finaldict[s][1]['c_cn'][atom]\n",
    "        CNi1 = xyz1.finaldict[step[s_index-1]][1]['c_cn'][atom]\n",
    "        CNi2 = xyz1.finaldict[step[s_index-2]][1]['c_cn'][atom]\n",
    "        \n",
    "        xf0 = xyz1.finaldict[s][1]['xu'][atom]\n",
    "        yf0 = xyz1.finaldict[s][1]['yu'][atom]\n",
    "        zf0 = xyz1.finaldict[s][1]['z'][atom]\n",
    "        rf0 = np.array([xf0, yf0, zf0])\n",
    "        nf0 = np.dot(rf0,a3)\n",
    "        xf0_prj = xf0 - yf0/tan    # Tilted projection of x\n",
    "        Nxf0 = math.floor(xf0_prj/Lx)    # x image number\n",
    "        Nyf0 = math.floor(yf0/Ly)    # y image number\n",
    "    \n",
    "        xf1 = xyz1.finaldict[step[s_index+1]][1]['xu'][atom]\n",
    "        yf1 = xyz1.finaldict[step[s_index+1]][1]['yu'][atom]\n",
    "        zf1 = xyz1.finaldict[step[s_index+1]][1]['z'][atom]\n",
    "        rf1 = np.array([xf1, yf1, zf1])\n",
    "        nf1 = np.dot(rf1,a3)\n",
    "        xf1_prj = xf1 - yf1/tan\n",
    "        Nxf1 = math.floor(xf1_prj/Lx)\n",
    "        Nyf1 = math.floor(yf1/Ly)\n",
    "        \n",
    "        xf2 = xyz1.finaldict[step[s_index+2]][1]['xu'][atom]\n",
    "        yf2 = xyz1.finaldict[step[s_index+2]][1]['yu'][atom]\n",
    "        zf2 = xyz1.finaldict[step[s_index+2]][1]['z'][atom]\n",
    "        rf2 = np.array([xf2, yf2, zf2])\n",
    "        nf2 = np.dot(rf2,a3)\n",
    "        xf2_prj = xf2 - yf2/tan\n",
    "        Nxf2 = math.floor(xf2_prj/Lx)\n",
    "        Nyf2 = math.floor(yf2/Ly)\n",
    "        \n",
    "        xi1 = xyz1.finaldict[step[s_index-1]][1]['xu'][atom]\n",
    "        yi1 = xyz1.finaldict[step[s_index-1]][1]['yu'][atom]\n",
    "        zi1 = xyz1.finaldict[step[s_index-1]][1]['z'][atom]\n",
    "        ri1 = np.array([xi1, yi1, zi1])\n",
    "        ni1 = np.dot(ri1,a3)\n",
    "        xi1_prj = xi1 - yi1/tan\n",
    "        Nxi1 = math.floor(xi1_prj/Lx)\n",
    "        Nyi1 = math.floor(yi1/Ly)\n",
    "        \n",
    "        xi2 = xyz1.finaldict[step[s_index-2]][1]['xu'][atom]\n",
    "        yi2 = xyz1.finaldict[step[s_index-2]][1]['yu'][atom]\n",
    "        zi2 = xyz1.finaldict[step[s_index-2]][1]['z'][atom]\n",
    "        ri2 = np.array([xi2, yi2, zi2])\n",
    "        ni2 = np.dot(ri2,a3)\n",
    "        xi2_prj = xi2 - yi2/tan\n",
    "        Nxi2 = math.floor(xi2_prj/Lx)\n",
    "        Nyi2 = math.floor(yi2/Ly)\n",
    "        \n",
    "        xi3 = xyz1.finaldict[step[s_index-3]][1]['xu'][atom]\n",
    "        yi3 = xyz1.finaldict[step[s_index-3]][1]['yu'][atom]\n",
    "        zi3 = xyz1.finaldict[step[s_index-3]][1]['z'][atom]\n",
    "        ri3 = np.array([xi3, yi3, zi3])\n",
    "        ni3 = np.dot(ri3,a3)\n",
    "        xi3_prj = xi3 - yi3/tan\n",
    "        Nxi3 = math.floor(xi3_prj/Lx)\n",
    "        Nyi3 = math.floor(yi3/Ly)\n",
    "    \n",
    "        dni1 = nf0 - ni1    # Change in normal from the past\n",
    "        dni2 = nf0 - ni2\n",
    "        dni3 = nf0 - ni3\n",
    "        \n",
    "        dnf1 = nf1 - nf0    # Change in normal in the future\n",
    "        dnf2 = nf2 - nf0\n",
    "        \n",
    "        # Test 1 step in the past\n",
    "        # Criterion: Change by more than 90% of interlayer distance\n",
    "        # History check: Must not revert back 1 or 2 steps in the past & in the future\n",
    "        \n",
    "        if (abs(dni1) > 0.90*dn_avg) and (abs(dni2) > 0.90*dn_avg) and (abs(dnf1) < 0.10*dn_avg) and (abs(dnf2) < 0.10*dn_avg):\n",
    "            \n",
    "            ev_atom.append(atom)\n",
    "            dt = 1\n",
    "            \n",
    "            xf0 -= xy * Nyf0    # Wrap y first\n",
    "            yf0 -= Ly * Nyf0\n",
    "            \n",
    "            xi1 -= xy * Nyi1\n",
    "            yi1 -= Ly * Nyi1\n",
    "            \n",
    "            if (Nxf0 != 0) or (Nxi1 != 0):\n",
    "            \n",
    "                xf0 -= Lx * np.maximum(Nxf0,Nxi1)    # Shift initial & final x together toward the unit cell\n",
    "                rf0 = np.array([xf0, yf0, zf0])\n",
    "                nf0 = np.dot(rf0,a3)\n",
    "            \n",
    "                xi1 -= Lx * np.maximum(Nxf0,Nxi1)\n",
    "                ri1 = np.array([xi1, yi1, zi1])\n",
    "                ni1 = np.dot(ri1,a3)\n",
    "            \n",
    "            line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni1, nf0, CNi1, CNf0)\n",
    "            out_ev.write(line)\n",
    "            \n",
    "            for aux in range(0,N_dump):    # Detect any auxiliary atom that may move together\n",
    "                if aux not in ev_atom:    # Test only if not flagged before\n",
    "                    \n",
    "                    # Auxiliary atom labeled with \"i\" & \"f\" only\n",
    "                    \n",
    "                    ID = xyz1.finaldict[s][1]['id'][aux]\n",
    "                    Type = xyz1.finaldict[s][1]['type'][aux]\n",
    "                    \n",
    "                    CNf = xyz1.finaldict[s][1]['c_cn'][aux]\n",
    "                    CNi = xyz1.finaldict[step[s_index-1]][1]['c_cn'][aux]\n",
    "    \n",
    "                    xf = xyz1.finaldict[s][1]['xu'][aux]\n",
    "                    yf = xyz1.finaldict[s][1]['yu'][aux]\n",
    "                    zf = xyz1.finaldict[s][1]['z'][aux]\n",
    "                    rf = np.array([xf, yf, zf])\n",
    "                    nf = np.dot(rf,a3)\n",
    "                    xf_prj = xf - yf/tan\n",
    "                    Nxf = math.floor(xf_prj/Lx)\n",
    "                    Nyf = math.floor(yf/Ly)\n",
    "                    \n",
    "                    xi = xyz1.finaldict[step[s_index-1]][1]['xu'][aux]\n",
    "                    yi = xyz1.finaldict[step[s_index-1]][1]['yu'][aux]\n",
    "                    zi = xyz1.finaldict[step[s_index-1]][1]['z'][aux]\n",
    "                    ri = np.array([xi, yi, zi])\n",
    "                    ni = np.dot(ri,a3)\n",
    "                    xi_prj = xi - yi/tan\n",
    "                    Nxi = math.floor(xi_prj/Lx)\n",
    "                    Nyi = math.floor(yi/Ly)\n",
    "    \n",
    "                    dx = xf - xi\n",
    "                    dy = yf - yi\n",
    "                    dz = zf - zi                        \n",
    "                    dr = math.sqrt(dx**2 + dy**2 + dz**2)    # Change in position, rather than normal\n",
    "                    \n",
    "                    if dr > 0.80*dr_avg:    # Criterion: Change by more than 80% of the NN distance\n",
    "                        ev_atom.append(aux)\n",
    "                        \n",
    "                        # Auxiliary atom must be wrapped toward the active atom\n",
    "\n",
    "                        xf -= xy * Nyf    # Wrap y first\n",
    "                        yf -= Ly * Nyf\n",
    "                        \n",
    "                        xi -= xy * Nyi\n",
    "                        yi -= Ly * Nyi\n",
    "                        \n",
    "                        dxf = (xf - xf0)/Lx    # Extent of x deviation from the active atom\n",
    "                        dxi = (xi - xi1)/Lx\n",
    "                        \n",
    "                        if (abs(dxf) > 0.70) or (abs(dxi) > 0.70):    # Wrap if x deviation is more than 70% of the unit cell x dimension\n",
    "                                \n",
    "                            xf -= Lx * round(dxf)\n",
    "                            rf = np.array([xf, yf, zf])\n",
    "                            nf = np.dot(rf,a3)\n",
    "                            \n",
    "                            xi -= Lx * round(dxi)\n",
    "                            ri = np.array([xi, yi, zi])\n",
    "                            ni = np.dot(ri,a3)\n",
    "\n",
    "                        line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ev.write(line)\n",
    "        \n",
    "        # If no change 1 step in the past, test 2 steps in the past (change may be gradual)\n",
    "        elif (abs(dni2) > 0.90*dn_avg) and (abs(dni3) > 0.90*dn_avg) and (abs(dnf1) < 0.10*dn_avg) and (abs(dnf2) < 0.10*dn_avg): \n",
    "            \n",
    "            ev_atom.append(atom)\n",
    "            dt = 2\n",
    "            \n",
    "            xf0 -= xy * Nyf0\n",
    "            yf0 -= Ly * Nyf0\n",
    "            \n",
    "            xi2 -= xy * Nyi2\n",
    "            yi2 -= Ly * Nyi2\n",
    "    \n",
    "            if (Nxf0 != 0) or (Nxi2 != 0):\n",
    "            \n",
    "                xf0 -= Lx * np.maximum(Nxf0,Nxi2)\n",
    "                rf0 = np.array([xf0, yf0, zf0])\n",
    "                nf0 = np.dot(rf0,a3)\n",
    "            \n",
    "                xi2 -= Lx * np.maximum(Nxf0,Nxi2)\n",
    "                ri2 = np.array([xi2, yi2, zi2])\n",
    "                ni2 = np.dot(ri2,a3)\n",
    "            \n",
    "            line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni2, nf0, CNi2, CNf0)\n",
    "            out_ev.write(line)\n",
    "    \n",
    "            for aux in range(0,N_dump):\n",
    "                if aux not in ev_atom:\n",
    "                    \n",
    "                    ID = xyz1.finaldict[s][1]['id'][aux]\n",
    "                    Type = xyz1.finaldict[s][1]['type'][aux]\n",
    "                    \n",
    "                    CNf = xyz1.finaldict[s][1]['c_cn'][aux]\n",
    "                    CNi = xyz1.finaldict[step[s_index-1]][1]['c_cn'][aux]\n",
    "    \n",
    "                    xf = xyz1.finaldict[s][1]['xu'][aux]                            \n",
    "                    yf = xyz1.finaldict[s][1]['yu'][aux]\n",
    "                    zf = xyz1.finaldict[s][1]['z'][aux]\n",
    "                    rf = np.array([xf, yf, zf])\n",
    "                    nf = np.dot(rf,a3)\n",
    "                    xf_prj = xf - yf/tan\n",
    "                    Nxf = math.floor(xf_prj/Lx)\n",
    "                    Nyf = math.floor(yf/Ly)\n",
    "      \n",
    "                    xi = xyz1.finaldict[step[s_index-2]][1]['xu'][aux]\n",
    "                    yi = xyz1.finaldict[step[s_index-2]][1]['yu'][aux]\n",
    "                    zi = xyz1.finaldict[step[s_index-2]][1]['z'][aux]\n",
    "                    ri = np.array([xi, yi, zi])\n",
    "                    ni = np.dot(ri,a3)\n",
    "                    xi_prj = xi - yi/tan\n",
    "                    Nxi = math.floor(xi_prj/Lx)\n",
    "                    Nyi = math.floor(yi/Ly)\n",
    "    \n",
    "                    dx = xf - xi\n",
    "                    dy = yf - yi\n",
    "                    dz = zf - zi                        \n",
    "                    dr = math.sqrt(dx**2 + dy**2 + dz**2)\n",
    "                    \n",
    "                    if dr > 0.80*dr_avg:\n",
    "                        ev_atom.append(aux)\n",
    "                        \n",
    "                        xf -= xy * Nyf\n",
    "                        yf -= Ly * Nyf\n",
    "                        \n",
    "                        xi -= xy * Nyi\n",
    "                        yi -= Ly * Nyi\n",
    "                        \n",
    "                        dxf = (xf - xf0)/Lx\n",
    "                        dxi = (xi - xi2)/Lx\n",
    "                        \n",
    "                        if (abs(dxf) > 0.70) or (abs(dxi) > 0.70):\n",
    "\n",
    "                            xf -= Lx * round(dxf)\n",
    "                            rf = np.array([xf, yf, zf])\n",
    "                            nf = np.dot(rf,a3)\n",
    "\n",
    "                            xi -= Lx * round(dxi)\n",
    "                            ri = np.array([xi, yi, zi])\n",
    "                            ni = np.dot(ri,a3)\n",
    "                            \n",
    "                        line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ev.write(line)\n",
    "                            \n",
    "out_ev.close()\n",
    "\n",
    "t_end = time.perf_counter()\n",
    "print('Event detection done in %.2f' % (t_end - t_begin) + 's > event_flag.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = open(file_ev,'r')    # Load list of time steps flagged as event points\n",
    "ev_lines = ev.readlines()\n",
    "ev_lines = [line.split() for line in ev_lines]\n",
    "\n",
    "for line_index, line in enumerate(ev_lines):\n",
    "    if line_index > 7:\n",
    "        line[0:4] = np.array(line[0:4]).astype(int)    # Time step, duration, atom ID, atom type\n",
    "        line[4:6] = np.array(line[4:6]).astype(float)    # Initial normal, final normal\n",
    "        line[6:8] = np.array(line[6:8]).astype(int)    # Initial coordination, final coordination\n",
    "\n",
    "ev.close()\n",
    "\n",
    "file_cluster = pwd + '/event_cluster.txt'    # List of events\n",
    "out_cluster = open(file_cluster,'w')\n",
    "\n",
    "# List active atoms first, followed by event #, final time, and duration\n",
    "header1 = ' tf = final time (flagged step) \\n dt = event duration \\n ni = initial normal component \\n nf = final normal component \\n CNi = initial coordination number \\n CNf = final coordination number \\n \\n'\n",
    "header2 = '\\t ID \\t type \\t ni \\t nf \\t CNi \\t CNf \\n event # \\t tf (ps) \\t dt (ps) \\n \\n'\n",
    "out_cluster.write(header1)\n",
    "out_cluster.write(header2)\n",
    "\n",
    "ev_no = 0\n",
    "ev_atom = []    # List of active atoms for each event\n",
    "ev_flag = []    # List of flags for each event\n",
    "\n",
    "for line_index, line in enumerate(ev_lines):\n",
    "\n",
    "    if line_index < 8:\n",
    "        continue\n",
    "    \n",
    "    elif line_index == 8:\n",
    "        \n",
    "        ev_flag.append(line)\n",
    "        ti = line[0] - line[1]    # ti = initial time\n",
    "        ID = line[2]\n",
    "        ev_atom.append(ID)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        # If a step is within 3 steps of the previously flagged step, it is part of the same event\n",
    "        if line[0] - ev_lines[line_index-1][0] < 3:\n",
    "            \n",
    "            ev_flag.append(line)\n",
    "            ID = line[2]\n",
    "            \n",
    "            if ID not in ev_atom:\n",
    "                ev_atom.append(ID)\n",
    "\n",
    "            if line_index == len(ev_lines)-1:    # If at the last entry, finish clustering\n",
    "                \n",
    "                ev_no += 1\n",
    "                tf = ev_lines[line_index-1][0]\n",
    "\n",
    "                if tf == ti:\n",
    "                    dt = ev_lines[line_index-1][1]\n",
    "                else:\n",
    "                    dt = tf - ti\n",
    "\n",
    "                for atom in ev_atom:\n",
    "                    atom_flag = []    # List of flags for each atom\n",
    "                    \n",
    "                    for flag in ev_flag:\n",
    "                        ID = flag[2]\n",
    "                        \n",
    "                        if ID == atom:\n",
    "                            atom_flag.append(flag)\n",
    "                    \n",
    "                    for flag_index, flag in enumerate(atom_flag):\n",
    "                        \n",
    "                        ID = flag[2]\n",
    "                        Type = flag[3]\n",
    "                        ni = flag[4]\n",
    "                        nf = flag[5]\n",
    "                        CNi = flag[6]\n",
    "                        CNf = flag[7]\n",
    "                        dn = nf - ni\n",
    "                        \n",
    "                        if abs(dn) < 0.10*dn_avg:    # If normal does not change, print only if at the last flag\n",
    "                            \n",
    "                            if flag_index != len(atom_flag)-1:\n",
    "                                continue\n",
    "                            \n",
    "                            else:\n",
    "                                out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                                out_cluster.write(out)\n",
    "                        \n",
    "                        elif abs(dn) > 0.90*dn_avg:    # If normal changes, print that flag and move on to the next atom\n",
    "                            out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                            out_cluster.write(out)\n",
    "                            break\n",
    "\n",
    "                out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)    # Print the event number, time step, and duration\n",
    "                out_cluster.write(out)\n",
    "                    \n",
    "        # If a step is after 3 steps or more from the previously flagged step, it is part of a different event, so finish clustering\n",
    "        else:\n",
    "\n",
    "            ev_no += 1\n",
    "            tf = ev_lines[line_index-1][0]\n",
    "\n",
    "            if tf == ti:\n",
    "                dt = ev_lines[line_index-1][1]\n",
    "            else:\n",
    "                dt = tf - ti\n",
    "\n",
    "            for atom in ev_atom:\n",
    "                atom_flag = []\n",
    "                \n",
    "                for flag in ev_flag:\n",
    "                    ID = flag[2]\n",
    "                    \n",
    "                    if ID == atom:\n",
    "                        atom_flag.append(flag)\n",
    "                \n",
    "                for flag_index, flag in enumerate(atom_flag):\n",
    "                    \n",
    "                    ID = flag[2]\n",
    "                    Type = flag[3]\n",
    "                    ni = flag[4]\n",
    "                    nf = flag[5]\n",
    "                    CNi = flag[6]\n",
    "                    CNf = flag[7]\n",
    "                    dn = nf - ni\n",
    "                    \n",
    "                    if abs(dn) < 0.10*dn_avg:\n",
    "                        \n",
    "                        if flag_index != len(atom_flag)-1:\n",
    "                            continue\n",
    "                        \n",
    "                        else:\n",
    "                            out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                            out_cluster.write(out)\n",
    "                    \n",
    "                    elif abs(dn) > 0.90*dn_avg:\n",
    "                        out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_cluster.write(out)\n",
    "                        break\n",
    "            \n",
    "            out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)\n",
    "            out_cluster.write(out)\n",
    "            \n",
    "            ev_atom = []    # Reset for the next event\n",
    "            ev_flag = []\n",
    "            \n",
    "            ev_flag.append(line)\n",
    "            ti = line[0] - line[1]\n",
    "            ID = line[2]\n",
    "            \n",
    "            if ID not in ev_atom:\n",
    "                ev_atom.append(ID)\n",
    "\n",
    "out_cluster.close()\n",
    "\n",
    "print('Event clustering done > event_cluster.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event localization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pbc = PBCCalculator(lat)    # Periodic boundary condition\n",
    "\n",
    "cluster = open(file_cluster,'r')    # Load event clusters\n",
    "cluster_lines = cluster.readlines()\n",
    "cluster_lines = [line.split() for line in cluster_lines]\n",
    "\n",
    "file_ls = pwd + '/event_list.txt'    # List of localized events\n",
    "out_ls = open(file_ls,'w')\n",
    "\n",
    "out_ls.write(header1)\n",
    "out_ls.write(header2)\n",
    "\n",
    "for line_index, line in enumerate(cluster_lines):\n",
    "    if line and (line_index > 10):\n",
    "        if len(line) == 6:\n",
    "            line[0:2] = np.array(line[0:2]).astype(int)    # Atom ID, atom type\n",
    "            line[2:4] = np.array(line[2:4]).astype(float)    # Initial normal, final normal\n",
    "            line[4:6] = np.array(line[4:6]).astype(int)    # Initial coordination, final coordination\n",
    "        else:\n",
    "            line[0:3] = np.array(line[0:3]).astype(int)    # event number, time step, duration\n",
    "\n",
    "line_cont = 10    # Line number continuation\n",
    "N_splice = 0    # Number of events spliced\n",
    "while line_cont < len(cluster_lines)-2:\n",
    "    \n",
    "    ev_atom = []    # List of event atoms\n",
    "    for line_index, line in enumerate(cluster_lines):\n",
    "        if line and (line_index > line_cont):\n",
    "            if len(line) == 6:\n",
    "                ev_atom.append(int(line[0]))\n",
    "            else:\n",
    "                ev_no = line[0]\n",
    "                tf = line[1]\n",
    "                dt = line[2]\n",
    "                break\n",
    "    \n",
    "    ev_coord = []    # List of Cartesian coordinates for event atoms\n",
    "    for ID in ev_atom:\n",
    "        for atom in range(0,N_dump):\n",
    "            if xyz1.finaldict[0][1]['id'][atom] == ID:\n",
    "                index = atom\n",
    "                x = xyz1.finaldict[step[tf]][1]['xu'][index]\n",
    "                y = xyz1.finaldict[step[tf]][1]['yu'][index]\n",
    "                z = xyz1.finaldict[step[tf]][1]['z'][index]\n",
    "                r = np.array([x, y, z])\n",
    "                ev_coord.append(r)\n",
    "                break\n",
    "    ev_coord = np.array(ev_coord)\n",
    "        \n",
    "    ev_rij = pbc.pairwise_distances(ev_coord)    # Pairwise distance matrix for event atoms\n",
    "    ev_connectivity = (ev_rij <= 2.0*dr_avg)    # Connectivity matrix (T/F) within a distance cutoff\n",
    "    \n",
    "    from scipy.sparse.csgraph import connected_components\n",
    "    N_group, ls_group = connected_components(ev_connectivity)    # N_group = number of groups; ls_group = list of groups for each event atom\n",
    "    \n",
    "    if N_group == 1:    # No splicing needed\n",
    "        \n",
    "        for line_index, line in enumerate(cluster_lines):\n",
    "            if line and (line_index > line_cont):\n",
    "                if len(line) == 6:\n",
    "                    out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(line[0], line[1], line[2], line[3], line[4], line[5])\n",
    "                    out_ls.write(out)\n",
    "                else:\n",
    "                    out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)\n",
    "                    out_ls.write(out)\n",
    "                    line_cont = line_index\n",
    "                    break\n",
    "    \n",
    "    else:    # Splicing needed\n",
    "        \n",
    "        N_splice += 1\n",
    "        for group_index, group in enumerate(range(0,N_group)):\n",
    "            \n",
    "            ev_group = []    # List of event atoms for each group\n",
    "            for atom_index, atom_group in ls_group:\n",
    "                if atom_group == group:\n",
    "                    ev_group.append(ev_atom[atom_index])\n",
    "            \n",
    "            for line_index, line in enumerate(cluster_lines):    # Check whether group is z-active\n",
    "                if line and (line_index > line_cont):\n",
    "                    if len(line) == 6:\n",
    "                        if line[0] in ev_group:\n",
    "                            if abs(atom[2] - atom[3]) > 0.50*dn_avg:\n",
    "                                z_active = True\n",
    "                                break\n",
    "                    else:\n",
    "                        z_active = False\n",
    "                        break\n",
    "            \n",
    "            if z_active:    # Write only if group is z-active\n",
    "            \n",
    "                for line_index, line in enumerate(cluster_lines):\n",
    "                    if line and (line_index > line_cont):\n",
    "                        if len(line) == 6:\n",
    "                            if line[0] in ev_group:\n",
    "                                out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(line[0], line[1], line[2], line[3], line[4], line[5])\n",
    "                                out_ls.write(out)\n",
    "                        else:\n",
    "                            ev_no_new = str(ev_no) + '-' + str(group_index+1)\n",
    "                            out = '%s\\t%i\\t%i\\n\\n'%(ev_no_new, tf, dt)\n",
    "                            out_ls.write(out)\n",
    "                            if group_index == N_group-1:\n",
    "                                line_cont = line_index\n",
    "                            break                \n",
    "\n",
    "cluster.close()\n",
    "out_ls.close()\n",
    "\n",
    "print('Event localization done: %i ' %(N_splice) + 'events spliced for simultaneous events > event_list.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_ls = pwd + '/event_list.txt'    # List of events\n",
    "ls = open(file_ls,'r')    # Load list of events\n",
    "ls_lines = ls.readlines()\n",
    "ls_lines = [line.split() for line in ls_lines]\n",
    "\n",
    "file_cls = pwd + '/event_class.txt'    # List of event classes\n",
    "out_cls = open(file_cls,'w')\n",
    "\n",
    "disclaimer = 'Intended only for estimation of the event statistics. Please confirm with visual inspection as needed.\\n'\n",
    "header = 'ID \\t Type \\t Site_i \\t Site_f \\n\\n'\n",
    "out_cls.write(disclaimer)\n",
    "out_cls.write(header)\n",
    "\n",
    "Lv1_1 = 'Isolated adatom'\n",
    "Lv1_2 = 'Edge/kink/vacancy'\n",
    "Lv1_3 = 'Deposit intralayer'\n",
    "Lv0_1 = 'Surface intralayer'\n",
    "Lv0_2 = 'Subdeposit'\n",
    "Lv_sub = 'Subsurface'\n",
    "\n",
    "for line_index, line in enumerate(ls_lines):\n",
    "    if line and (line_index > 9):\n",
    "        if len(line) == 6:\n",
    "            line[0:2] = np.array(line[0:2]).astype(int)    # Atom ID, atom type\n",
    "            line[2:4] = np.array(line[2:4]).astype(float)    # Initial normal, final normal\n",
    "            line[4:6] = np.array(line[4:6]).astype(int)    # Initial coordination, final coordination\n",
    "        else:\n",
    "            line[0:3] = np.array(line[0:3]).astype(int)   # Event #, time step, duration\n",
    "    \n",
    "line_cont = 9    # Line number continuation\n",
    "while line_cont < len(ls_lines)-2:\n",
    "    \n",
    "    ev_atom = []    # List of atom attributes for each event\n",
    "    for line_index, line in enumerate(ls_lines):\n",
    "        if line and (line_index > line_cont):\n",
    "            if len(line) == 6:\n",
    "                ev_atom.append(line)\n",
    "            else:\n",
    "                ev_no = line[0]\n",
    "                tf = line[1]\n",
    "                dt = line[2]\n",
    "                line_cont = line_index\n",
    "                break\n",
    "    \n",
    "    out = '%s%i\\n%s%i%s\\n%s%i%s\\n'%('Event #', ev_no, 'tf = ', tf, ' ps', 'dt = ', dt, ' ps')\n",
    "    out_cls.write(out)\n",
    "    \n",
    "    N = len(ev_atom)    # Number of active atoms\n",
    "    index_z = []    # Line indices of z-active atoms\n",
    "    index_aux = []    # Line indices of auxiliary atoms\n",
    "    \n",
    "    for atom_index, atom in enumerate(ev_atom):\n",
    "        if abs(atom[2] - atom[3]) > 0.50*dn_avg:\n",
    "            index_z.append(atom_index)\n",
    "        else:\n",
    "            index_aux.append(atom_index)\n",
    "    \n",
    "    N_z = len(index_z)    # Number of z-active atoms\n",
    "    N_aux = len(index_aux)    # Number of auxiliary atoms\n",
    "    \n",
    "    site = []    # Initial & final site & type matrix\n",
    "    for atom in ev_atom:\n",
    "        \n",
    "        ID = atom[0]\n",
    "        Type = atom[1]\n",
    "        n = [atom[2], atom[3]]\n",
    "        CN = [atom[4], atom[5]]\n",
    "        atom_site = [[],[]]\n",
    "        \n",
    "        for i in range(0,2):\n",
    "            \n",
    "            if n[i] < z_layer[N_slab-1]:\n",
    "                atom_site[i] = Lv_sub    # Subsurface\n",
    "            \n",
    "            elif n[i] == z_layer[N_slab-1]:\n",
    "                \n",
    "                if 11 <= CN[i] <= 12:\n",
    "                    atom_site[i] = Lv0_2    # Subdeposit\n",
    "                \n",
    "                elif 8 <= CN[i] <= 10:\n",
    "                    atom_site[i] = Lv0_1    # Surface intralayer\n",
    "            \n",
    "            elif n[i] > z_layer[N_slab-1]:\n",
    "                \n",
    "                if 8 <= CN[i] <= 12:\n",
    "                    atom_site[i] = Lv1_3    # Deposit intralayer\n",
    "                \n",
    "                elif 4 <= CN[i] <= 7:\n",
    "                    atom_site[i] = Lv1_2    # Edge/kink/vacancy\n",
    "            \n",
    "                elif 1 <= CN[i] <= 3:\n",
    "                    atom_site[i] = Lv1_1    # Isolated adatom\n",
    "        \n",
    "        site.append(atom_site)\n",
    "        out = '%i\\t%i\\t%s %s %s\\n'%(ID, Type, atom_site[0], '==>', atom_site[1])\n",
    "        out_cls.write(out)\n",
    "\n",
    "    # More than 5 active atoms = Indeterminate\n",
    "    if N > 5:\n",
    "        out = 'Visual inspection needed: More than 5 active atoms. Most likely conjoined events or intralayer shift.\\n'\n",
    "        out_cls.write(out)\n",
    "\n",
    "    # 1 z-active atom = Popout vs. insertion vs. step motion\n",
    "    elif N_z == 1:\n",
    "        \n",
    "        Type = ev_atom[index_z[0]][1]\n",
    "        \n",
    "        if Lv_sub in site[index_z[0]]:\n",
    "            out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif ((site[index_z[0]][0] == Lv0_1) or (site[index_z[0]][0] == Lv0_2)) and ((site[index_z[0]][1] == Lv1_1) or (site[index_z[0]][1] == Lv1_2) or (site[index_z[0]][1] == Lv1_3)):\n",
    "            out = 'Class = Popout\\nAtom type = ' + str(Type) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif ((site[index_z[0]][0] == Lv1_1) or (site[index_z[0]][0] == Lv1_2) or (site[index_z[0]][0] == Lv1_3)) and ((site[index_z[0]][1] == Lv0_1) or (site[index_z[0]][1] == Lv0_2)):\n",
    "            out = 'Class = Vacancy insertion\\nAtom type = ' + str(Type) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif ((site[index_z[0]][0] == Lv1_1) or (site[index_z[0]][0] == Lv1_2) or (site[index_z[0]][0] == Lv1_3)) and ((site[index_z[0]][1] == Lv1_1) or (site[index_z[0]][1] == Lv1_2) or (site[index_z[0]][1] == Lv1_3)):\n",
    "            \n",
    "            dn = ev_atom[index_z[0]][3] - ev_atom[index_z[0]][2]\n",
    "            \n",
    "            if N_aux != 0:\n",
    "                if dn > 0:\n",
    "                    out = 'Class = Exchange ascent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "                else:\n",
    "                    out = 'Class = Exchange descent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "            \n",
    "            else:\n",
    "                if dn > 0:\n",
    "                    out = 'Class = Hopping ascent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "                else:\n",
    "                    out = 'Class = Hopping descent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "    \n",
    "    # 2 z-active atoms = Popout vs. insertion vs. adatom exchange\n",
    "    elif N_z == 2:\n",
    "        \n",
    "        ID1 = ev_atom[index_z[0]][0]\n",
    "        ID2 = ev_atom[index_z[1]][0]\n",
    "        Type1 = ev_atom[index_z[0]][1]\n",
    "        Type2 = ev_atom[index_z[1]][1]\n",
    "        dn1 = ev_atom[index_z[0]][3] - ev_atom[index_z[0]][2]\n",
    "        dn2 = ev_atom[index_z[1]][3] - ev_atom[index_z[1]][2]\n",
    "\n",
    "        if Lv_sub in site[index_z[0]]:\n",
    "            out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type1) + '\\n'\n",
    "            out_cls.write(out)\n",
    "            \n",
    "        elif Lv_sub in site[index_z[1]]:\n",
    "            out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type2) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif (dn1 > 0) and (dn2 > 0):\n",
    "            \n",
    "            if ((site[index_z[0]][0] == Lv0_1) or (site[index_z[0]][0] == Lv0_2)):\n",
    "                out = 'Class = Popout\\nAtom type = ' + str(Type1) + '\\n'\n",
    "                out_cls.write(out)\n",
    "                \n",
    "            elif ((site[index_z[1]][0] == Lv0_1) or (site[index_z[1]][0] == Lv0_2)):\n",
    "                out = 'Class = Popout\\nAtom type = ' + str(Type2) + '\\n'\n",
    "                out_cls.write(out)\n",
    "            \n",
    "            else:\n",
    "                out = 'Visual inspection needed: Most likely simultaneous ascent.\\n'\n",
    "                out_cls.write(out)\n",
    "        \n",
    "        elif (dn1 < 0) and (dn2 < 0):\n",
    "            \n",
    "            if ((site[index_z[0]][1] == Lv0_1) or (site[index_z[0]][1] == Lv0_2)):\n",
    "                out = 'Class = Vacancy insertion\\nAtom type = ' + str(Type1) + '\\n'\n",
    "                out_cls.write(out)\n",
    "\n",
    "            elif ((site[index_z[1]][1] == Lv0_1) or (site[index_z[1]][1] == Lv0_2)):\n",
    "                out = 'Class = Vacancy insertion\\nAtom type = ' + str(Type2) + '\\n'\n",
    "                out_cls.write(out)\n",
    "            \n",
    "            else:\n",
    "                out = 'Visual inspection needed: Most likely simultaneous descent.\\n'\n",
    "                out_cls.write(out)\n",
    "        \n",
    "        else:\n",
    "            \n",
    "            if N_aux != 0:\n",
    "                out = 'Class = Indirect exchange\\nAtom type = ' + str(Type1) + '+' + str(Type2) + '\\n'\n",
    "                out_cls.write(out)\n",
    "            \n",
    "            else:    # Distance check\n",
    "                \n",
    "                for atom in range(0,N_dump):\n",
    "                    if xyz1.finaldict[0][1]['id'][atom] == ID1:\n",
    "                        index1 = atom\n",
    "                    elif xyz1.finaldict[0][1]['id'][atom] == ID2:\n",
    "                        index2 = atom\n",
    "                \n",
    "                x1 = xyz1.finaldict[step[tf]][1]['xu'][index1]\n",
    "                y1 = xyz1.finaldict[step[tf]][1]['yu'][index1]\n",
    "                z1 = xyz1.finaldict[step[tf]][1]['z'][index1]\n",
    "                r1 = np.array([x1, y1, z1])\n",
    "                \n",
    "                x2 = xyz1.finaldict[step[tf]][1]['xu'][index2]\n",
    "                y2 = xyz1.finaldict[step[tf]][1]['yu'][index2]\n",
    "                z2 = xyz1.finaldict[step[tf]][1]['z'][index2]\n",
    "                r2 = np.array([x2, y2, z2])\n",
    "                \n",
    "                dr = np.subtract(r1,r2)\n",
    "                \n",
    "                if np.linalg.norm(dr) > 1.10*dr_avg:\n",
    "                    out = 'Class = Indirect exchange\\nAtom type = ' + str(Type1) + '+' + str(Type2) + '\\n' + 'Visual inspection advised: Missing auxiliary atom.\\n'\n",
    "                    out_cls.write(out)\n",
    "                else:\n",
    "                    out = 'Class = Direct exchange\\nAtom type = ' + str(Type1) + '+' + str(Type2) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "    \n",
    "    # More than 2 z-active atoms = Interlayer exchange vs. indeterminate\n",
    "    elif N_z > 2:\n",
    "        \n",
    "        N_sub = 0\n",
    "        for index, i in enumerate(index_z):\n",
    "            \n",
    "            if Lv_sub in site[i]:\n",
    "                N_sub += 1\n",
    "                index_last = index\n",
    "            \n",
    "            if index == len(index_z)-1:\n",
    "                \n",
    "                if N_sub == 0:\n",
    "                    out = 'Visual inspection needed: More than two z-active atoms. Most likely conjoined events, or extra fluctuating atoms.\\n'\n",
    "                    out_cls.write(out)\n",
    "                \n",
    "                elif N_sub == 1:                    \n",
    "                    Type = ev_atom[index_last][1]\n",
    "                    out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "\n",
    "                else:\n",
    "                    out = 'Class = Interlayer transfer\\nAtom type = Multi\\n'\n",
    "                    out_cls.write(out)\n",
    "\n",
    "    out_cls.write('\\n')\n",
    "    \n",
    "out_cls.close()\n",
    "ls.close()\n",
    "\n",
    "print('Event classification done > event_class.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_cls = pwd + '/event_class.txt'    # List of event classes\n",
    "cls = open(file_cls,'r')    # Load list of event classes\n",
    "cls_lines = cls.readlines()\n",
    "\n",
    "file_ind = pwd + '/event_indeterminate.txt'    # Indeterminate events\n",
    "out_ind = open(file_ind,'w')\n",
    "\n",
    "file_stat = pwd + '/event_stat.txt'    # Event statistics\n",
    "out_stat = open(file_stat,'w')\n",
    "\n",
    "N_event = 0    # Number of events\n",
    "N_indeterminate = 0    # Number of indeterminate events\n",
    "for line_index, line in enumerate(cls_lines):\n",
    "    if 'Event' in line:\n",
    "        N_event += 1\n",
    "    elif 'Visual inspection needed:' in line:\n",
    "        N_indeterminate += 1\n",
    "        \n",
    "        end = line_index\n",
    "        for i in reversed(range(end)):\n",
    "            if 'Event' in cls_lines[i]:\n",
    "                start = i\n",
    "                break\n",
    "        \n",
    "        for i in range(start,end+1):\n",
    "            out_ind.write(cls_lines[i])\n",
    "        out_ind.write('\\n')\n",
    "        \n",
    "indeterminacy = 100*N_indeterminate/N_event\n",
    "header1 = 'Indeterminacy = %.2f' % (indeterminacy) + '% (' + str(N_indeterminate) + ' events)\\n'\n",
    "out_stat.write(header1)\n",
    "\n",
    "Class = ['Popout', 'Vacancy insertion', 'Exchange descent', 'Exchange ascent', 'Hopping descent', 'Hopping ascent', 'Indirect exchange', 'Direct exchange', 'Interlayer transfer']\n",
    "if N_type == 1:\n",
    "    hist = np.zeros((9,), dtype=int)\n",
    "else:\n",
    "    hist = np.zeros((21,), dtype=int)\n",
    "\n",
    "for line_index, line in enumerate(cls_lines):    # Count number of events in each class\n",
    "    if 'Class' in line:\n",
    "        \n",
    "        for C_index, C in enumerate(Class):\n",
    "            if C in line:\n",
    "                \n",
    "                for Type_index, Type in enumerate(range(1,N_type+1)):\n",
    "                    \n",
    "                    if C_index <= 5:\n",
    "                        if str(Type) in cls_lines[line_index+1]:\n",
    "                            hist[2*C_index + Type_index] += 1\n",
    "                            \n",
    "                    elif 6 <= C_index <= 7:\n",
    "                        \n",
    "                        Atom_type = str(Type) + '+' + str(Type)\n",
    "                        if Atom_type in cls_lines[line_index+1]:\n",
    "                            hist[2*C_index + Type_index + (C_index-6)] += 1\n",
    "                        \n",
    "                        if Type_index == 1:\n",
    "                            Atom_type1 = str(Type-1) + '+' + str(Type)\n",
    "                            Atom_type2 = str(Type) + '+' + str(Type-1)\n",
    "                            if (Atom_type1 in cls_lines[line_index+1]) or (Atom_type2 in cls_lines[line_index+1]):\n",
    "                                hist[2*C_index + Type_index + (C_index-6) + 1] += 1\n",
    "                    \n",
    "                    else:\n",
    "                        \n",
    "                        if str(Type) in cls_lines[line_index+1]:\n",
    "                            hist[2*C_index + Type_index + 2] += 1\n",
    "                        \n",
    "                        if Type_index == 1:\n",
    "                            if 'Multi' in cls_lines[line_index+1]:\n",
    "                                hist[-1] += 1\n",
    "                    \n",
    "total = np.sum(hist)\n",
    "header2 = str(total) + ' out of ' + str(N_event) + ' events classified.\\n\\n'\n",
    "out_stat.write(header2)\n",
    "\n",
    "header3 = ['Event class', 'Atom type', '# of occurrence', '% frequency']\n",
    "types = ['%%%is', '%%%is', '%%%ii', '%%%i.2f']\n",
    "colwidth = max(max(len(h) for h in header3), max(len(c) for c in Class))\n",
    "out_stat.write((('\\t'.join([\"%%%is\" % colwidth] * len(header3))) % tuple(header3)) + \"\\n\")\n",
    "\n",
    "for C_index, C in enumerate(Class):\n",
    "    for Type_index, Type in enumerate(range(1,N_type+1)):\n",
    "        \n",
    "        if C_index <= 5:\n",
    "            number = hist[2*C_index + Type_index]\n",
    "            freq = float(100*number/total)\n",
    "            out = ('\\t'.join(t % colwidth for t in types)) % (C, str(Type), number, freq)\n",
    "            out += '\\n'\n",
    "            out_stat.write(out)\n",
    "            \n",
    "        elif 6 <= C_index <= 7:\n",
    "            \n",
    "            Atom_type = str(Type) + '+' + str(Type)\n",
    "            number = hist[2*C_index + Type_index + (C_index-6)]\n",
    "            freq = float(100*number/total)\n",
    "            out = ('\\t'.join(t % colwidth for t in types)) % (C, Atom_type, number, freq)\n",
    "            out += '\\n'\n",
    "            out_stat.write(out)\n",
    "\n",
    "            if Type_index == 1:\n",
    "                Atom_type = str(Type-1) + '+' + str(Type)\n",
    "                number = hist[2*C_index + Type_index + (C_index-6) + 1]\n",
    "                freq = float(100*number/total)\n",
    "                out = ('\\t'.join(t % colwidth for t in types)) % (C, Atom_type, number, freq)\n",
    "                out += '\\n'\n",
    "                out_stat.write(out)\n",
    "        \n",
    "        else:\n",
    "            \n",
    "            number = hist[2*C_index + Type_index + 2]\n",
    "            freq = float(100*number/total)\n",
    "            out = ('\\t'.join(t % colwidth for t in types)) % (C, str(Type), number, freq)\n",
    "            out += '\\n'\n",
    "            out_stat.write(out)\n",
    "\n",
    "            if Type_index == 1:\n",
    "                Atom_type = 'Multi'\n",
    "                number = hist[-1]\n",
    "                freq = float(100*number/total)\n",
    "                out = ('\\t'.join(t % colwidth for t in types)) % (C, Atom_type, number, freq)\n",
    "                out += '\\n'\n",
    "                out_stat.write(out)\n",
    "\n",
    "out_stat.close()\n",
    "out_ind.close()\n",
    "cls.close()\n",
    "\n",
    "print('Event statistics done > event_stat.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LAMMPS NEB preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timestep = float(input('Enter the timestep used for the simulation, in ps : '))\n",
    "dmp_step = 0.1    # Dump frequency in ps\n",
    "d_step = int(dmp_step/timestep)    # Number of steps in a 0.1 ps interval\n",
    "\n",
    "# Header for LAMMPS DUMP file of initial trajectory\n",
    "l_step = 'ITEM: TIMESTEP'\n",
    "l_atom = 'ITEM: NUMBER OF ATOMS'\n",
    "l_box = 'ITEM: BOX BOUNDS xy xz yz pp pp ff'\n",
    "l_xyz = 'ITEM: ATOMS id type x y z'\n",
    "\n",
    "print('Preparing LAMMPS NEB images in ./neb/ev_***/ for each event...')\n",
    "\n",
    "ls = open(file_ls,'r')    # Load list of events\n",
    "ls_lines = ls.readlines()\n",
    "ls_lines = [line.split() for line in ls_lines]\n",
    "\n",
    "log = open(file_eq,'r')\n",
    "log_lines = log.readlines()\n",
    "\n",
    "for line_index, line in enumerate(ls_lines):\n",
    "    if line:\n",
    "        if (line_index > 9) and (len(line) == 3):\n",
    "        \n",
    "            ev_no = line[0]\n",
    "            tf = line[1]\n",
    "            dt = line[2]\n",
    "\n",
    "            tf_step = int(tf/timestep)    # Number of 0.1 ps steps\n",
    "            ti_step = tf_step - int(dt/timestep)\n",
    "\n",
    "            ev_dir = pwd + '/ev_' + str(ev_no)    # Create NEB directory for each event\n",
    "            os.mkdir(ev_dir)\n",
    "\n",
    "            file_dmp = ev_dir + '/initial.dmp'    # LAMMPS DUMP file for visualizing initial trajectory\n",
    "            out_dmp = open(file_dmp,'w')\n",
    "\n",
    "            step = range(ti_step, tf_step+1, d_step)    # Event interval\n",
    "            for s_index, s in enumerate(step):\n",
    "\n",
    "                if s_index == 0:\n",
    "\n",
    "                    file_img0 = ev_dir + '/img0.dat'    # LAMMPS DATA file for initial structure (image 0)\n",
    "                    out_img0 = open(file_img0,'w')\n",
    "\n",
    "                    for l_index, l in enumerate(log_lines):\n",
    "\n",
    "                        if l_index <= head:\n",
    "                            out_img0.write(l)\n",
    "\n",
    "                    for atom in range(0,N_dump):\n",
    "\n",
    "                        ID = xyz2.finaldict[s][1]['id'][atom]\n",
    "                        Type = xyz2.finaldict[s][1]['type'][atom]\n",
    "                        x = xyz2.finaldict[s][1]['x'][atom]\n",
    "                        y = xyz2.finaldict[s][1]['y'][atom]\n",
    "                        z = xyz2.finaldict[s][1]['z'][atom]\n",
    "\n",
    "                        out = '%i\\t%i\\t%f\\t%f\\t%f\\n'%(ID, Type, x, y, z)\n",
    "                        out_img0.write(out)\n",
    "\n",
    "                    out_img0.close()\n",
    "\n",
    "                header = '%s\\n%i\\n%s\\n%i\\n%s\\n%f %f %f\\n%f %f %f\\n%f %f %f\\n%s\\n'%(l_step, s, l_atom, N_dump, l_box, xlo_ext, xhi_ext, xy, ylo, yhi, xz, zlo, zhi, yz, l_xyz)\n",
    "                out_dmp.write(header)\n",
    "\n",
    "                file_trj = ev_dir + '/coords.initial.' + str(s_index)    # Images 1 to N\n",
    "                out_trj = open(file_trj,'w')\n",
    "\n",
    "                header = '# ID \\t x \\t y \\t z \\n'\n",
    "                N_line = '%i\\n'%(N_dump)\n",
    "                out_trj.write(header)\n",
    "                out_trj.write(N_line)\n",
    "\n",
    "                for atom in range(0,N_dump):\n",
    "\n",
    "                    ID = xyz2.finaldict[s][1]['id'][atom]\n",
    "                    Type = xyz2.finaldict[s][1]['type'][atom]\n",
    "                    x = xyz2.finaldict[s][1]['x'][atom]\n",
    "                    y = xyz2.finaldict[s][1]['y'][atom]\n",
    "                    z = xyz2.finaldict[s][1]['z'][atom]\n",
    "\n",
    "                    out = '%i %i %f %f %f\\n'%(ID, Type, x, y, z)\n",
    "                    out_dmp.write(out)\n",
    "\n",
    "                    out = '%i\\t%f\\t%f\\t%f\\n'%(ID, x, y, z)\n",
    "                    out_trj.write(out)\n",
    "\n",
    "                out_trj.close()\n",
    "\n",
    "            out_dmp.close()\n",
    "\n",
    "ls.close()\n",
    "log.close()\n",
    "\n",
    "print('LAMMPS NEB images generated > ./ev_***/')\n",
    "print('Use initial.dmp to visualize the initial trajectory.')\n",
    "print('Use img0.dat as the initial structure.')\n",
    "print('Use coords.initial.** as the NEB images.\\n')\n",
    "\n",
    "print('Once the NEB calculations converge, NEB2DFT.ipynb or NEB2ASE.ipynb can be used to prepare DFT or ASE optimization, respectively.')"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}