{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align: center\"><a href=\"../index.ipynb\">LAMMPS Python Tutorials</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 5: Monte Carlo Relaxation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random, math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup perfect system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lammps import lammps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = lammps()\n",
    "cmd = L.cmd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd.units(\"lj\")\n",
    "cmd.atom_style(\"atomic\")\n",
    "cmd.atom_modify(\"map array sort\", 0, 0.0)\n",
    "\n",
    "cmd.dimension(2)\n",
    "\n",
    "cmd.lattice(\"hex\", 1.0)\n",
    "cmd.region(\"box block\", 0, 10, 0, 5, -0.5, 0.5)\n",
    "\n",
    "cmd.create_box(1, \"box\")\n",
    "cmd.create_atoms(1, \"box\")\n",
    "cmd.mass(1, 1.0)\n",
    "\n",
    "cmd.pair_style(\"lj/cut\", 2.5)\n",
    "cmd.pair_coeff(1, 1, 1.0, 1.0, 2.5)\n",
    "cmd.pair_modify(\"shift\", \"yes\")\n",
    "\n",
    "cmd.neighbor(0.3, \"bin\")\n",
    "cmd.neigh_modify(\"delay\", 0, \"every\", 1, \"check\", \"yes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L.ipython.image(zoom=1.6,size=[320,320])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd.run(0, \"post\", \"no\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emin = L.get_thermo(\"pe\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd.dump(\"3 all movie 25 movie.mp4 type type zoom 1.6 adiam 1.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Disorder system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(27848)\n",
    "deltaperturb = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = L.numpy.extract_atom(\"x\")\n",
    "for i in range(len(pos)):\n",
    "    x, y = pos[i][0], pos[i][1]\n",
    "    dx = deltaperturb * random.uniform(-1, 1)\n",
    "    dy = deltaperturb * random.uniform(-1, 1)\n",
    "    pos[i] = (x+dx, y+dy, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd.run(0, \"post\", \"no\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L.ipython.image(zoom=1.6,size=[320,320])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimize using Monte Carlo moves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estart = L.get_thermo(\"pe\")\n",
    "elast = estart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "naccept = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = [estart]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "niterations = 3000\n",
    "deltamove = 0.1\n",
    "kT = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "natoms = L.extract_global(\"natoms\")\n",
    "\n",
    "for i in range(niterations):\n",
    "    pos = L.numpy.extract_atom(\"x\")\n",
    "    iatom = random.randrange(0, natoms)\n",
    "    current_atom = pos[iatom]\n",
    "    \n",
    "    x0, y0 = current_atom[0], current_atom[1]\n",
    "    \n",
    "    dx = deltamove * random.uniform(-1, 1)\n",
    "    dy = deltamove * random.uniform(-1, 1)\n",
    "    \n",
    "    pos[iatom] = (x0+dx, y0+dy, 0)\n",
    "    \n",
    "    cmd.run(1, \"pre no post no\")\n",
    "    \n",
    "    e = L.get_thermo(\"pe\")\n",
    "    energies.append(e)\n",
    "    \n",
    "    if e <= elast:\n",
    "        naccept += 1\n",
    "        elast = e\n",
    "    elif random.random() <= math.exp(natoms*(elast-e)/kT):\n",
    "        naccept += 1\n",
    "        elast = e\n",
    "    else:\n",
    "        pos[iatom] = (x0, y0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xlabel('iteration')\n",
    "plt.ylabel('potential energy')\n",
    "plt.plot(energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L.get_thermo(\"pe\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "naccept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L.ipython.image(zoom=1.6, size=[320,320])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# close dump file to access it\n",
    "cmd.undump(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L.ipython.video(\"movie.mp4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}