Syncing codebases to update nucleosome wrapping abilities, internucleosome attraction calculation/scaling, steric checking, et al.

analysis/brownian_dynamics/00_Brownian_Dynamics_of_Rouse_Polymer.ipynb

@@ -0,0 +1,264 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Create Brownian Dyanmics Simulation of a Rouse Polymer"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import seaborn as sns\n",
+    "import numba\n",
+    "import matplotlib.pyplot as plt\n",
+    "import os"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import wlcsim\n",
+    "from wlcsim.bd import rouse, homolog\n",
+    "from wlcsim.bd import runge_kutta"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = int(1e2+1)\n",
+    "L = 17.475\n",
+    "b = 0.015\n",
+    "D = 2e1\n",
+    "R = 1\n",
+    "dt = rouse.recommended_dt(N, L, b, D)\n",
+    "Aex = 0; # no confinement\n",
+    "N_tot, loop_list = homolog.points_to_loops_list(N, []) # no loops\n",
+    "# tether_list = np.array([]).astype(int) # no tethered loci"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2.5196008663366337e-05"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "dt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "'/ru-auth/local/home/risc_soft/miniconda3/envs/bd2/lib/python3.7/site-packages/bruno_util/__init__.py'"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "import bruno_util\n",
+    "bruno_util.__file__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Nhat = L/b; L0 = L/(N-1); Dhat = D*N/Nhat; bhat = np.sqrt(L0*b)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Nt = 1e6; # total number of equi-space time steps\n",
+    "Nlin = 1e4; # see docstring of loglinsample\n",
+    "t = np.arange(0, Nt*dt, dt) # take Nt time steps\n",
+    "from bruno_util import numpy as bnp\n",
+    "t_i, i_i = bnp.loglinsample(Nt, Nlin, 0.6)\n",
+    "t_save = t[t_i]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(1000001,)"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.shape(t)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i_i = [np.round(i).astype(int) for i in i_i] # numba workaround"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Run Simulation with with_integrator\n",
+    "# x = rouse.with_integrator(runge_kutta.srk1/runge_kutta.lena4)\n",
+    "bead_coords= rouse.with_integrator(N, L, b, D, t, integrator = runge_kutta.rk4_thermal_lena)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note the dimensions of `bead_coords` are `len(t)` by `N`(number of beads) by 3 (x,y,z). `bead_coord` contains the x,y,z-coordinates for each of the 101 beads at each timestep "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(1000001, 101, 3)"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "bead_coords.shape"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Save bead_coord in your local directory:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "local_directory = '/ru-auth/local/home/abrenner/myscratch/bd/bead_coords/'\n",
+    "file_name = '2021_03_09_rouse_rk4lena.npy'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "/rugpfs/fs0/risc_lab/scratch/abrenner/bd/bead_coords\n"
+     ]
+    }
+   ],
+   "source": [
+    "cd $local_directory"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Subset to the first 100 snapshots"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "subset_bead_coords = bead_coords[:100,:,:]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 31,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.save(file_name, subset_bead_coords)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}