{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "63604ef2",
   "metadata": {},
   "source": [
    "# Custom Symmetry Group Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff87f25d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import symd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "246c1078",
   "metadata": {},
   "source": [
    "### Make custom group\n",
    "\n",
    "Rotational group with reflection boundary conditions. Asymmetric unit is just a guess - not actually correct, so we may have overlapping particles and need to re-run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bdbb677",
   "metadata": {},
   "outputs": [],
   "source": [
    "rot = 7\n",
    "genpos = []\n",
    "for ri in range(rot):\n",
    "    c = np.round(np.cos(ri * np.pi * 2 / rot), 4)\n",
    "    s = np.round(np.sin(ri * np.pi * 2 / rot), 4)\n",
    "    # convert numbers to affine notation\n",
    "    genpos.append(\n",
    "        f\"{c}x {-s:+}y {-0.5 * c + 0.5 * s + 0.5:+},{s}x {c:+}y {-0.5 * s - 0.5 * c + 0.5:+}\"\n",
    "    )\n",
    "specpos = [symd.Group(\"Oblique\", [\"1/2, 1/2\"], specpos=[])]\n",
    "asymm = f\"1/2≤x≤1.0;1/2≤y≤1.0\"\n",
    "my_group = symd.Group(\"Oblique\", genpos, specpos, \"quasi\", asymm)\n",
    "print(my_group)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cfe53fc",
   "metadata": {},
   "source": [
    "### Run\n",
    "\n",
    "Run the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fab73842",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_sim(\n",
    "    n, number_density, group, images, w=None, retries=50, pos_frames=0, steps=10**6\n",
    "):\n",
    "    for i in range(retries):\n",
    "        print(\"Trying on \", i)\n",
    "        try:\n",
    "            np.random.seed(i)\n",
    "            cell = symd.groups.get_cell(number_density, group, 2, n, w)\n",
    "            md = symd.Symd(\n",
    "                nparticles=n,\n",
    "                cell=cell,\n",
    "                ndims=2,\n",
    "                images=images,\n",
    "                force=\"lj\",\n",
    "                wyckoffs=w,\n",
    "                group=group,\n",
    "                steps=steps,\n",
    "                exeDir=\"quasi\",\n",
    "                start_temperature=0.5,\n",
    "                temperature=0.1,\n",
    "                pressure=0.25,\n",
    "            )\n",
    "            md.remove_overlap()\n",
    "            md.runParams[\"box_update_period\"] = 10\n",
    "            md.runParams[\"langevin_gamma\"] = 0.5\n",
    "            md.log_positions(frames=pos_frames)\n",
    "            try:\n",
    "                md.run()\n",
    "            except RuntimeError as e:\n",
    "                d = md.number_density()\n",
    "                if d < 0.6:\n",
    "                    print(\"Not dense enough, retrying\", d)\n",
    "                    continue\n",
    "\n",
    "            # Basically E-min\n",
    "            md.runParams[\"start_temperature\"] = 0.1\n",
    "            md.runParams[\"temperature\"] = 1e-2\n",
    "            md.runParams[\"langevin_gamma\"] = 0.5\n",
    "            md.runParams[\"Pressure\"] = None\n",
    "            md.runParams[\"box_update_period\"] = 0\n",
    "            md.runParams[\"steps\"] = steps // 10\n",
    "            if pos_frames > 0:\n",
    "                md.log_positions(filename=\"equil.xyz\", frames=pos_frames // 10)\n",
    "            try:\n",
    "                md.run()\n",
    "            except RuntimeError as e:\n",
    "                continue\n",
    "            config = md.positions[-1]\n",
    "            break\n",
    "        except RuntimeError as e:\n",
    "            print(e)\n",
    "            md = None\n",
    "    return md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b330a5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "md = run_sim(256, 0.1, my_group, [0, 0], pos_frames=1000, w=[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d84856",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(md.positions[-1, :, 0], md.positions[-1, :, 1], \".\")"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}