{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU",
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "This notebook can be run on Google Colab.\n",
        "\n",
        "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/nanotheorygroup/kaldo/blob/master/docs/docsource/crystal_presentation.ipynb)\n",
        "\n",
        "\n",
        "\n",
        "In Colab, you can enable the GPU acceleration from `Edit` > `Notebook Settings` > `Accelerator` > `GPU`."
      ],
      "metadata": {
        "id": "XBeaBnL82NOM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Silicon diamond tutorial\n",
        "\n",
        "## Complie [LAMMPS](https://github.com/lammps/lammps) as shared-library of python ($\\sim$ 8 min)"
      ],
      "metadata": {
        "id": "wW4-RGLk2QbJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Sudo apt-get softwares\n",
        "! apt-get update\n",
        "! apt install -y cmake build-essential git ccache openmpi-bin libopenmpi-dev python3.10-venv\n",
        "! echo \"Sudo apt-get finishes!\"\n",
        "\n",
        "# Upgrade python packages\n",
        "! pip install --upgrade pip\n",
        "! pip install numpy torch scipy virtualenv psutil pandas tabulate mpi4py Cython\n",
        "! echo \"Python packages installation finishes!\"\n",
        "\n",
        "# Build lammps with cmake\n",
        "%cd /content\n",
        "!rm -rf lammps\n",
        "! git clone https://github.com/lammps/lammps.git lammps\n",
        "%cd /content/lammps\n",
        "! rm -rf build\n",
        "! mkdir build\n",
        "%cd build\n",
        "! cmake ../cmake -DLAMMPS_EXCEPTIONS=yes \\\n",
        "               -DBUILD_SHARED_LIBS=yes \\\n",
        "               -DMLIAP_ENABLE_PYTHON=yes \\\n",
        "               -DPKG_PYTHON=yes \\\n",
        "               -DPKG_MANYBODY=yes \\\n",
        "               -DPKG_KSPACE=yes \\\n",
        "               -DPKG_PHONON=yes \\\n",
        "               -DPYTHON_EXECUTABLE:FILEPATH=`which python`\n",
        "\n",
        "# Complie lammps as share-libary of python\n",
        "! make -j 2\n",
        "! make install-python\n",
        "! echo \"LAMMPS compilation done!\"\n",
        "\n",
        "# Redirect back to main folder\n",
        "%cd /content/"
      ],
      "metadata": {
        "id": "OIAXE-KV6BF8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Remote fetch  and install source code from Github"
      ],
      "metadata": {
        "id": "xDuRPLMA5msm"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6G-DsQwX1nl4"
      },
      "outputs": [],
      "source": [
	"! git lfs install --skip-repo --skip-smudge\n",
        "! pip install git+https://github.com/nanotheorygroup/kaldo"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Remote fetch supplyment files\n"
      ],
      "metadata": {
        "id": "T2lEjnRw2dcB"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Remote fetch kaldo resources from drop box\n",
        "! wget https://www.dropbox.com/s/bvw0qcxy397g25q/kaldo_resources.zip?dl=0\n",
        "! mv kaldo_resources.zip?dl=0 kaldo_resources.zip\n",
        "! unzip kaldo_resources.zip\n",
        "\n",
        "# Unzip files\n",
        "!unzip forcefields.zip\n",
        "\n",
        "# Clean workspace\n",
        "! rm -r forcefields.zip\n",
        "! rm -r structure_a_si_512.zip\n",
        "! rm -r kaldo_resources.zip\n",
        "! rm -r sample_data"
      ],
      "metadata": {
        "id": "yddHjxOy2dv5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Thermal transport simulation for silicon-bulk"
      ],
      "metadata": {
        "id": "nvL1t2Wt5e4X"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from ase.build import bulk\n",
        "from ase.calculators.lammpslib import LAMMPSlib\n",
        "from kaldo.forceconstants import ForceConstants\n",
        "import numpy as np\n",
        "\n",
        "# We start from the atoms object\n",
        "atoms = bulk('Si', 'diamond', a=5.432)\n",
        "\n",
        "# Config super cell and calculator input\n",
        "supercell = np.array([3, 3, 3])\n",
        "lammps_inputs = {\n",
        "      'lmpcmds': [\n",
        "          'pair_style tersoff',\n",
        "          'pair_coeff * * forcefields/Si.tersoff Si'],\n",
        "\n",
        "      'log_file': 'lammps-si-bulk.log',\n",
        "      'keep_alive':True}\n",
        "\n",
        "# Create a finite difference object\n",
        "forceconstants_config  = {'atoms':atoms,'supercell': supercell,'folder':'fd'}\n",
        "forceconstants = ForceConstants(**forceconstants_config)\n",
        "\n",
        "# Compute 2nd and 3rd IFCs with the defined calculators\n",
        "forceconstants.second.calculate(LAMMPSlib(**lammps_inputs), delta_shift=1e-3)\n",
        "forceconstants.third.calculate(LAMMPSlib(**lammps_inputs), delta_shift=1e-3)"
      ],
      "metadata": {
        "id": "hXkcBeU5_QwG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Create phonons object\n"
      ],
      "metadata": {
        "id": "eqTobAF_CB-m"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from kaldo.phonons import Phonons\n",
        "\n",
        "# Define k-point grids, temperature\n",
        "# and the assumption for the \n",
        "# phonon poluation (i.e classical vs. quantum)\n",
        "k = 7\n",
        "kpts = [k, k, k]\n",
        "temperature = 300\n",
        "is_classic = False\n",
        "k_label = str(k) + '_' + str(k) + '_' + str(k)\n",
        "\n",
        "# Create a phonon object\n",
        "phonons = Phonons(forceconstants=forceconstants,\n",
        "                kpts=kpts,\n",
        "                is_classic=is_classic,\n",
        "                temperature=temperature,\n",
        "                folder='si-bulk-ald-' + k_label,\n",
        "                storage='numpy')"
      ],
      "metadata": {
        "id": "UYwM7PVXCEa0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Calculate conductivities for infinite-size sample\n"
      ],
      "metadata": {
        "id": "MLZRdZRwCIG0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from kaldo.conductivity import Conductivity\n",
        "\n",
        "# Calculate conductivity  with direct inversion approach (inverse)\n",
        "print('\\n')\n",
        "inv_cond_matrix = (Conductivity(phonons=phonons, method='inverse').conductivity.sum(axis=0))\n",
        "print('Inverse conductivity (W/mK): %.3f'%(np.mean(np.diag(inv_cond_matrix))))\n",
        "print(inv_cond_matrix)\n",
        "\n",
        "# Calculate conductivity  with  relaxation time approximation (rta)\n",
        "print('\\n')\n",
        "rta_cond_matrix = Conductivity(phonons=phonons, method='rta').conductivity.sum(axis=0)\n",
        "print('Rta conductivity (W/mK): %.3f'%(np.mean(np.diag(rta_cond_matrix))))\n",
        "print(rta_cond_matrix)\n",
        "# Calculate conductivity  with  self-consistent approach (sc)\n",
        "\n",
        "print('\\n')\n",
        "sc_cond_matrix = Conductivity(phonons=phonons, method='sc',n_iterations=20).conductivity.sum(axis=0)\n",
        "print('Self-consistent conductivity (W/mK): %.3f'%(np.mean(np.diag(sc_cond_matrix))))\n",
        "print(sc_cond_matrix)"
      ],
      "metadata": {
        "id": "QWs_1LgICIUd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Visualize harmonic properties using built-in plotter\n"
      ],
      "metadata": {
        "id": "FNMtnrSqCMel"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import kaldo.controllers.plotter as plotter\n",
        "import matplotlib.pyplot as plt\n",
        "plt.style.use('seaborn-poster')\n",
        "\n",
        "# Plot dispersion relation and group velocity in each direction\n",
        "plotter.plot_dispersion(phonons,n_k_points=int(k_label))\n",
        "print('\\n')"
      ],
      "metadata": {
        "id": "O-TT1ZgXCMvb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Access and visualize properties calculated during simulations"
      ],
      "metadata": {
        "id": "pncDcGFmCPLW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Direct access to properties\n",
        "# calculated during the simulation\n",
        "\n",
        "# Plot heat capacity vs frequency\n",
        "freq_full = phonons.frequency.flatten(order='C')\n",
        "cv_1d = phonons.heat_capacity.flatten(order='C')[3:]\n",
        "\n",
        "print('\\n')\n",
        "plt.figure()\n",
        "plt.scatter(freq_full[3:],1e23*cv_1d,s=15)\n",
        "plt.ylabel (r\"$C_{v}$ ($10^{23}$ J/K)\", fontsize=25, fontweight='bold')\n",
        "plt.xlabel (\"$\\\\nu$ (THz)\", fontsize=25, fontweight='bold')\n",
        "plt.ylim(0.9*1e23*cv_1d[cv_1d>0].min(), \n",
        "         1.05*1e23*cv_1d.max())\n",
        "plt.show()\n",
        "\n",
        "# Plot phonon bandwidth vs frequency\n",
        "band_width_flatten = phonons.bandwidth.flatten(order='C')\n",
        "freq = freq_full[band_width_flatten!=0]\n",
        "\n",
        "print('\\n')\n",
        "plt.figure()\n",
        "plt.scatter(freq,band_width_flatten[band_width_flatten!=0] ,s=15)\n",
        "plt.ylabel (r\"$\\Gamma$ (THz)\", fontsize=25, fontweight='bold')\n",
        "plt.xlabel (\"$\\\\nu$ (THz)\", fontsize=25, fontweight='bold')\n",
        "plt.ylim(0.95*band_width_flatten .min(), 1.05*band_width_flatten .max())\n",
        "plt.show()\n",
        "\n",
        "# Plot phase space vs frequency\n",
        "print('\\n')\n",
        "plt.figure()\n",
        "plt.scatter(freq_full[3:],phonons.phase_space.flatten(order='C')[3:],s=15)\n",
        "plt.ylabel (\"Phase space\", fontsize=25, fontweight='bold')\n",
        "plt.xlabel (\"$\\\\nu$ (THz)\", fontsize=25, fontweight='bold')\n",
        "plt.ylim(phonons.phase_space.min(), phonons.phase_space.max())\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "Tx5FIPnKCPZ1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Calculate and visualize $\\kappa_{per \\ mode}$ and $\\kappa_{cum}$"
      ],
      "metadata": {
        "id": "kq8Fg_H1CQu-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def cumulative_cond_cal(freq,full_cond,n_phonons):\n",
        "\n",
        "  conductivity = np.einsum('maa->m', 1/3 * full_cond)\n",
        "  conductivity = conductivity.reshape(n_phonons)\n",
        "  cumulative_cond = np.zeros_like(conductivity)\n",
        "  freq_reshaped = freq.reshape(n_phonons)\n",
        "\n",
        "  for mu in range(cumulative_cond.size):\n",
        "      single_cumulative_cond = conductivity[(freq_reshaped < freq_reshaped[mu])].sum()\n",
        "      cumulative_cond[mu] = single_cumulative_cond\n",
        "  \n",
        "  return cumulative_cond"
      ],
      "metadata": {
        "id": "km6Vvh8lCTkt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Compute conductivity with per phonon mode basis using different methods\n",
        "kappa_rta_per_mode = np.einsum('maa->m',1/3*Conductivity(phonons=phonons, method='rta').conductivity)\n",
        "kappa_inv_per_mode = np.einsum('maa->m',1/3*Conductivity(phonons=phonons, method='inverse').conductivity)\n",
        "kappa_sc_per_mode = np.einsum('maa->m',1/3*Conductivity(phonons=phonons, method='sc',n_iterations=20).conductivity)\n",
        "\n",
        "# Compute cumulative conductivity by frequency using different methods\n",
        "kappa_rta_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='rta').conductivity,phonons.n_phonons)\n",
        "kappa_sc_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='sc',n_iterations=20).conductivity,phonons.n_phonons)\n",
        "kappa_inv_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='inverse').conductivity,phonons.n_phonons)\n",
        "kappa_qhgk_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='qhgk').conductivity,phonons.n_phonons)\n",
        "print('\\n')\n",
        "\n",
        "# Visualize the cumulative conductivity vs frequency\n",
        "plt.figure()\n",
        "plt.plot(freq_full,kappa_rta_per_mode,'r.',label='RTA')\n",
        "plt.plot(freq_full,kappa_sc_per_mode,'mo',label='Self Consistent',ms=8)\n",
        "plt.plot(freq_full,kappa_inv_per_mode,'k.',label='Direct Inversion')\n",
        "plt.xlabel (\"$\\\\nu$ (THz)\", fontsize=25, fontweight='bold')\n",
        "plt.ylabel(r'$\\kappa(W/m/K)$', fontsize=25, fontweight='bold')\n",
        "plt.legend(loc=1,frameon=False)\n",
        "#plt.grid()\n",
        "plt.show()\n",
        "print('\\n')\n",
        "\n",
        "# Visualize the cumulative conductivity vs frequency\n",
        "plt.figure()\n",
        "plt.plot(freq_full,kappa_rta_cum_freq,'r.',label='RTA')\n",
        "plt.plot(freq_full,kappa_sc_cum_freq,'mo',label='Self Consistent',ms=8)\n",
        "plt.plot(freq_full,kappa_inv_cum_freq,'k.',label='Direct Inversion')\n",
        "plt.xlabel (\"$\\\\nu$ (THz)\", fontsize=25, fontweight='bold')\n",
        "plt.ylabel(r'$\\kappa_{cum}(W/m/K)$', fontsize=25, fontweight='bold')\n",
        "plt.legend(loc=4,frameon=False)\n",
        "plt.grid()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "HZSC-vLJCVI0"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}