sample.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Welcome to the OpenKIM tutorial at MACH 2023\n",
    "\n",
    "Let's take a look at what's inside the LAMMPS file. Run the cell below to see:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cat in.lammps"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!lmp_mpi < in.lammps"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's do a calculation in ASE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.calculators.kim import KIM\n",
    "from ase.lattice.cubic import FaceCenteredCubic\n",
    "from kim_query import get_lattice_constant_cubic\n",
    "from ase.optimize import LBFGS\n",
    "# Define KIM model and get Al fcc lattice parameter for this potential\n",
    "model = \"EAM_Dynamo_ErcolessiAdams_1994_Al__MO_123629422045_005\"\n",
    "calc = KIM(model)\n",
    "a0 = get_lattice_constant_cubic([model], [\"fcc\"], [\"Al\"], [\"angstrom\"])[0]\n",
    "# Set up fcc crystal\n",
    "atoms = FaceCenteredCubic(\"Al\", latticeconstant=a0,size=(9,9,9),pbc=(1,1,0))\n",
    "atoms.calc = calc\n",
    "# Compute energy\n",
    "dyn = LBFGS(atoms)\n",
    "dyn.run(fmax=1e-8)\n",
    "print(atoms.get_potential_energy())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convex hull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.spatial import ConvexHull\n",
    "from scipy.spatial import ConvexHull, convex_hull_plot_2d\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "from curses.ascii import isalpha, isdigit\n",
    "import requests\n",
    "from typing import List\n",
    "\n",
    "#\"\"\"\n",
    "\n",
    "def get_stoich_reduced_list_from_prototype(prototype_label: str) -> List[int]:\n",
    "    \"\"\"\n",
    "    Get numerical list of stoichiometry from prototype label, i.e. \"AB3\\_....\" -> [1,3]\n",
    "\n",
    "    Args:\n",
    "        prototype_label:\n",
    "            AFLOW prototype label\n",
    "\n",
    "    Returns:\n",
    "        List of reduced stoichiometric numbers\n",
    "    \"\"\"                        \n",
    "    stoich_reduced_formula = prototype_label.split(\"_\")[0]\n",
    "    stoich_reduced_list=[]\n",
    "    stoich_reduced_curr = None\n",
    "    for char in stoich_reduced_formula:\n",
    "        if isalpha(char):\n",
    "            if stoich_reduced_curr is not None:\n",
    "                if stoich_reduced_curr == 0:\n",
    "                    stoich_reduced_curr = 1\n",
    "                stoich_reduced_list.append(stoich_reduced_curr)\n",
    "            stoich_reduced_curr = 0\n",
    "        else:\n",
    "            assert isdigit(char)                            \n",
    "            stoich_reduced_curr*=10 # will throw an error if we haven't encountered an alphabetical letter, good\n",
    "            stoich_reduced_curr+=int(char)\n",
    "    # write final number                    \n",
    "    if stoich_reduced_curr == 0:\n",
    "        stoich_reduced_curr = 1\n",
    "    stoich_reduced_list.append(stoich_reduced_curr)    \n",
    "    return stoich_reduced_list\n",
    "\n",
    "response = requests.post(\"https://query.openkim.org/api\",data={'query': '{\"property-id\":\"tag:staff@noreply.openkim.org,2023-02-21:property/cohesive-potential-energy-crystal\"}', 'fields': '{\"stoichiometric-species.source-value\":1,\"prototype-label.source-value\":1,\"formation-potential-energy-per-atom.source-value\":1}', 'database': 'data'}).json()\n",
    "unary_results = {}\n",
    "compound_results = {}\n",
    "\n",
    "for result in response:\n",
    "    species_combo = tuple(result[\"stoichiometric-species\"][\"source-value\"])\n",
    "    result_dict = {\n",
    "                \"prototype-label\": result[\"prototype-label\"][\"source-value\"],\n",
    "                \"formation-potential-energy-per-atom\": result[\"formation-potential-energy-per-atom\"][\"source-value\"],\n",
    "                \"stoichiometric-species\": result[\"stoichiometric-species\"][\"source-value\"]\n",
    "            }\n",
    "    if len(species_combo)==1:\n",
    "        if species_combo[0] in unary_results:\n",
    "            unary_results[species_combo[0]].append(result_dict)\n",
    "        else:\n",
    "            unary_results[species_combo[0]] = [result_dict]\n",
    "    else:\n",
    "        if species_combo in compound_results:\n",
    "            compound_results[species_combo].append(result_dict)\n",
    "        else:\n",
    "            compound_results[species_combo] = [result_dict]\n",
    "\n",
    "#loop over all compound species_combos\n",
    "vertices_dict = {}\n",
    "for species_combo in compound_results:\n",
    "    arity = len(species_combo)\n",
    "    compounds = compound_results[species_combo]\n",
    "    num_compounds = len(compounds)\n",
    "    # get reference energies\n",
    "    elemental_reference_list = []    \n",
    "    elemental_energies = []\n",
    "    for element in species_combo:\n",
    "        elemental_reference = min(unary_results[element],key=lambda tr:tr[\"formation-potential-energy-per-atom\"])\n",
    "        elemental_reference_list.append(elemental_reference)\n",
    "        elemental_energies.append(elemental_reference[\"formation-potential-energy-per-atom\"])\n",
    "\n",
    "    \n",
    "\n",
    "    stoich_energy_array = np.zeros((num_compounds+2,arity))\n",
    "    stoich_energy_array[-1,-2]=1\n",
    "\n",
    "    for (i,compound) in enumerate(compounds):\n",
    "        stoich_reduced_list = get_stoich_reduced_list_from_prototype(compound[\"prototype-label\"])\n",
    "        stoich_fractional_list = [elem_num/sum(stoich_reduced_list) for elem_num in stoich_reduced_list]\n",
    "        reference_energy_per_atom = np.dot(stoich_fractional_list,elemental_energies)\n",
    "        #set the coordinates to the stoichiometry, skipping first element\n",
    "        stoich_energy_array[i,:-1]=stoich_fractional_list[1:]\n",
    "        stoich_energy_array[i,-1]=compound[\"formation-potential-energy-per-atom\"]-reference_energy_per_atom        \n",
    "\n",
    "    \n",
    "        \n",
    "    hull = ConvexHull(stoich_energy_array)\n",
    "    vertices_dict[species_combo] = hull.vertices.shape[0]\n",
    "    _ = convex_hull_plot_2d(hull)\n",
    "    plt.show()\n",
    "    plt.close(\"all\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.10.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}