gamma_surfaces.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from gamma_surfaces import get_gamma_surface, show_master_gamma_surface, show_gamma_surface_fit\n",
    "from utilities import make_structure, get_lammps_parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pot_base_dir = Path('.').resolve()\n",
    "common_lammps_params = get_lammps_parameters(base_path=pot_base_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For conversion from eV / Ang^2 --> J / m^2\n",
    "UNIT_CONV = 16.02176565"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: if examining any of the twist GB gamma surfaces, the energy landscape is shown on orthogonal axes (when in fact the axes are non-orthogonal)\n",
    "gs = get_gamma_surface('s7-tlA', 'dft')\n",
    "show_master_gamma_surface(gs, 'grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs = get_gamma_surface('s13-tlA', 'eam')\n",
    "show_master_gamma_surface(gs, 'grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualise gamma surface data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below uses JSON data files that are stored within this repository to visualise the fitted gamma surfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "structure_code = 's7-tlA'\n",
    "method = 'dft'\n",
    "\n",
    "gs = get_gamma_surface(structure_code, method)\n",
    "show_master_gamma_surface(gs, data_name='grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "structure_code = 's13-tlA'\n",
    "method = 'eam'\n",
    "\n",
    "gs = get_gamma_surface(structure_code, method)\n",
    "show_master_gamma_surface(gs, data_name='grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the translation with the minimum fitted grain boundary energy:\n",
    "min_shift = gs.get_minimum_fitted_shift('grain_boundary_energy')[0].shift\n",
    "min_shift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show quadratic fitting at a given shift:\n",
    "show_gamma_surface_fit(gs, shift=min_shift, data_name='grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show quadratic fitting at a given shift:\n",
    "show_gamma_surface_fit(gs, shift=[0.25, 0.5], data_name='grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the fitted minimum expansion and grain boundary energy:\n",
    "gs.fitted_data['grain_boundary_energy']['minimum'][\n",
    "    np.argmin(gs.fitted_data['grain_boundary_energy']['minimum'][:, 1]), :\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true,
    "toc-nb-collapsed": true
   },
   "source": [
    "# Method of gamma surface generation and analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section breaks down the steps used to:\n",
    "\n",
    "- generate gamma surface input files (for LAMMPS)\n",
    "- run the simulations\n",
    "- collate the results\n",
    "- perform quadratic fitting\n",
    "- visualise the master gamma surface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Note**: you cannot run simulations on Binder!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from subprocess import run, PIPE\n",
    "import numpy as np\n",
    "\n",
    "from atomistic.bicrystal import GammaSurface\n",
    "from lammps_parse import write_lammps_inputs, read_lammps_output\n",
    "from plotly import graph_objects\n",
    "\n",
    "# Change this variable to the directory where LAMMPS gamma surface simulation input files should be generated:\n",
    "LAMMPS_SIMS_DIR = Path('/path/to/simulation/directory')\n",
    "\n",
    "# Change this to the LAMMPS executable name (note, you cannot run simulations on Binder!):\n",
    "LAMMPS_EXECUTABLE = 'lmp_serial'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generate a bicrystal whose gamma surface is to be explored**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "structure_key = 's7-tlA'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bicrystal = make_structure(f'{structure_key}-gb', method='eam', configuration='sized')\n",
    "bicrystal.show(\n",
    "    include={'points': ['atoms']},\n",
    "    layout_args={'height': 800},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generate a corresponding bulk model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bulk = make_structure(f'{structure_key}-b', method='eam', configuration='sized')\n",
    "bulk.show(\n",
    "    include={'points': ['atoms']},\n",
    "    layout_args={'height': 600},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generate a gamma surface from a grid**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_surface = GammaSurface.from_grid(bicrystal, grid=[12, 18], expansions=[-1, 0, 1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Iterate over all coordinates in the gamma surface to generate simulation input files**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pot_base_dir = Path('.').resolve()\n",
    "GB_input_paths = []\n",
    "\n",
    "common_lammps_params = get_lammps_parameters(base_path=pot_base_dir)\n",
    "\n",
    "for i in gamma_surface.all_coordinates():\n",
    "    \n",
    "    # Make a directory for this sim:\n",
    "    sim_path = LAMMPS_SIMS_DIR.joinpath(f'{structure_key}-gb', i.coordinate_fmt)\n",
    "    sim_path.mkdir(parents=True)\n",
    "    \n",
    "    # Write simulation inputs:\n",
    "    lammps_params_gb_i = {\n",
    "        'supercell': i.structure.supercell,\n",
    "        'atom_sites': i.structure.atoms.coords,\n",
    "        'species': i.structure.atoms.labels['species'].unique_values,\n",
    "        'species_idx': i.structure.atoms.labels['species'].values_idx,\n",
    "        'dir_path': sim_path,\n",
    "        **common_lammps_params,\n",
    "    }\n",
    "    GB_input_paths.append(write_lammps_inputs(**lammps_params_gb_i))\n",
    "    \n",
    "# Make a directory for the bulk sim:\n",
    "sim_path = LAMMPS_SIMS_DIR.joinpath(f'{structure_key}-b')\n",
    "sim_path.mkdir()\n",
    "\n",
    "# Write bulk simulation inputs:\n",
    "lammps_params_bulk = {\n",
    "    'supercell': bulk.supercell,\n",
    "    'atom_sites': bulk.atoms.coords,\n",
    "    'species': bulk.atoms.labels['species'].unique_values,\n",
    "    'species_idx': bulk.atoms.labels['species'].values_idx,\n",
    "    'dir_path': sim_path,\n",
    "    **common_lammps_params,\n",
    "}\n",
    "bulk_input_path = write_lammps_inputs(**lammps_params_bulk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Run simulations**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i in GB_input_paths + [bulk_input_path]:\n",
    "    cmd = '{} < {}'.format(LAMMPS_EXECUTABLE, i.name)\n",
    "    proc = run(cmd, shell=True, cwd=i.parent , stdout=PIPE, stderr=PIPE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Collate simulation outputs**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_gamma_surface_params = {\n",
    "    'shifts': [],\n",
    "    'expansions': [],\n",
    "    'data': {\n",
    "        'energy': [],\n",
    "        'grain_boundary_energy': [],\n",
    "    },\n",
    "    'metadata': {},\n",
    "}\n",
    "# First get the bulk sim data:\n",
    "lammps_out_bulk = read_lammps_output(dir_path=bulk_input_path.parent)\n",
    "E_tot_bulk = lammps_out_bulk['final_energy'][-1]\n",
    "simulated_gamma_surface_params['metadata'].update({\n",
    "    'E_tot_bulk': E_tot_bulk,\n",
    "    'grain_boundary_area': bicrystal.boundary_area,\n",
    "    'num_atoms_bulk': bulk.num_atoms,\n",
    "    'num_atoms_grain_boundary': bicrystal.num_atoms,\n",
    "})\n",
    "\n",
    "# Now iterate over GB sims:\n",
    "for i in GB_input_paths:\n",
    "        \n",
    "    shift_str, exp_str = i.parent.name.split('__')\n",
    "    shift = []\n",
    "    for j in shift_str.split('_'):\n",
    "        num, denom = j.split('.')\n",
    "        shift.append(int(num) / int(denom))\n",
    "        \n",
    "    simulated_gamma_surface_params['shifts'].append(shift)\n",
    "    simulated_gamma_surface_params['expansions'].append(float(exp_str))\n",
    "    \n",
    "    lammps_out_GB_i = read_lammps_output(dir_path=i.parent)    \n",
    "    E_tot_GB_i = lammps_out_GB_i['final_energy'][-1]\n",
    "    simulated_gamma_surface_params['data']['energy'].append(E_tot_GB_i)    \n",
    "    \n",
    "    E_GB_i = (1 / (2 * bicrystal.boundary_area)) * (\n",
    "        E_tot_GB_i - (bicrystal.num_atoms / bulk.num_atoms) * E_tot_bulk\n",
    "    ) * UNIT_CONV\n",
    "    simulated_gamma_surface_params['data']['grain_boundary_energy'].append(E_GB_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_gamma_surface = GammaSurface(bicrystal, **simulated_gamma_surface_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot a slice of the gamma surface**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_label = 'grain_boundary_energy'\n",
    "plot_dat = simulated_gamma_surface.get_surface_plot_data(data_label, expansion=0, xy_as_grid=False)\n",
    "grid_dat = simulated_gamma_surface.get_xy_plot_data()\n",
    "graph_objects.FigureWidget(\n",
    "    data=[\n",
    "        {\n",
    "            'type': 'contour',\n",
    "            'colorscale': 'viridis',\n",
    "            'colorbar': {\n",
    "                'title': data_label,\n",
    "            },\n",
    "            **plot_dat,            \n",
    "        },\n",
    "        {\n",
    "            **grid_dat,\n",
    "            'mode': 'markers',\n",
    "            'marker': {\n",
    "                'size': 2,\n",
    "            },\n",
    "        },\n",
    "    ],\n",
    "    layout={\n",
    "        'xaxis': {\n",
    "            'scaleanchor': 'y',            \n",
    "        },        \n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Fit at each shift to find the master gamma surface, and plot**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_gamma_surface.add_fit(data_label, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "master_plot_data = simulated_gamma_surface.get_fitted_surface_plot_data(data_label, xy_as_grid=False)\n",
    "graph_objects.FigureWidget(\n",
    "    data=[\n",
    "        {\n",
    "            'type': 'contour',\n",
    "            'colorscale': 'viridis',\n",
    "            **master_plot_data,            \n",
    "        },\n",
    "        {\n",
    "            **grid_dat,\n",
    "            'mode': 'markers',\n",
    "            'marker': {\n",
    "                'size': 2,\n",
    "            },\n",
    "        },        \n",
    "    ],\n",
    "    layout={\n",
    "        'xaxis': {\n",
    "            'scaleanchor': 'y',\n",
    "        }\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plot the fit at a given shift**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_plot_dat = simulated_gamma_surface.get_fit_plot_data(data_label, [0, 0])\n",
    "fig = graph_objects.FigureWidget(\n",
    "    data=[\n",
    "        {\n",
    "            **fit_plot_dat['fitted_data'],\n",
    "            'name': 'Fit',\n",
    "        },\n",
    "        {\n",
    "            **fit_plot_dat['data'],\n",
    "            'name': data_label,\n",
    "        },\n",
    "        {\n",
    "            **fit_plot_dat['minimum'],\n",
    "            'name': 'Fit min.',\n",
    "        },\n",
    "    ],\n",
    "    layout={\n",
    "        'xaxis': {\n",
    "            'title': 'Expansion',\n",
    "        },\n",
    "        'yaxis': {\n",
    "            'title': data_label,\n",
    "        },\n",
    "        'width': 500,\n",
    "        'height': 400,\n",
    "    }\n",
    ")\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true,
    "toc-nb-collapsed": true
   },
   "source": [
    "# Using wrapper functions in `gamma_surfaces.py`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above workflow is wrapped up into a few wrapper functions in the python file `gamma_surfaces.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gamma_surfaces import compute_master_gamma_surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "structure_code = 's7-tlA'\n",
    "gamma_surface = compute_master_gamma_surface(structure_code, LAMMPS_SIMS_DIR) # 5-10 mins run time for the STGBs (much longer for the twists)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_master_gamma_surface(gamma_surface, data_name='grain_boundary_energy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LAMMPS simulations to find grain boundary energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the minimum-energy configurations (shift and expansion) from the EAM γ-surfaces to do full-atom relaxations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "common_lammps_params = copy.deepcopy(common_lammps_params)\n",
    "common_lammps_params['atom_constraints']['fix_xy_idx'] = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "structure_codes = [\n",
    "    's7-tlA',\n",
    "    's13-tlA',\n",
    "    's19-tlA',\n",
    "    's31-tlA',\n",
    "    's7-tw',\n",
    "    's13-tw',\n",
    "    's19-tw',\n",
    "    's7-tlB',    \n",
    "]\n",
    "E_GB_EAM = {}\n",
    "for structure_code in structure_codes:    \n",
    "    bicrystal = make_structure(\n",
    "        f'{structure_code}-gb',\n",
    "        method='eam',\n",
    "        configuration='minimum_energy',\n",
    "    )\n",
    "    bulk = make_structure(\n",
    "        f'{structure_code}-b',\n",
    "        method='eam',\n",
    "        configuration='minimum_energy'\n",
    "    )\n",
    "    \n",
    "    sim_path_GB = LAMMPS_SIMS_DIR.joinpath('full_relaxations', f'{structure_code}-gb')\n",
    "    sim_path_GB.mkdir(parents=True)\n",
    "    \n",
    "    # Write simulation inputs:\n",
    "    lammps_params_gb = {\n",
    "        'supercell': bicrystal.supercell,\n",
    "        'atom_sites': bicrystal.atoms.coords,\n",
    "        'species': bicrystal.atoms.labels['species'].unique_values,\n",
    "        'species_idx': bicrystal.atoms.labels['species'].values_idx,\n",
    "        'dir_path': sim_path_GB,\n",
    "        **common_lammps_params,\n",
    "    }\n",
    "    input_path_GB = write_lammps_inputs(**lammps_params_gb)    \n",
    "    \n",
    "    sim_path_bulk = LAMMPS_SIMS_DIR.joinpath('full_relaxations', f'{structure_code}-b', )\n",
    "    sim_path_bulk.mkdir(parents=True)\n",
    "    \n",
    "    # Write simulation inputs:\n",
    "    lammps_params_bulk = {\n",
    "        'supercell': bulk.supercell,\n",
    "        'atom_sites': bulk.atoms.coords,\n",
    "        'species': bulk.atoms.labels['species'].unique_values,\n",
    "        'species_idx': bulk.atoms.labels['species'].values_idx,\n",
    "        'dir_path': sim_path_bulk,\n",
    "        **common_lammps_params,\n",
    "    }\n",
    "    input_path_bulk = write_lammps_inputs(**lammps_params_bulk)      \n",
    "    \n",
    "    for i in [input_path_GB, input_path_bulk]:\n",
    "        cmd = '{} < {}'.format(LAMMPS_EXECUTABLE, i.name)\n",
    "        proc = run(cmd, shell=True, cwd=i.parent , stdout=PIPE, stderr=PIPE)\n",
    "    \n",
    "    lammps_out_bulk = read_lammps_output(dir_path=input_path_bulk.parent)\n",
    "    lammps_out_GB = read_lammps_output(dir_path=input_path_GB.parent)\n",
    "    \n",
    "    E_tot_bulk = lammps_out_bulk['final_energy'][-1]\n",
    "    E_tot_GB = lammps_out_GB['final_energy'][-1]\n",
    "    \n",
    "    E_GB = (1 / (2 * bicrystal.boundary_area)) * (\n",
    "        E_tot_GB - (bicrystal.num_atoms / bulk.num_atoms) * E_tot_bulk\n",
    "    ) * UNIT_CONV\n",
    "    print(f'{structure_code}: {E_GB}')\n",
    "    E_GB_EAM.update({structure_code: E_GB})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "interpreter": {
   "hash": "c775c1e6b3eb5e71267bf5e16cde787343ccf7d7636bba966b077c32ed40a4ad"
  },
  "kernelspec": {
   "display_name": "Python [conda env:zr_gb_papers_env]",
   "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.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}