EXAMPLE SIM TEST.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "2d3829ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd\n",
    "import gsd.hoomd\n",
    "import numpy\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "import itertools\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "bf1207b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_particles=200\n",
    "\n",
    "spacing = 2\n",
    "K = math.ceil(N_particles**(1 / 3))\n",
    "L = K * spacing\n",
    "x = numpy.linspace(-L / 2, L / 2, K, endpoint=False)\n",
    "position = list(itertools.product(x, repeat=3))\n",
    "\n",
    "frame = gsd.hoomd.Frame()\n",
    "frame.particles.N = N_particles\n",
    "frame.particles.position = position[0:N_particles]\n",
    "frame.particles.typeid = [0] * N_particles\n",
    "frame.particles.types = ['A']\n",
    "frame.particles.mass=[5]*N_particles\n",
    "frame.configuration.box = [L, L, L, 0, 0, 0]\n",
    "a=frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ce53f4d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "ac52297e",
   "metadata": {},
   "outputs": [],
   "source": [
    "test = [np.random.uniform(-10,10,size=3) for _ in range(200)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "0637f8fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau=1\n",
    "kT=.1\n",
    "dt=.0005\n",
    "epsilon=2\n",
    "frames=3e5\n",
    "trigger=1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c78dd84b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "5df56397",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.8, 0. , 0. ],\n",
       "       [0. , 1.8, 0. ],\n",
       "       [0. , 0. , 1.8]])"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Script in order to account for Intertia, which hopefully solves the line 9 problem of providing no torques\n",
    "#The problem: adding the mass property to the particles breaks the line 4 command\n",
    "general_positions = numpy.array([[0, 0, 0], [1, 0, 0]])\n",
    "\n",
    "I_ref = numpy.array([[0.4, 0, 0],\n",
    "                   [0, 0.4, 0],\n",
    "                   [0, 0, 0.4]])\n",
    "I_general = numpy.zeros(shape=(3,3))\n",
    "for r in general_positions:\n",
    "    I_general += I_ref + 1 * (numpy.dot(r, r) * numpy.identity(3) - numpy.outer(r, r))\n",
    "    \n",
    "I_diagonal, E_vec = numpy.linalg.eig(I_general)\n",
    "\n",
    "R = E_vec.T\n",
    "\n",
    "diagonal_positions = numpy.dot(R, general_positions.T).T\n",
    "\n",
    "I_check = numpy.zeros(shape=(3,3))\n",
    "for r in diagonal_positions:\n",
    "    I_check += I_ref + 1 * (numpy.dot(r, r) * numpy.identity(3) - numpy.outer(r, r))\n",
    "\n",
    "#mass = 1\n",
    "#I = numpy.zeros(shape=(3,3))\n",
    "#for r in central_position:\n",
    "#    I += mass * (numpy.dot(r, r) * numpy.identity(3) - numpy.outer(r, r))\n",
    "I_check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "d1c64cf8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.8, 1.7999999999999998]\n"
     ]
    }
   ],
   "source": [
    "a.particles.moment_inertia = [I_general[0, 0], I_general[1, 1], I_general[2, 2]] * 200\n",
    "print(a.particles.moment_inertia[:2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "9cfe72d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "device = hoomd.device.auto_select()\n",
    "sim = hoomd.Simulation(device=device, seed=1)\n",
    "sim.create_state_from_snapshot(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "89a786ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "rigid = hoomd.md.constrain.Rigid()\n",
    "\n",
    "rigid.body['dimer'] = {\n",
    "    \"constituent_types\": ['A', 'A'],\n",
    "    \"positions\": [[0, 0, -0.5], [0, 0, 0.5]],\n",
    "    \"orientations\": [(1.0, 0.0, 0.0, 0.0), (1.0, 0, 0, 1.0)],\n",
    "    \"charges\": [0.0, 0.0],\n",
    "    \"diameters\": [1.0, 1.0]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "a77a76d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator = hoomd.md.Integrator(dt=dt, integrate_rotational_dof=True)\n",
    "\n",
    "cell=hoomd.md.nlist.Cell(buffer=0.4)\n",
    "\n",
    "lperp=0.3\n",
    "lpar=1.0\n",
    "sigmin=2*min(lperp,lpar)\n",
    "sigmax=2*max(lperp,lpar)\n",
    "\n",
    "gay_berne=hoomd.md.pair.aniso.GayBerne(nlist=cell, default_r_cut=2.5)\n",
    "gay_berne.params[('A','A')]=dict(epsilon=epsilon, lperp=lperp, lpar=lpar)\n",
    "gay_berne.r_cut[('A','A')]=2.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "d8fa0dcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.operations.integrator=integrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "9ce6df32",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator.rigid=rigid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d277ef",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "ff02cb02",
   "metadata": {},
   "outputs": [],
   "source": [
    "rigid_centers_and_free=hoomd.filter.Rigid((\"center\",\"free\"))\n",
    "nvt=hoomd.md.methods.NVT(filter=hoomd.filter.All(),kT=kT,tau=tau)\n",
    "#nve=hoomd.md.methods.NVE(filter=hoomd.filter.All())\n",
    "integrator.methods.append(nvt)\n",
    "integrator.forces.append(gay_berne)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "54a1306e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=kT)\n",
    "sim.run(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "1206e448",
   "metadata": {},
   "outputs": [],
   "source": [
    "logger=hoomd.logging.Logger(categories=['scalar'])\n",
    "logger+=gay_berne\n",
    "logger.add(gay_berne,quantities=['energy'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "fe1c64ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "Table = hoomd.write.Table(output=open('gay_berne_log_different_directions.txt', mode='w', newline='\\n'),\n",
    "                          trigger = hoomd.trigger.Periodic(1), logger=logger)\n",
    "sim.operations.writers.append(Table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "06744524",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.run(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "f1836148",
   "metadata": {},
   "outputs": [],
   "source": [
    "#g=sim.state.get_snapshot()\n",
    "#g.particles.orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "447e9f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "thermodynamic_properties=hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())\n",
    "gsd_writer=hoomd.write.GSD(filename='traj_different_directions1.gsd',trigger=hoomd.trigger.Periodic(trigger),mode='wb',filter=hoomd.filter.All())\n",
    "sim.operations.writers.append(gsd_writer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "1766dcb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ramp = hoomd.variant.Ramp(A=0, B=1, t_start=sim.timestep, t_ramp=20000)\n",
    "#initial_box = sim.state.box\n",
    "#final_box = sim.state.box\n",
    "#final_box.Lx *= 0.80\n",
    "#final_box.Ly *= 0.80\n",
    "#final_box.Lz *= 0.80\n",
    "#box_resize_trigger = hoomd.trigger.Periodic(trigger)\n",
    "#box_resize = hoomd.update.BoxResize(box1=initial_box,\n",
    "#                                    box2=final_box,\n",
    "#                                    variant=ramp,\n",
    "#                                    trigger=box_resize_trigger)\n",
    "#sim.operations.updaters.append(box_resize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "00b456d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.state.box.Lx *= 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "c7ca8449",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.operations.computes.append(thermodynamic_properties)\n",
    "sim.run(frames,write_at_start=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "5c9a2222",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300000"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "id": "05acb17f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2978.140596925596"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.tps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "67a937df",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = gsd.hoomd.open(\"traj_different_directions1.gsd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "26f756d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "301"
      ]
     },
     "execution_count": 121,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b491645d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec9fb1a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdb6d93a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "551c9fcc",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89a9704a",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}