Attempting_to_create_the_desired_two_particle_environment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "29600823",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hoomd\n",
    "import gsd.hoomd\n",
    "import numpy\n",
    "import math\n",
    "import itertools "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4eb4c297",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "fn = os.path.join(os.getcwd(), 'dimer_centers.gsd')\n",
    "![ -e \"$fn\" ] && rm \"$fn\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0daedfe7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#dimer_positions = [[-0.5, 0, 0], [0.5, 0, 0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "81d06902",
   "metadata": {},
   "outputs": [],
   "source": [
    "central_rotation = 0.9\n",
    "central_position = [2, 1, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b3a0c53b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/noah/miniconda3/envs/try2/lib/python3.10/site-packages/gsd/hoomd.py:420: FutureWarning: Snapshot is deprecated, use Frame\n",
      "  warnings.warn(\"Snapshot is deprecated, use Frame\", FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "s = gsd.hoomd.Snapshot()\n",
    "s.particles.N = 2\n",
    "s.particles.types = ['A']\n",
    "s.particles.typeid = [0,0]\n",
    "s.particles.position = [[1,1,1], [2,1,2]]\n",
    "s.particles.orientation = [[1, 0, 0, 0], [0, 0.70710678118, 0.0, 0.70710678118]]\n",
    "s.configuration.box = [8, 8, 8, 0, 0, 0]\n",
    "s.particles.mass = [2] * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ffe7e67d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0., -5., -5.],\n",
       "       [-5.,  0., -5.],\n",
       "       [-5., -5.,  0.]])"
      ]
     },
     "execution_count": 6,
     "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e162566a",
   "metadata": {},
   "outputs": [],
   "source": [
    "s.particles.moment_intertia = [I[0, 0], I[1, 1], I[2, 2]] * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1cfd0413",
   "metadata": {},
   "outputs": [],
   "source": [
    "cpu = hoomd.device.CPU()\n",
    "sim = hoomd.Simulation(device=cpu, seed=1)\n",
    "sim.create_state_from_snapshot(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2447cced",
   "metadata": {},
   "outputs": [],
   "source": [
    "rigid = hoomd.md.constrain.Rigid()\n",
    "\n",
    "#rigid.body['dimer'] = {\n",
    "#    \"constituent_types\": ['A', 'A'],\n",
    "#    \"positions\": [[-0.5, 0, 0], [0.5, 0, 0]],\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": 10,
   "id": "04e55ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator = hoomd.md.Integrator(dt=0.005, integrate_rotational_dof=True)\n",
    "\n",
    "cell = hoomd.md.nlist.Cell(buffer=0.4)\n",
    "\n",
    "gay_berne = hoomd.md.pair.aniso.GayBerne(nlist=cell, default_r_cut=2.5)\n",
    "gay_berne.params[('A', 'A')] = dict(epsilon=1.0, lperp=0.2, lpar=0.7)\n",
    "gay_berne.r_cut[('A', 'A')] = 2.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d4d241fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.operations.integrator = integrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "13fa3bca",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator.rigid = rigid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bc1165d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/noah/miniconda3/envs/try2/lib/python3.10/site-packages/hoomd/md/methods/methods.py:102: FutureWarning: NVT is deprecated and wil be removed in hoomd 4.0. In version 4.0, use the ConstantVolume method with the desired thermostat from hoomd.md.methods.thermostats.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1.5)\n",
    "\n",
    "kT = 1.5\n",
    "rigid_centers_and_free = hoomd.filter.Rigid((\"center\", \"free\"))\n",
    "nvt = hoomd.md.methods.NVT(kT=1.5, filter=hoomd.filter.All(), tau=1.0)\n",
    "integrator.methods.append(nvt)\n",
    "integrator.forces.append(gay_berne)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "91ec2d42",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.state.thermalize_particle_momenta(filter=rigid_centers_and_free, kT=kT)\n",
    "sim.run(0)\n",
    "nvt.thermalize_thermostat_dof()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a6dbdb83",
   "metadata": {},
   "outputs": [],
   "source": [
    "thermodynamic_properties = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())\n",
    "gsd_writer = hoomd.write.GSD(filename = 'traj.gsd', trigger = hoomd.trigger.Periodic(10), mode = 'wb', filter = hoomd.filter.All())\n",
    "sim.operations.writers.append(gsd_writer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31807383",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "eec805c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.operations.computes.append(thermodynamic_properties)\n",
    "sim.run(5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1c3f1d96",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sim.state.get_snapshot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "2a575d35",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        , 0.        , 0.        ],\n",
       "       [0.        , 0.70710678, 0.        , 0.70710678]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.particles.orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "3679dd72",
   "metadata": {},
   "outputs": [],
   "source": [
    "thermodynamic_quantities = hoomd.md.compute.ThermodynamicQuantities(filter=hoomd.filter.All())\n",
    "sim.operations.computes.append(thermodynamic_quantities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "6f97a852",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Translational degrees of freedon: 6.0\n"
     ]
    }
   ],
   "source": [
    "translational_degrees_of_freedom = thermodynamic_quantities.translational_degrees_of_freedom\n",
    "print('Translational degrees of freedon:', translational_degrees_of_freedom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b16d5a93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thermodynamic_quantities.rotational_kinetic_energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "93562c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Current State: Everything \"works\" The only problem is that the particles are not rotating, and don't have any\n",
    "# rotational kinetic energy, the source of the problem is likely found 06 01 which provides the initial system,\n",
    "# which I believe provides the rotational kinetic energy, thus before writing the state I need to find that"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b213c451",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}