1-SEI_Builder-step1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "separated-mobile",
   "metadata": {},
   "source": [
    "# SEI (Solid Electrolyte Interface) Builder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "effective-environment",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <table style=\"border-collapse:collapse; vertical-align:top\">\n",
    "    <col style=\"width:35%; min-width:150px;\"/>\n",
    "    <col style=\"width:65%; min-width:200px;\"/>\n",
    "    <tfoot>\n",
    "        <tr>\n",
    "            <td style=\"border-right:2px solid #ffffff; min-width:300px\">\n",
    "                <a href=\"https://www.big-map.eu/\">\n",
    "                    <img style=\"width:12%; min-width:25px\" src=\"img/logo-bigmap.png\" alt=\"BIG MAP site\" >\n",
    "                </a>\n",
    "                <a href=\"https://areeweb.polito.it/ricerca/small/\">\n",
    "                    <img style=\"width:50%; min-width:100px\" src=\"img/logo-small.png\" alt=\"SMALL site\" >\n",
    "                </a>\n",
    "                <a href=\"https://www.polito.it/\">\n",
    "                    <img style=\"width:30%; min-width:60px\" src=\"img/logo-polito.png\" alt=\"POLITO site\" >\n",
    "                </a>\n",
    "            </td>\n",
    "            <td style=\"border-collapse:collapse;color:#0097a7; font-size:11pt; text-align:right; min-width:200px\">\n",
    "                <b>Autors:</b><br>            \n",
    "                Paolo De Angelis (<a href=\"mailto:paolodeangelis@polito.it\">paolodeangelis@polito.it</a>)<br>  \n",
    "                Roberta Cappabianca (<a href=\"mailto:roberta.cappabianca@polito.it\">roberta.cappabianca@polito.it</a>)<br>  \n",
    "                Eliodoro Chiavazzo (<a href=\"mailto:eliodoro.chiavazzo@polito.it\">eliodoro.chiavazzo@polito.it</a>)<br>  \n",
    "                Pietro Asinari (<a href=\"mailto:pietro.asinari@polito.it\">pietro.asinari@polito.it</a>)\n",
    "            </td>\n",
    "        </tr>\n",
    "    </tfoot>\n",
    "    </table>\n",
    "</div>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "wicked-wrestling",
   "metadata": {},
   "source": [
    "**Base Idea**\n",
    "\n",
    "Reconstruction of SEI by randomly positioning of monocrystal grains in mass proportion from the previous step and from a size distribution (from experiments?).\n",
    "\n",
    "**Assumptions**\n",
    "* SEI is made by packing monocrystal grains.\n",
    "* The grain size follows a given size distribution.\n",
    "* Each grain shape can be obtained by cutting the supercell with different crystal planes.\n",
    "* By annealing the system, we will allow the atoms on the grain surfaces to form the \"grain bonds.\"\n",
    "\n",
    "![Sei builder sketch](img/sei_builder_sketch.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "complicated-prison",
   "metadata": {},
   "source": [
    "## Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alternate-thesis",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "\n",
    "import matplotlib.gridspec as gridspec\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "from ase import Atoms\n",
    "from ase.io import write\n",
    "from ase.lattice.cubic import DiamondFactory\n",
    "from ase.lattice.hexagonal import HexagonalFactory\n",
    "from ase.visualize import view\n",
    "from ase.visualize.plot import plot_atoms\n",
    "from IPython.display import Image, display\n",
    "from numpy.random import PCG64, Generator\n",
    "from scipy.spatial.transform import Rotation\n",
    "\n",
    "from seibuilder import (\n",
    "    get_bulk_atoms,\n",
    "    get_gcd_pedices,\n",
    "    get_stable_crystal,\n",
    "    random_sei_grains,\n",
    ")\n",
    "from seibuilder.packmol import PackmolStructure, packmol\n",
    "from seibuilder.utils import makedir\n",
    "\n",
    "plt.style.use(\"default\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7b9555d",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\">\n",
    "<b>⚠️ Danger:</b> if the previus <i>cell code</i> where we import all the useful libraries you get the \n",
    "following message\n",
    "    \n",
    "<code>.../SEI_builder/seibuilder/mpinterfaces/\\_\\_init\\_\\_.py:35: UserWarning: mpint_config.yaml file not configured.</code>\n",
    "\n",
    "You have to get and store the Material Project API Key. Read the relative section in the README file or \n",
    "the documentation.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0cc47c3",
   "metadata": {},
   "source": [
    "## Settings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "953b0b7b",
   "metadata": {},
   "source": [
    "In the **code cell** below, we set the main constants.\n",
    "\n",
    "* **main**:\n",
    "    * `PACKMOL_BIN`: is the path to packmol binary.\n",
    "    * `SEED`: seed number for the random generator.\n",
    "    * `SEI_COMPONENTS_INFO`: is a `pandas.DataFrame` containing all the information about the inorganic salts \n",
    "        forming the SEI.\n",
    "    * `ANODE_SHAPE`: the numerical box for the anode.\n",
    "* *optional*:\n",
    "    * `NATOMS`: number of atoms inside the numerical box.\n",
    "    * `NATOMS_GRAIN`: max number of atoms per grain.\n",
    "    * `CUTTING_PLANES`: list of Miller indices with all possible cutting surfaces to use to get an SEI grain.\n",
    "    * `PLOT_DIR`: plots storing directory.\n",
    "    * `CONF_DIR`: configurations storing directory, i.e., where all the \".XYZ\" files are stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a64d24d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "PACKMOL_BIN = \"./packmol/packmol\"  # Optional, if packmol is bin folder is in the PATH you can set it as `None`\n",
    "SEED = 42\n",
    "ANODE_SHAPE = [35.0, 35.0, 20.0]  # Anode shape to use\n",
    "# (the number are indicative since the shape will adjust to the closest supercell)\n",
    "\n",
    "SEI_COMPONENTS_INFO = pd.DataFrame(\n",
    "    index=[\"LiF\", \"Li2CO3\", \"Li2O\"],  # Empircal chemical formulat of the components (rows)\n",
    "    columns=[  # name of each data (columns)\n",
    "        \"molar fraction\",  # The cristals molar fraction\n",
    "        \"size distribution (mean)\",  # The mean grain diameters\n",
    "        \"size distribution (std)\",  # The std grain diameters\n",
    "    ],\n",
    "    data=[\n",
    "        [24 / 35, 15.0, 3.0],  # LiF\n",
    "        [7 / 35, 13.0, 3.0],  # Li2CO3\n",
    "        [4 / 35, 10.0, 2.0],  # Li2O\n",
    "    ],\n",
    ")\n",
    "\n",
    "NATOMS = 5000  # Number of atoms inside the numerical box\n",
    "NATOMS_GRAIN = 700  # Max number of atoms per grain\n",
    "CUTTING_PLANES = [  # List of Miller indices with all possible cutting surfaces to use to get an SEI grain\n",
    "    (1, 0, 0),\n",
    "    (1, 1, 0),\n",
    "    (1, 1, 1),\n",
    "    (1, 2, 0),\n",
    "    (1, 3, 0),\n",
    "    (1, 2, 1),\n",
    "    (1, 3, 1),\n",
    "]\n",
    "# Outputs\n",
    "PLOT_DIR = \"out/plots\"\n",
    "makedir(PLOT_DIR)\n",
    "CONF_DIR = \"out/sei_conf\"\n",
    "makedir(CONF_DIR)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1c857b3",
   "metadata": {},
   "source": [
    "showing the data for the SEI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4106cc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "SEI_COMPONENTS_INFO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abstract-thumbnail",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Get Salts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "255bab67",
   "metadata": {},
   "source": [
    "First thing first, we are querying the [Materials Project](https://materialsproject.org/) database to get a more stable(i.e., the one with minimum formation energy) crystal unit for each salt in the DataFrame `SEI_COMPONENTS_INFO`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d4378d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "unitcell_mp = {}\n",
    "unitcell_ase = {}\n",
    "\n",
    "for formula in SEI_COMPONENTS_INFO.index:\n",
    "    uc_mp, uc_ase = get_stable_crystal(formula)\n",
    "    unitcell_mp[formula] = uc_mp\n",
    "    unitcell_ase[formula] = uc_ase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54adaaab",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nuc = len(unitcell_ase)\n",
    "fig = plt.figure(figsize=[4 * min(Nuc, 3), 4 * np.floor(Nuc / 3)], dpi=100, facecolor=\"white\")\n",
    "\n",
    "with plt.style.context(\"seaborn-whitegrid\"):\n",
    "    spec = gridspec.GridSpec(ncols=min(Nuc, 3), nrows=int(np.floor(Nuc / 3)), figure=fig, wspace=0.3)\n",
    "    ax = []\n",
    "    for i, (formula, atoms) in enumerate(unitcell_ase.items()):\n",
    "        ax.append(fig.add_subplot(spec[i]))\n",
    "        plot_atoms(atoms, ax[-1], radii=0.4, rotation=\"0x,0y,0z\")\n",
    "        ax[-1].set_title(f\"{formula}\")\n",
    "        ax[-1].set_xlabel(r\"x [$\\mathrm{\\AA}$]\")\n",
    "        ax[-1].set_ylabel(r\"y [$\\mathrm{\\AA}$]\")\n",
    "\n",
    "fig.savefig(os.path.join(PLOT_DIR, \"1-crystals-unitcell.png\"), dpi=300, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "marine-yukon",
   "metadata": {},
   "source": [
    "## Random sempling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abandoned-china",
   "metadata": {},
   "source": [
    "### Set up sampling space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11c84363",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "pcg64 = PCG64(seed=SEED)  # Set random number generator (PCG 64) and random state seed\n",
    "rg = Generator(pcg64)\n",
    "rg.standard_normal()  # Since we do not know the experimental size distribution, we use a gues normal PDF\n",
    "\n",
    "\n",
    "def grain_d(d_mean: float, d_std: float) -> float:\n",
    "    \"\"\"Get an SEI's grain size randomly from a gaussian distribution.\n",
    "\n",
    "    This is an auxiliary function and can be replaced with a function that sample\n",
    "    from an experimental data or different PDF.\n",
    "    Args:\n",
    "        d_mean (float): mean size grain size.\n",
    "        d_std (float): grain size standard deviation.\n",
    "\n",
    "    Returns:\n",
    "        float: random sampled size.\n",
    "    \"\"\"\n",
    "    d = rg.normal(loc=d_mean, scale=d_std, size=1)[0]\n",
    "    return d\n",
    "\n",
    "\n",
    "semplings_functions = []\n",
    "d_mean = SEI_COMPONENTS_INFO[\"size distribution (mean)\"]\n",
    "d_std = SEI_COMPONENTS_INFO[\"size distribution (std)\"]\n",
    "for d_mean_, d_std_ in zip(d_mean, d_std):\n",
    "    semplings_functions.append(partial(grain_d, d_mean_, d_std_))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "variable-spray",
   "metadata": {},
   "source": [
    "### SEI's grains sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84b25848",
   "metadata": {},
   "source": [
    "With the function `random_sei_grains`, we build all the grains (nanoparticles) forming the SEI using the Wulff construction. The sampling is constrained, forcing the final molar fraction and the total number of atoms.\n",
    "\n",
    "Below are all the function inputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db6b7f93",
   "metadata": {},
   "outputs": [],
   "source": [
    "help(random_sei_grains)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adapted-posting",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "species_unitcell = list(unitcell_mp.values())\n",
    "species_fractions = SEI_COMPONENTS_INFO[\"molar fraction\"]\n",
    "\n",
    "out_species, out_d, out_vol, out_species_fraction, out_grains = random_sei_grains(\n",
    "    NATOMS,\n",
    "    species_unitcell,\n",
    "    species_fractions,\n",
    "    semplings_functions,\n",
    "    species_fraction_tol=0.05,\n",
    "    Ngrains_max=NATOMS_GRAIN,\n",
    "    cutting_planes=CUTTING_PLANES,\n",
    "    n_planes=2,\n",
    "    seed=SEED,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bigger-thought",
   "metadata": {},
   "source": [
    "Store the results (as pickle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "compound-joint",
   "metadata": {
    "vscode": {
     "languageId": "markdown"
    }
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(CONF_DIR, \"0-grain_list.pkl\"), \"wb\") as file_:\n",
    "    pickle.dump([out_species, out_d, out_vol, out_species_fraction, out_grains], file_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "existing-radius",
   "metadata": {},
   "source": [
    "### Final results analisy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "reserved-municipality",
   "metadata": {},
   "source": [
    "Load result from previus step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "victorian-judge",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(CONF_DIR, \"0-grain_list.pkl\"), \"rb\") as file_:\n",
    "    out_species, out_d, out_vol, out_species_fraction, out_grains = pickle.load(file_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18489e79",
   "metadata": {},
   "source": [
    "Plotting particle size theoretical and computed distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "composite-appraisal",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[7, 4], dpi=100)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "labels = [r\"$\\mathrm{LiF}$\", r\"$\\mathrm{Li_2CO_3}$\", r\"$\\mathrm{Li_2O}$\"]\n",
    "bins = np.linspace(out_d.min(), out_d.max(), 8)\n",
    "d_x = np.linspace(out_d.min(), out_d.max(), 200)\n",
    "\n",
    "with plt.style.context(\"seaborn-whitegrid\"):\n",
    "    hist_d, _ = np.histogram(out_d, bins=bins, density=True)\n",
    "    ax.bar(\n",
    "        bins[:-1], np.round(hist_d * len(out_species)), width=bins[1] - bins[0], align=\"edge\", label=\"total\", alpha=0.6\n",
    "    )\n",
    "    tot_dist = np.zeros(len(d_x))\n",
    "    for i in np.unique(out_species):\n",
    "        tot_dist += stats.norm.pdf(d_x, d_mean[i], d_std[i]) * np.sum(out_species == i)\n",
    "    ax.plot(d_x, tot_dist)\n",
    "\n",
    "    for i in np.unique(out_species):\n",
    "        hist_d, _ = np.histogram(out_d[out_species == i], bins=bins, density=True)\n",
    "        ax.bar(\n",
    "            bins[:-1],\n",
    "            np.round(hist_d * np.sum(out_species == i)),\n",
    "            width=bins[1] - bins[0],\n",
    "            align=\"edge\",\n",
    "            label=labels[i],\n",
    "            alpha=0.6,\n",
    "        )\n",
    "        ax.plot(d_x, stats.norm.pdf(d_x, d_mean[i], d_std[i]) * np.sum(out_species == i))\n",
    "\n",
    "    ax.legend()\n",
    "    ax.grid()\n",
    "    ax.set_xlabel(r\"Characteristic length $d$ [$\\mathrm{\\AA}$]\")\n",
    "    ax.set_ylabel(\"Number of grains [$-$]\")\n",
    "plt.show()\n",
    "fig.savefig(os.path.join(PLOT_DIR, \"2-grain_distributions.png\"), bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38ebfdef",
   "metadata": {},
   "source": [
    "Plotting particle molar fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "international-column",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[4, 7], dpi=100)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "labels = [r\"$\\mathrm{LiF}$\", r\"$\\mathrm{Li_2CO_3}$\", r\"$\\mathrm{Li_2O}$\"]\n",
    "x = np.arange(len(out_species_fraction))\n",
    "\n",
    "with plt.style.context(\"seaborn-whitegrid\"):\n",
    "    ax.bar(x + 0.15, species_fractions / species_fractions.sum() * 100, width=0.3, label=\"target\", alpha=0.6)\n",
    "    ax.bar(x - 0.15, out_species_fraction * 100, width=0.3, label=\"final\", alpha=0.6)\n",
    "\n",
    "    ax.legend()\n",
    "    ax.grid()\n",
    "    ax.set_xticks(x)\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.set_xlabel(\"species\")\n",
    "    ax.set_ylabel(\"mole fraction [%]\")\n",
    "plt.show()\n",
    "fig.savefig(os.path.join(PLOT_DIR, \"3-grain_mole_fraction.png\"), bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c91d98f",
   "metadata": {},
   "source": [
    "Plotting particle mass fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "protective-infrastructure",
   "metadata": {},
   "outputs": [],
   "source": [
    "mass = np.array([np.sum(atoms.get_masses()) for atoms in unitcell_ase.values()])\n",
    "molecules_in_unit_cell = np.array([get_gcd_pedices(system.formula) for system in species_unitcell])\n",
    "mass = mass / molecules_in_unit_cell\n",
    "atoms_in_unit_cell = np.array([len(system.sites) for system in species_unitcell])\n",
    "out_atoms = np.zeros(len(np.unique(out_species)))\n",
    "\n",
    "for i in np.unique(out_species):\n",
    "    for particle in out_grains[out_species == i]:\n",
    "        out_atoms[i] += len(particle)\n",
    "        out_mass_fraction = out_atoms / atoms_in_unit_cell * molecules_in_unit_cell\n",
    "out_mass_fraction = out_mass_fraction / out_mass_fraction.sum() * mass\n",
    "\n",
    "fig = plt.figure(figsize=[4, 7], dpi=100)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "labels = [r\"$\\mathrm{LiF}$\", r\"$\\mathrm{Li_2CO_3}$\", r\"$\\mathrm{Li_2O}$\"]\n",
    "x = np.arange(len(out_species_fraction))\n",
    "\n",
    "with plt.style.context(\"seaborn-whitegrid\"):\n",
    "    ax.bar(x + 0.15, species_fractions / species_fractions.sum() * mass, width=0.3, label=\"target\", alpha=0.6)\n",
    "    ax.bar(x - 0.15, out_mass_fraction, width=0.3, label=\"final\", alpha=0.6)\n",
    "\n",
    "    ax.legend()\n",
    "    ax.grid()\n",
    "    ax.set_xticks(x)\n",
    "    ax.set_xticklabels(labels)\n",
    "    ax.set_xlabel(\"species\")\n",
    "    ax.set_ylabel(\"weight fraction [g/mol]\")\n",
    "plt.show()\n",
    "fig.savefig(os.path.join(PLOT_DIR, \"4-grain_weight_fraction.png\"), bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78b167d0-14e4-4b6e-bea3-38dc08f7764f",
   "metadata": {},
   "source": [
    "## Get additional characteristics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8711a06b-1ea7-42b8-94eb-26618c4db726",
   "metadata": {},
   "source": [
    "### Surface Atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14ec6490",
   "metadata": {},
   "source": [
    "Computing the coordination number, we identify the atoms in bulk or on the surface, which will be helpful for studying the *grain boundaries*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d834298f-c326-40ba-b24c-4efd661c4a78",
   "metadata": {},
   "outputs": [],
   "source": [
    "for grain_ in out_grains:\n",
    "    bulk_atoms, score = get_bulk_atoms(grain_, threshold=0.5, cutoff=5.0)\n",
    "    grain_.info[\"bulk_atoms\"] = bulk_atoms\n",
    "    grain_.info[\"bulk_atoms_score\"] = score"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dae0a5aa-b793-4133-b007-ab3a69cca80d",
   "metadata": {},
   "source": [
    "To visualize the result, we change the symbols for the atoms in bulk as \"C\" (Carbon) and on the surface as \"O\" (Oxgine). In this way, the first will show as gray and the second as red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aaa52e5-c1aa-4cc8-8d7c-c33d90653f64",
   "metadata": {},
   "outputs": [],
   "source": [
    "grain_id = 11\n",
    "clust_ = out_grains[grain_id]\n",
    "positions = clust_.get_positions()\n",
    "symbols = [\"C\" if x else \"O\" for x in clust_.info[\"bulk_atoms\"]]\n",
    "col_clust = Atoms(symbols=symbols, positions=positions)\n",
    "\n",
    "v = view(col_clust, viewer=\"ngl\")\n",
    "v.view.clear()\n",
    "v.view.add_ball_and_stick(col_clust)\n",
    "v.view.add_unitcell(col_clust)\n",
    "rot = Rotation.from_euler(\"xyz\", [-20, 10, -10], degrees=True)\n",
    "v.view.control.rotate(rot.as_quat().tolist())\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f2c6afb",
   "metadata": {},
   "source": [
    "Rendering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66993b22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# render\n",
    "img = v.view.render_image(factor=7, trim=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a2fbd73",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save image\n",
    "img_path = os.path.join(PLOT_DIR, f\"render1-grain{grain_id}_bulk_surf.png\")\n",
    "with open(img_path, \"wb\") as fh:\n",
    "    fh.write(img.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee286812",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display\n",
    "display(Image(img_path, width=400))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13f780bf-3d0f-4842-9d98-3cb0b75fb95e",
   "metadata": {},
   "source": [
    "Seve the results (as pickle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "509df46f-e9a1-4d77-ae52-089eb887c4ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.path.join(CONF_DIR, \"1-grain_list_with_info.pkl\"), \"wb\") as file_:\n",
    "    pickle.dump([out_species, out_d, out_vol, out_species_fraction, out_grains], file_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "liberal-baptist",
   "metadata": {},
   "source": [
    "## Pack all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "grateful-carry",
   "metadata": {},
   "source": [
    "### Build Anode"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d131b112",
   "metadata": {},
   "source": [
    "Since we are going to add the anode to the simulation box, the length and width of the simulation box are constrained by this latter.\n",
    "Thus first, we build two possible anodes, the graphite one and the silicon one.\n",
    "The user can later choose which to use."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b566cc0-9465-4b3f-9509-90aff680d5a0",
   "metadata": {},
   "source": [
    "#### $C_6$ Anode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "japanese-breast",
   "metadata": {},
   "outputs": [],
   "source": [
    "class GrFactory(HexagonalFactory):  # make a new crsital object from the HexagonalFactory\n",
    "    \"A factory for creating Graphite lattices.\"\n",
    "    xtal_name = \"graphite\"\n",
    "    bravais_basis = [[0, 0, 0], [1.0 / 3.0, 2.0 / 3.0, 0], [1.0 / 3.0, 2.0 / 3.0, 0.5], [2.0 / 3.0, 1.0 / 3.0, 0.5]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "better-assist",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gr = GrFactory()\n",
    "a = 2.461  # [Å]\n",
    "c = 6.708  # [Å]\n",
    "Lx, Ly, Lz = ANODE_SHAPE  # [Å]\n",
    "Nx, Ny, Nz = [int(np.round(Lx / c)), int(np.round(Ly / a)), int(np.round(Lz / a / 3**0.5))]\n",
    "GrAnode = Gr(\n",
    "    latticeconstant={\"a\": 2.461, \"c\": 6.708},\n",
    "    symbol=\"C\",\n",
    "    size=(Nx, Ny, Nz),\n",
    "    directions=[[0, 0, 0, 1], [1, -2, 1, 0], [1, 0, -1, 0]],\n",
    "    pbc=True,\n",
    ")\n",
    "GrAnode.info[\"name\"] = \"Anode\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e85c34a1-5679-4326-8f54-e2ac19db53f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = view(GrAnode, viewer=\"ngl\")\n",
    "v.view.clear()\n",
    "v.view.add_ball_and_stick(GrAnode)\n",
    "v.view.add_unitcell(GrAnode)\n",
    "# rot = Rotation.from_euler('xyz', [90-35.264, 45, -90], degrees=True)\n",
    "rot = Rotation.from_euler(\"xyz\", [90, 0, -90], degrees=True)\n",
    "v.view.control.rotate(rot.as_quat().tolist())\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ea48a99",
   "metadata": {},
   "source": [
    "Rendering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48dc24f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# render\n",
    "img = v.view.render_image(factor=7, trim=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c9dfe2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save image\n",
    "img_path = os.path.join(PLOT_DIR, \"render2-C6_anode.png\")\n",
    "with open(img_path, \"wb\") as fh:\n",
    "    fh.write(img.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f5a6576",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display\n",
    "display(Image(img_path, width=200))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80be41e5-7445-4915-9458-e38d3c50d364",
   "metadata": {},
   "source": [
    "#### $SI$ Anode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef565ec9-45db-48ab-99bd-c1ccda37c45d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Diamond = DiamondFactory()\n",
    "a = 5.43  # [Å]\n",
    "Nslab = 18\n",
    "da_z = a / 4\n",
    "Lx, Ly, Lz = ANODE_SHAPE  # [Å]\n",
    "Nx, Ny, Nz = [int(np.round(Lx / a * 2**0.5)), int(np.round(Ly / a * 2**0.5)), int(np.round(Lz / a))]\n",
    "# Nx, Ny, Nz = [9, 9, 5]\n",
    "SiAnode = Diamond(\n",
    "    latticeconstant=a, symbol=\"Si\", size=(Nx, Ny, Nz), directions=[[1, 1, 0], [-1, 1, 0], [0, 0, 1]], pbc=True\n",
    ")\n",
    "SiAnode.info[\"name\"] = \"Anode\"\n",
    "SiAnode.cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "431a642d",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = view(SiAnode, viewer=\"ngl\")\n",
    "v.view.clear()\n",
    "v.view.add_ball_and_stick(SiAnode)\n",
    "v.view.add_unitcell(SiAnode)\n",
    "# rot = Rotation.from_euler('xyz', [90-35.264, 45, -90], degrees=True)\n",
    "rot = Rotation.from_euler(\"xyz\", [90, 0, -90], degrees=True)\n",
    "v.view.control.rotate(rot.as_quat().tolist())\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ded76048",
   "metadata": {},
   "source": [
    "Rendering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0929aa17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# render\n",
    "img = v.view.render_image(factor=7, trim=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9f2f2f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save image\n",
    "img_path = os.path.join(PLOT_DIR, \"render3-Si_anode.png\")\n",
    "with open(img_path, \"wb\") as fh:\n",
    "    fh.write(img.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2c7c51d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display\n",
    "display(Image(img_path, width=200))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "atlantic-missile",
   "metadata": {},
   "source": [
    "### Pack SEI grains"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3542acb1",
   "metadata": {},
   "source": [
    "We are placing using the code [`packmol`](http://leandro.iqm.unicamp.br/m3g/packmol/home.shtml) stochastically and \n",
    "rotating and translating each grain to find the optima packing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1020a9d-39d1-4e72-ad45-15e256b3068f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a Anode\n",
    "ANODE = SiAnode  # or GrAnode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "convenient-nigeria",
   "metadata": {},
   "outputs": [],
   "source": [
    "porosity = 2.0\n",
    "tot_vol = out_vol.sum() * (1 + porosity)\n",
    "lc, wc, hc = ANODE.get_cell().diagonal()  # [Å]\n",
    "clearance = 3.0  # [Å]\n",
    "box = [lc, wc, tot_vol / (lc - clearance) / (wc - clearance) + clearance]\n",
    "print(f\"Intial SEI's box size {box}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "extended-lending",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 2.5\n",
    "SEI_box = packmol(\n",
    "    [\n",
    "        PackmolStructure(\n",
    "            grain_,\n",
    "            number=1,\n",
    "            instructions=[\n",
    "                \"inside box {:.16f} {:.16f} {:.16f} {:.16f} {:.16f} {:.16f}\".format(\n",
    "                    clearance / 2,\n",
    "                    clearance / 2,\n",
    "                    clearance / 2,\n",
    "                    box[0] - clearance / 2,\n",
    "                    box[1] - clearance / 2,\n",
    "                    box[2] - clearance / 2,\n",
    "                )\n",
    "            ],\n",
    "        )\n",
    "        for grain_ in out_grains\n",
    "    ],\n",
    "    packmol_bin=PACKMOL_BIN,\n",
    "    verbose=2,\n",
    "    n_lines=20,\n",
    "    cell=box,\n",
    "    tolerance=tolerance,\n",
    "    nloop=200,\n",
    "    seed=SEED,  # maxit=50,#, maxit=100,  cell=box_size,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nuclear-indie",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = view(SEI_box, viewer=\"ngl\")\n",
    "v.view.clear()\n",
    "v.view.add_ball_and_stick(SEI_box)\n",
    "v.view.add_unitcell(SEI_box)\n",
    "rot = Rotation.from_euler(\"xyz\", [90, 0, -90], degrees=True)\n",
    "v.view.control.rotate(rot.as_quat().tolist())\n",
    "v.view.display(gui=False)\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21fd820c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# render\n",
    "img = v.view.render_image(factor=7, trim=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac9dcd4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save image\n",
    "img_path = os.path.join(PLOT_DIR, \"render4-SEI_step0_packmol.png\")\n",
    "with open(img_path, \"wb\") as fh:\n",
    "    fh.write(img.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abfba333",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display\n",
    "display(Image(img_path, height=400))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9194e9be",
   "metadata": {},
   "source": [
    "Storing only the SEI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "040119f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "write(os.path.join(CONF_DIR, \"sei_box.xyz\"), SEI_box)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "749b2498",
   "metadata": {},
   "source": [
    "### add Anode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4daf353",
   "metadata": {},
   "outputs": [],
   "source": [
    "h_z = ANODE.get_cell_lengths_and_angles()[2]\n",
    "\n",
    "sei_tmp = SEI_box.copy()\n",
    "sei_tmp.translate([0, 0, h_z + 1.0])\n",
    "SEI_Anode = sei_tmp + ANODE\n",
    "new_box = SEI_box.get_cell_lengths_and_angles()[:3]\n",
    "new_box[2] += h_z + 1.0\n",
    "SEI_Anode.set_cell(new_box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6125fc0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = view(SEI_Anode, viewer=\"ngl\")\n",
    "v.view.clear()\n",
    "v.view.add_ball_and_stick(SEI_Anode)\n",
    "v.view.add_unitcell(SEI_Anode)\n",
    "rot = Rotation.from_euler(\"xyz\", [90, 0, -90], degrees=True)\n",
    "v.view.control.rotate(rot.as_quat().tolist())\n",
    "v.view.display(gui=False)\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71fdaa5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# render\n",
    "img = v.view.render_image(factor=7, trim=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85bdbcfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save image\n",
    "img_path = os.path.join(PLOT_DIR, \"render5-SEI_step1_addanode.png\")\n",
    "with open(img_path, \"wb\") as fh:\n",
    "    fh.write(img.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cd1dd22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display\n",
    "display(Image(img_path, height=400))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11282319",
   "metadata": {},
   "source": [
    "Storing only the SEI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "500b9048",
   "metadata": {},
   "outputs": [],
   "source": [
    "write(os.path.join(CONF_DIR, \"sei_anode.xyz\"), SEI_Anode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c73c596-bc49-4b77-9eaa-51048c71be6b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}