diff --git a/05_msd_analysis_and_example/msd_analysis.ipynb b/05_msd_analysis_and_example/msd_analysis.ipynb new file mode 100644 index 0000000..ba27c6f --- /dev/null +++ b/05_msd_analysis_and_example/msd_analysis.ipynb @@ -0,0 +1,425 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calculate MSDs" + ] + }, + { + "cell_type": "code", + "execution_count": 425, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Populating the interactive namespace from numpy and matplotlib\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/research/samuel.cohen/.conda/envs/pappuPython/lib/python3.9/site-packages/IPython/core/magics/pylab.py:159: UserWarning: pylab import has clobbered these variables: ['matrix']\n", + "`%matplotlib` prevents importing * from pylab and numpy\n", + " warn(\"pylab import has clobbered these variables: %s\" % clobbered +\n" + ] + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import seaborn as sns\n", + "%pylab inline\n", + "sns.set()\n", + "sns.set_style(\"white\")\n", + "sns.set_style(\"ticks\")\n", + "sns.set_context(\"poster\")\n", + "pylab.rcParams['font.family'] = 'STIXGeneral'\n", + "pylab.rcParams['mathtext.fontset'] = 'stix'\n", + "pylab.rcParams['xtick.labelsize'] = 24\n", + "pylab.rcParams['ytick.labelsize'] = 24\n", + "import sys\n", + "myList = [[]]\n", + "from importlib import reload\n", + "import os\n", + "import subprocess as sproc\n", + "import pickle\n", + "plEKW = {'lw':2, 'alpha':0.7}\n", + "erEKW = {'fmt':'.-', 'lw':2, 'capsize':2, 'capthick':2, 'elinewidth':2, 'alpha':0.7}" + ] + }, + { + "cell_type": "code", + "execution_count": 426, + "metadata": {}, + "outputs": [], + "source": [ + "import MDAnalysis as mdan\n", + "from MDAnalysis.analysis import rdf as mAnRDF\n", + "from MDAnalysis.analysis import distances as mAnDist\n", + "import networkx as nx\n", + "sns.set_style(\"white\")\n", + "sns.set_style(\"ticks\")\n", + "sns.set_context(\"paper\")" + ] + }, + { + "cell_type": "code", + "execution_count": 427, + "metadata": {}, + "outputs": [], + "source": [ + "def gen_random_residue_selection(num_of_sample, tot_mols=1728, mol_cycle=16):\n", + " dum_res_list = np.arange(0, tot_mols+1)\n", + " dum_res_list = dum_res_list[np.where(dum_res_list % mol_cycle != 0)[0]]+1\n", + " dum_res_list = np.random.choice(dum_res_list, num_of_sample, replace=False)\n", + " dum_res_list.sort()\n", + " dum_res_list = [\"resid \" + str(a_res) for a_res in dum_res_list]\n", + " dum_res_list = \" or \".join(dum_res_list)\n", + " return dum_res_list\n", + "\n", + "def Read_N130_Trajectory(top_file, trj_file, num_of_reps=108, verbose=False):\n", + " dum_universe = mdan.Universe(top_file, trj_file, format=\"LAMMPS\")\n", + " if verbose:\n", + " print(\"The simulation box is\", dum_universe.dimensions[:3])\n", + " #Adding the names of the molecules\n", + " names_of_res = ['N130']\n", + " [names_of_res.append('rpL5') for i in range(15)]\n", + " names_of_res = names_of_res*num_of_reps\n", + " dum_universe.add_TopologyAttr(topologyattr='resnames', values=names_of_res)\n", + " \n", + " #Getting the lengths of the different molecules (or bits) we care about\n", + " n130_len = len(dum_universe.residues[0].atoms)\n", + " pept_len = len(dum_universe.residues[1].atoms)\n", + " n130_arm_len = int((n130_len - 1)/5)\n", + " mol_lens = [n130_len, n130_arm_len, pept_len]\n", + " if verbose:\n", + " print(\"N130 has {:} beads; each arm has {:} beads.\".format(n130_len, n130_arm_len))\n", + " print(\"The peptide has {:} beads.\".format(pept_len))\n", + " pept_num = int(len(dum_universe.select_atoms('resname rpL5'))/pept_len)\n", + " n130_num = 108\n", + " mol_nums = [n130_num, pept_num]\n", + " return [dum_universe, mol_lens, mol_nums]" + ] + }, + { + "cell_type": "code", + "execution_count": 428, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The simulation box is [338.5646 338.5646 338.5646]\n", + "N130 has 211 beads; each arm has 42 beads.\n", + "The peptide has 17 beads.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/research/samuel.cohen/.conda/envs/pappuPython/lib/python3.9/site-packages/MDAnalysis/coordinates/DCD.py:165: DeprecationWarning: DCDReader currently makes independent timesteps by copying self.ts while other readers update self.ts inplace. This behavior will be changed in 3.0 to be the same as other readers. Read more at https://github.com/MDAnalysis/mdanalysis/issues/3889 to learn if this change in behavior might affect you.\n", + " warnings.warn(\"DCDReader currently makes independent timesteps\"\n" + ] + } + ], + "source": [ + "# Run this for each replica\n", + "dum_dat = 'lk_8_Init_0.data'\n", + "dum_trj = '/project/fava/work/furqan.dar/2020/1_N130/5_NVT/1_AllSys/lk_8/Run_5/lk_8_5_postNPT_FULL.dcd'\n", + "dum_un, dum_lens, dum_nums = Read_N130_Trajectory(dum_dat, dum_trj, 108, True)" + ] + }, + { + "cell_type": "code", + "execution_count": 429, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000\n" + ] + } + ], + "source": [ + "from MDAnalysis.coordinates.memory import MemoryReader\n", + "\n", + "sel = dum_un.trajectory[1001:]\n", + "positions = np.array([ts.positions for ts in sel])\n", + "new_universe = mdan.Universe(dum_dat, positions, format=MemoryReader)\n", + "print(new_universe.trajectory.n_frames)" + ] + }, + { + "cell_type": "code", + "execution_count": 430, + "metadata": {}, + "outputs": [], + "source": [ + "#Only the charged residues from each block\n", + "\n", + "A0_LIST = [6, 7, 10, 12, 17]\n", + "A1_LIST = [20, 22, 24, 25, 27]\n", + "A2_LIST = [33, 34, 35, 37, 39, 40, 41, 42, 43]\n", + "PE_LIST = [44, 45, 46, 47, 48, 50, 52, 56, 57, 58]" + ] + }, + { + "cell_type": "code", + "execution_count": 431, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "108\n" + ] + } + ], + "source": [ + "def _get_only_FDs(mdan_universe):\n", + " \"\"\"\n", + " Given this MDAnalysis Universe object for the N130 CG simulations, we extract only the Folded Domains (FD's)\n", + " of the N130 molecules. This can easily be modified, however.\n", + " \"\"\"\n", + " \n", + " return mdan_universe.select_atoms(\"type 1\")\n", + "\n", + "dum_sel = _get_only_FDs(new_universe)\n", + "print(dum_sel.n_atoms)\n", + "\n", + "def _get_distance_map(mdan_universe, atom_sel):\n", + " \"\"\"\n", + " Given a frame from an MDAnalysis Universe object, we calculate the distance map for the given atom_selection.\n", + " \"\"\"\n", + " \n", + " dist_map = mAnDist.distance_array(atom_sel.positions, atom_sel.positions,\n", + " box=mdan_universe.dimensions)\n", + " \n", + " return dist_map\n", + "\n", + "dum_map = _get_distance_map(new_universe, dum_sel)\n", + "\n", + "def _filter_distance_map(this_map, r_lo=50, r_hi=100):\n", + " \"\"\"\n", + " Given this Distance Map, we create a sort of adjacency matrix where only the distances between r_lo and r_hi\n", + " define an unweighted edge.\n", + " \"\"\"\n", + " \n", + " map_copy = this_map[:]\n", + " map_copy[ map_copy > r_hi] = 0\n", + " map_copy[ map_copy < r_lo] = 0\n", + " map_copy[map_copy != 0] = 1\n", + " \n", + " return np.array(map_copy, int)\n", + "\n", + "dum_map = _filter_distance_map(dum_map)" + ] + }, + { + "cell_type": "code", + "execution_count": 432, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate windowed MSD by averaging over all possible lag times\n", + "\n", + "# Adapts code from https://colab.research.google.com/github/kaityo256/zenn-content/blob/main/articles/msd_fft_python/msd_fft_python.ipynb#scrollTo=CuYbCG8DcQfH\n", + "\n", + "def calc_msd_np(x,y,z):\n", + " n = len(x)\n", + " msd = []\n", + " for s in range(1,n//4):\n", + " dx = x[s:] - x[:-s]\n", + " dy = y[s:] - y[:-s]\n", + " dz = z[s:] - z[:-s]\n", + " r2 = dx**2 + dy**2 + dz**2\n", + " msd.append(np.average(r2))\n", + " return msd" + ] + }, + { + "cell_type": "code", + "execution_count": 433, + "metadata": {}, + "outputs": [], + "source": [ + "matrix = np.array([[[dum_un.dimensions[:1]]]])\n", + "L = matrix.item()\n", + "def adjust_periodic(x):\n", + " for i in range(len(x)-1):\n", + " if x[i+1] - x[i] > L/2:\n", + " x[i+1] -= (x[i+1] - x[i]+L/2)//L*L\n", + " if x[i+1] - x[i] < -L/2:\n", + " x[i+1] += (x[i] - x[i+1]+L/2)//L*L" + ] + }, + { + "cell_type": "code", + "execution_count": 434, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"For MSD analysis\"\"\"\n", + "\"\"\"Extract coordinates\"\"\"\n", + "\n", + "x = []\n", + "y = []\n", + "z = []\n", + "\n", + "for a_frame in new_universe.trajectory: # for each snapshot\n", + " \n", + " FD = new_universe.select_atoms(\"type 1\")\n", + "\n", + " # Positions of all molecules at each snapshot\n", + " r = FD.positions\n", + " r = list(r)\n", + "\n", + " x_new = [coordinates[0] for coordinates in r]\n", + " y_new = [coordinates[1] for coordinates in r]\n", + " z_new = [coordinates[2] for coordinates in r]\n", + " \n", + " # x, y, z coordinates for all molecules over the entire trajectory\n", + " x.append(x_new)\n", + " y.append(y_new)\n", + " z.append(z_new)" + ] + }, + { + "cell_type": "code", + "execution_count": 435, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"Extract coordinates of *each* molecule over entire trajectory\"\"\"\n", + "\n", + "x_newest = [] \n", + "y_newest = [] \n", + "z_newest = [] \n", + "\n", + "msd_tot = [] \n", + "\n", + "for i in range(len(FD)): # for each N130 domain\n", + " x_new = [coordinates[i] for coordinates in x]\n", + " y_new = [coordinates[i] for coordinates in y]\n", + " z_new = [coordinates[i] for coordinates in z]\n", + " \n", + " # Combine x, y, z coordinates over entire trajectory\n", + " x_newest.append(x_new)\n", + " y_newest.append(y_new)\n", + " z_newest.append(z_new)" + ] + }, + { + "cell_type": "code", + "execution_count": 436, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"Calculate windowed, time-averaged MSD for each molecule\"\"\" \n", + "msd_tot = [] \n", + "for i in range(len(FD)):\n", + " \n", + " # Convert list to array\n", + " x_dum = np.array(x_newest[i])\n", + " y_dum = np.array(y_newest[i])\n", + " z_dum = np.array(z_newest[i])\n", + "\n", + " adjust_periodic(x_dum)\n", + " adjust_periodic(y_dum)\n", + " adjust_periodic(z_dum)\n", + " \n", + " msd_np = calc_msd_np(x_dum,y_dum,z_dum)\n", + " # print(msd_np) # uncomment to output individual MSDs\n", + " msd_tot.append(msd_np)" + ] + }, + { + "cell_type": "code", + "execution_count": 437, + "metadata": {}, + "outputs": [], + "source": [ + "msd_avg = np.average(msd_tot, axis=0) # average MSD for all molecules" + ] + }, + { + "cell_type": "code", + "execution_count": 438, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAak0lEQVR4nO3deXhTZb4H8G/TlKW06d6ytaWFNm0hFLBQkV1QdFC2GbQjXBkdcGYcM4LOHerV63N9rlwQrUJZROcqKiCtTEGUQXQQBxDRUqALTRe6JF2hG226kjQ59w+UO8jakPScnHw//+W0SX/lhS/v865ugiAIICIi2VCIXQAREdkXg52ISGYY7EREMsNgJyKSGQY7EZHMKMUuQK1Wi10CEZFTKioquu5z0YMduHFxRER0fTfrFHMohohIZhjsREQyI8pQTHp6OtLT08X40UREsucm9pECarWaY+xERD10s+zkUAwRkcww2ImIZIbBTkQkMwx2IqJe1nWpG6VVzTB3Wxzy+ZLYoEREJFcWq4Dy6hbklTbgXGUzyqpbUNPQBoWbG/779/dAMzzQ7j+TwU5EZEdWqwDDeSPyShqQW9KAs2WN6OwyY/hQX6jD/fCre0cgcogvQkO84KF0d0gNDHYiojvwU5CfLW3E2bIGnC1tRGuHCRGDfKAZEYj7E8MxMjIAA/p79FpNDHYioh7otlhRVt1yJch15U3o6DIjYpAP4iL98cyieIwaHghvzz6i1chgJyK6CZPZgnOVzVd644X6Jpi6rYga6ouRkQF4cOIwxEYEwKsXe+S3wmAnIvoXXZe6UWho+rFH3ojiiosAgOgwP4waHoBfzhgBdbg/+veVbnxKtzIiol7Q1mmGrrwR+aWNyC9rRElVMzyUCsQM88dYdRD+7cFYRIX6oo+HYyY6HYHBTkQupbn1EvLLL4d4fmkjymtb4NnPAyMjAnDP6MF4aoEGkUN8oHR33m0+DHYiki1BEFDb2I6C8iYU6JugK29E5YU2+Hr1xcjIAMyaEIZRwwMQNlAFd4Wb2OXaDYOdiGTD3G1FaXXzlSAvKG9Cc9slDAzwROwwfzw8ZThGRQZgaLAX3NzkE+Q/x2AnIqfV2mG6EuAF+iacq7gIi1XA8KE+iB0WgGnjhiJ2mD/8Vf3ELrVXMdiJyCkIgoDahnbofuqN6y8Pqwzo74HYYf4Ypw7G4gdiEBXqi359XDvaXPu3JyLJ6rZYUVLVDF3Z5RAv0Dehpc2EQQEDEBvhj7lThiM2wh+hwd5QyGh83B4Y7EQkCeZuC4ormnG29PL5KoX6Jpi7rRgx1BexEf6YcVcoYof5w8/FhlVsYddg37dvHzIzM9He3o5169ahTx/xttQSkbRZrQLKalpwpqgO2cX1KNQ3QcDVG4Fiwv3RT8IbgaTqln9iJpMJjz/+OEpKSvDqq6/igQceAACkpaXh008/hVKpxOrVqxEeHo558+Zh3rx5eO+991BVVYXIyEiH/wJE5DwuGrtwprgOpwvrkX2uDq3tJkSH+WGsOhi/vl+N6DA/p9oIJFW3DHalUonU1FSkp6dfedbc3IyMjAykpaVBp9MhJSUFqampAIC2tjYYDAZEREQ4rmoicgomswW68kacKarH6aI66GuNCPLrj3HqYPxhYTziowLhJeJhWXJ1y2BXKBQIDg6+6llOTg4SExPh7u4OjUYDvV4P4HLvft26dVixYsV114hu3LgRmzZtsk/lRCQ5giCgqq4Np4vqcKaoDnmljXBzAzTDA3FfYhjGRgfLfg25FNg0eGU0GqFSqa68FgQBALB582ZUV1dj/fr1WLZsGcLCwq56n1arhVarveqZWq22pQQikoi2TjOyi+twurAOZ4rr0dDcicjBPhirDsKC6SMQF+HvsAsl6PpsCnaVSoXi4uIrrxWKy2cqrFy50j5VEZFkCYKAigutyNJdwMmCCyjQN0Hl2Qdj1EF4/BexGBMVxJUrIrMp2OPj47FlyxZYLBYUFhYiPDy8R+9PT0+/asyeiKTtktmC3HP1yCq4gKyCC6i72IkRob4YHxuCJx8eiRFDfbmWXELchJ/GUW5Cq9VCp9PB09MTkyZNQnJyMnbt2oV9+/ZdtSrGFmq1GkVFRTa9l4gcp7GlEz/kn8dJ3QXknquHu7sCY9VBGB8bgrtiQtgrF9nNsvO2gt2RGOxE0lFT34bv8mrxfV4tiiouYkjQAIyPG4iE2BDERQTAQ+m8R9nKzc2ykyv/iVyYIAgoq27BibO1OJFXi4rzrRgx1Ad3awbhT4+OQWiIN1ewOCFRgp1j7ETiEQQB5TVGHD1ThW9zalB/sQOxEQGYnRiOu0cNQrC/p9gl0h3iUAyRi6iqa8XRM9U4eqYatQ1t0IwIxJQxQ3D3qEHw8eordnnUQxyKIXJRdU0dOJZ9OczLaloQE+6HOZMiMDl+MCc/ZYzBTiQzHV1mfJtTg8NZlcgva0TEYBWmjh2K/3hiAkI4zOISOMZOJAMWq4Ccc/U4fLISJ87WwtvTAzPuCsXTvxyNsIGqW38AyQrH2ImcWOWFVhzOqsQ3pyrR2mHGPZpBuDchFKOjgmR1OTNdi2PsRDJyyWzB8ZwaHDyhR4G+CSMjA7B4dgwmxQ+GZz8PscsjCWCwEzkJw3kjvvzegMNZlfBwV2DWhDCs/PU4DAocIHZpJDEMdiIJM5ktOJZdjS+/N6DQ0IT4qCBoF43BhJEDuQuUboiTp0QS1NjSib8fL8eX3xugcHPDrAlheO6xcRgYwN453RonT4kkpNDQhM+PluF4bg0ihvhg7pRITI4fwt45XYOTp0QSZu624nhONT7/tgwlVS24RzMIa56ejJhhfjynhWzCYCcSSXunGQe+K8f+b8tg7hbwwMRwJD8+AUF+/cUujZwcg52ol100dmHf0VJ8cUIPX6++eGx2DKbfFYq+Hrw+juyDwU7US843tmPPP0twKLMCoSHe0D4yBhM1g7mRiOyOq2KIHKzivBGfHDqHYznVGBkRgJeeTMTY6CCOn5PDcFUMkYPU1Lfh4y+LcCy7CgmxA7FoVhRiwv3FLotkgqtiiHrRhaYOpP+jCF9nVWJsdBDeeHYqokL9xC6LXAiDnchOGls6kX6oGP/4wYC4iACsfXoyYiPYQ6fex2AnukPNrZew+3AxvvhOjxFDffHKUxMxekSQ2GWRC2OwE9moy9SNfUdKkfHNOQwJ9sZLTyRirJqToiQ+BjtRD1msAr7JqsSOgwVQuiugfWQsJscPZqCTZHC5I1EPZBfX4f3P81F/sROP3heNOZMi4KHkxiKSFi53JLoN1fVt+N99Z5FdXIeHJkfikVnR8PbsI3ZZ5MK43JHIRh1dZnxyqBj7jpYiITYEW/4ykxdbkOQx2ImuQxAEHDldhW378+HZzwP/+du7MU4dLHZZRLeFwU70M6VVzXhnbx70tUb8+n41HpocyfPQyakw2Il+1N5pxkcHdDj4vQHTxw3FC0vHw0/VT+yyiHqMwU4uTxAEHM+twbt78+Dj1Zc7RsnpMdjJpV1o6sDWPbnILWnA4tlqzJ06HEp3DruQc2Owk0uyWKzYd7QMH39ViJGRAdj87zN4UTTJBjcokcsprriIzbtz0NTahWcfGYvJY7hrlOSFG5TIZXSZurH9iwLs/7Yc9yeGY+mcOHj19xC7LCKbcIMSubz8skZsSDsDhQL4nz9MwsjIALFLInIYBjvJWpepG9sPFGD/8XLMnRKJJQ/G8tJokj0GO8nWv/bSuYSRXAmDnWSHvXRydQx2khX20okY7CQTJrMFHx7QYf+37KUTMdjJ6ZXXtOCNnadgMluw5ulJiIvgihdybQx2clpWq4DPjpXiw78XYPq4oVg+fxQ8+3FdOhGDnZxSY0sn1u86g9LqZvx5yV2YNHqw2CURSQaDnZzO8dwabN6djeFDfLHxzzMQ4NNf7JKIJIXBTk6jo8uMdz/Nw5HT1Vg6Jw5zp0RCoeAZL0Q/x0PAyCkU6puQ8vEp9PVwx5srpiJisI/YJRFJFg8BI0mzWKxIP1SMTw4VY86kCCydE4c+XMZIxEPAyDldaOrAGzuycKGpAy8v42XSRLeLwU6SdCKvBhvSsxE7zB8vPZkIH6++YpdE5DQY7CQp5m4L3v88HwdP6LF0ThzmTR3OSzCIeojBTpJR09CGdduz0NphxmvPTEF0mJ/YJRE5JQY7ScLRM1XYtDsHY6KD8Orvx/JmI6I7wGAnUV0yW/DXT/NwOKsSv507Cr+4ZxiHXojuEIOdRFN5oRWvfXQS5m4rXtdOwfChvmKXRCQLDHYSxdcnK/D2nlzcPXIQnv7VaB7eRWRHDHbqVZ2XurF1Ty6+zanB7xZocN+EMA69ENkZg516TXlNC177KAsKhRvefHYqwgepxC6JSJYY7ORwgiDgy+8NePfTPEwbOxS/W6BBv778q0fkKPzXRQ7VZerGlr/l4Lu8WjyzaAzuTQgVuyQi2WOwk8NU17dhzQeZ6LYISHl2KsIHcuiFqDcw2MkhjufWYEPaGdwVEwztI2O46oWoFzHYya66LVZ8sF+Hvx8vwxMPj8TDkyO56oWolzHYyW4aWzrx2kdZqLvYgTVPT0bMMH+xSyJySQx2sovcknq8vv0Uhg1SYcNz03nMLpGIFPb8sLy8PCxfvhwHDx6058eShFmtAnZ/XYyX3zmB2RPD8V9PTWSoE4nslsFuMpmQlJSEhISEqwI7LS0NSUlJWLJkCQwGAwBAo9Fg2bJljquWJKWtw4RXt/2Avf8swUtPJmLJA7Fw5+XSRKK75VCMUqlEamrqVZdPNzc3IyMjA2lpadDpdEhJSUFqaqpDCyVpKalqxtoPT8J7QB+8tXI6Qvw9xS6JiH50y2BXKBQIDr76rsmcnBwkJibC3d0dGo0Ger3+tn7Yxo0bsWnTJpsKJWkQBAFf/WDAO3vzMGt8GJbPHwUPJS+XJpISmyZPjUYjVKr/32wiCAIAoKKiAtu2bYPJZEJYWBji4uKuep9Wq4VWq73qmVqttqUEEkGXqRtvZ1w+wEv7yBjMuIu7SImkyKZgV6lUKC4uvvJaobg8VB8WFoatW7fapzKSlJr6Nqz58CTM3RYe4EUkcTatiomPj0dmZiYsFgvy8/MRHh5u77pIQk7k1WDl+iMYEuSFN1dMY6gTSdxt9di1Wi10Oh08PT2RnZ2N5ORkzJ8/H4sXL4ZSqcTq1at79EPT09Ovmowlaeq2WPHRgQJ8fqwUv3loJOZO4S5SImfgJvw0QC4StVqNoqIiMUug62gydmHd9izUNrRj1eMJiIsIELskIvoXN8tO7jyla+SVNGDdjiyEhXhjw3PT4evNDUdEzoTBTlcIgoCMb0qw44sC/PLeKDw2O4YbjoickCjBzjF26WnrNGP9rtPIL2vEi09MwPi4gWKXREQ24hg7oay6BWs+zIRXfw8kL53AXaREToBj7HRDX/1gwDt7cnHv+DAsnzcKfTy4i5TI2THYXdQlswXv7MnFkTPV+OOieNybECZ2SURkJxxjd0G1De1Y++FJdJm6kfLsVAzjhiMiWeEYu4s5kVeLDWmnMToqCCuSxvIuUiInxTF2QrfFiu0HCrDvaCl+81Ac5k0dzl2kRDLFYHcBjS2dWLc9C+cbO7D6D5MwMpK7SInkjMEuc2eK6vDGzlMYPsQHqc/zLlIiV8DJU5myWAWkfVWE3V8XI+l+NRbNjOYuUiIXwclTGbrY2oWUnadgqG3FnxffhfjoILFLIiI74+SpCzlb2oDXd2RhUKAXNjw/Hf6qfmKXRES9jMEuE1argIxvzmHnwUIsnDECi2fHwN3dpntUiMjJMdhlwNhuwlu7TqPI0ISXnkxEQmyI2CURkYgY7E6u0NCE1z7KQoBPP6x/bjqC/XiAF5Gr46oYJyUIAj47VoYP9udjzqRILJ0TBw8lh16IiKtinFJ7pxkb0s8g51w9ViSNxUTNYLFLIqJexlUxMlJa1Yy1H52EZz8PrF85HYMCB4hdEhFJDIPdSQiCgAPf6fHeZ2cxa3wYlvHsdCK6AQa7E2hpu4SNn2Qjr7QBzz46FtPGDRW7JCKSMAa7xGUX1+GtXacR7OeJDc9Nx8AADr0Q0c0x2CXK3G3Fji8uH7O7aGY0ku6L5oYjIrotXO4oQdX1bXh9RxaM7SYes0tEPcbljhIiCAL+kVmBdz/Nw/jYEPxx0Rh49ecNR0R0LS53dAJtHSZs2p2DU4UX8PsFozFzfChvOCIimzDYJeBsaQNSPj4NX+++2PDcdAwO8hK7JCJyYgx2EXVbrEj7qgh/O3wOC6aPwGOzY3gsABHdMQa7SM43tuONHadQ39yJV56aiPgoXoZBRPbBYBfBN6cq8XZGLuKjAvHysruhGtBH7JKISEYY7L2ovdOMrXty8V1eLZbPG4XZd4dzgpSI7I7B3ksK9U14fecpDOinxPqV0xAa4i12SUQkU9yg5GAWq4DdXxcj7asiPDQ5EkvnxMJDycO7iMhxuEHJgeouduDNj0+jur4NK5PGYVxMsNglEZFMcIOSCI5lV2Pz7mzERgRg4/Mz4OvdV+ySiMhFMNjtzNhuwnufncWx7Go88dBIPDQ5ghOkRNSrGOx2YrUKOJxVgW37dQj264+UZ6ciYrCP2GURkQtisNuBodaILRk50Nca8fiDsXjgngi4K9hLJyJxMNjvQOelbuz6qgifHS3F5PghSH58PPxU/cQui4hcHIPdBoIg4PuztXj307Po66HAK8snIj6aRwIQkTQw2HvofGM73tmbh9xz9XhkVjQWzhjBdelEJCkM9ttk7rZi7z9LkH6oGKNHBGLzX+7l/aNEJEkM9tuQW1KPtzNy0XWpG88/Ng4TNYO4hJGIJIvBfhMXW7vw/mf5OJZdjYenROKx2THo35d/ZEQkbUyp67BYBRw8ocf2AzqEDVThrZXTuCadiJwGDwH7mZLKZmzJyMH5xg78du4ozBwfBgXXpBORE+EhYD9q7zRjxxcFOHBCj5kJoVg6Jw4+XjzfhYikiYeA3YQgCDh6phrvfXYWPl59sebpSYiLCBC7LCIim7l0sFfVtWLrnlwUGS7isdkxeHhKJJTuvEyaiJybSwb7JbMFuw8VI+ObEoyPC8GWv8xEkF9/scsiIrILlwv2rIILeGdvLgQBePGJCUiIDRG7JCIiu3KZYG9o7sRf9+UhM/88Fs6IwqKZUejXx2V+fSJyIbJPNovFis+/LcPHXxYiKtQPqc/P4EXSRCRrsg72gvImbMnIQXPbJTz9y3hMGzeURwEQkezJMtiN7SZ8sD8fX2dV4sGJw7DkwVh49fcQuywiol4hq2C3WgV8ffLy9XQh/v3xxp+mICrUT+yyiIh6lWyCXV9rxJa/5cBwntfTEZFrc/pg//n1dC8s5fV0ROTaJBHs+lrjdZ/f6hibqgtteP/zs+jbR4lXnpqI+CheT0dEJIlg177xjU3v6+PhjkdmRWHhdF5PR0T0E0kE+57XHrbpfQqFG8fRiYh+RhLB7qHkwVtERPbCRCUikhkGOxGRzNh1KKayshKbNm2Cm5sbXn75ZXh6etrz44mI6DbcssduMpmQlJSEhIQEHDx48MrztLQ0JCUlYcmSJTAYDACAvXv34oUXXsCCBQtw5MgRx1VNREQ3dMtgVyqVSE1NxdKlS688a25uRkZGBnbu3IlVq1YhJSUFAGA0GuHr64uQkBDU19c7rmoiIrqhWw7FKBQKBAcHX/UsJycHiYmJcHd3h0ajgV6vBwB4e3ujpaUFdXV1CAwMvOazNm7ciE2bNtmnciIiui6bxtiNRiNUKtWV1z/tEF24cCHWrl0LAHjppZeueZ9Wq4VWq73qmVqttqUEIiK6AZuCXaVSobi4+MprheLyiE5oaCjWrFljn8qIiMgmNi13jI+PR2ZmJiwWC/Lz8xEeHm7vuoiIyEa31WPXarXQ6XTw9PREdnY2kpOTMX/+fCxevBhKpRKrV6/u0Q9NT09Henq6TQUTEdHNuQm3OkLRwdRqNYqKisQsgYjI6dwsO7nzlIhIZkQ5BIxDMUREjsOhGCIiJ8ShGCIiF8JgJyKSGQY7EZHMcPKUiEhmOHlKROSEOHlKRORCGOxERDLDYCcikhlOnhIRyQwnT4mInBAnT4mIXAiDnYhIZhjsREQyw8lTIiKZ4eQpEZET4uQpEZELYbATEckMg52ISGYY7EREMsNgJyKSGQY7EZHMcB07EZHMcB07EZET4jp2IiIXwmAnIpIZBjsRkcww2ImIZIbBTkQkMwx2IiKZYbATEckMNygREckMNygRETkhblAiInIhDHYiIplhsBMRyYwok6c/p1arxS6BiEg+BAlZsGBBj74WHR3tyHLuyM1+FzE/t6fvv93vv53v62n73uprUm1/tn3Pvodtb//P5VAMEZHMSCrYH330UZu+JkWOqvdOP7en77/d77+d77O1fdn29vlctr3jSaXtRV/Hfie4Bt61sf1dF9v+5iTVY++pZ555RuwSSERsf9fFtr85p+6xExHRtZy6x05ERNdisBMRyQyDnYhIZhjsREQyI6tgr6ysxKpVq5CcnIyOjg6xy6FelJeXh+XLl+PgwYNil0K9bN++fXjxxRexYsUKmEwmscuRBMkHu8lkQlJSEhISEq76R5uWloakpCQsWbIEBoMBALB371688MILWLBgAY4cOSJWyWQnPWl7jUaDZcuWiVUq2VlP2n7evHlYvXo1NBoNqqqqxCpZUiQf7EqlEqmpqVi6dOmVZ83NzcjIyMDOnTuxatUqpKSkAACMRiN8fX0REhKC+vp6sUomO+lJ25O89LTt29raYDAYEBERIUa5kiP5YFcoFAgODr7qWU5ODhITE+Hu7g6NRgO9Xg8A8Pb2RktLC+rq6hAYGChCtWRPPWl7kpeetL3JZMK6deuwYsUKuLm5iVCt9Eg+2K/HaDRCpVJdef3THquFCxdi7dq12Lt3L6ZNmyZWeeRAN2r7iooKbNu2DZ988gl0Op1Y5ZED3ajtN2/ejOrqaqxfvx4VFRVilScpkjiPvadUKhWKi4uvvFYoLv//FBoaijVr1ohVFvWCG7V9WFgYtm7dKlZZ1Atu1PYrV64UqyTJcsoee3x8PDIzM2GxWJCfn4/w8HCxS6JewrZ3XWz72+cUPXatVgudTgdPT09kZ2cjOTkZ8+fPx+LFi6FUKrF69WqxSyQHYdu7Lra97XgIGBGRzDjlUAwREd0Yg52ISGYY7EREMsNgJyKSGQY7EZHMMNiJiGSGwU5EJDMMdiIimfk/zKuT6t7NLpUAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\"\"\"Plot MSD\"\"\"\n", + "plt.loglog(msd_avg)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 439, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.34813812, 1.2741015, 2.647874, 4.3709273, 6.3643317, 8.562542, 10.910605, 13.362199, 15.878576, 18.427135, 20.981758, 23.52136, 26.029282, 28.492788, 30.90219, 33.2509, 35.534565, 37.75038, 39.89707, 41.974033, 43.98095, 45.918133, 47.78734, 49.5914, 51.333023, 53.015324, 54.641724, 56.21574, 57.74034, 59.21842, 60.652515, 62.04567, 63.400356, 64.71856, 66.00155, 67.250824, 68.46969, 69.660065, 70.823746, 71.96182, 73.07512, 74.16329, 75.22746, 76.267624, 77.284134, 78.27746, 79.24828, 80.197464, 81.12534, 82.03311, 82.92249, 83.79491, 84.651184, 85.491936, 86.317314, 87.12684, 87.92105, 88.701515, 89.46894, 90.22471, 90.96953, 91.70468, 92.43193, 93.152596, 93.86661, 94.57387, 95.27456, 95.96885, 96.65563, 97.33437, 98.00493, 98.66652, 99.319984, 99.967476, 100.609375, 101.245804, 101.87765, 102.50451, 103.125534, 103.74044, 104.34754, 104.94508, 105.53214, 106.10796, 106.672676, 107.225426, 107.766594, 108.296616, 108.815346, 109.32294, 109.81912, 110.30355, 110.777275, 111.24183, 111.698845, 112.14981, 112.595505, 113.03822, 113.48041, 113.92206, 114.364006, 114.80576, 115.24795, 115.69098, 116.13387, 116.57511, 117.013855, 117.45001, 117.88383, 118.315056, 118.74291, 119.16574, 119.5837, 119.99798, 120.40952, 120.818596, 121.22681, 121.63534, 122.04548, 122.45892, 122.87619, 123.29721, 123.72154, 124.146904, 124.57166, 124.99492, 125.41684, 125.836136, 126.25196, 126.66288, 127.068794, 127.469345, 127.86187, 128.24596, 128.62286, 128.99413, 129.36261, 129.72816, 130.09143, 130.45364, 130.81555, 131.17735, 131.53787, 131.89755, 132.25598, 132.61476, 132.97379, 133.33316, 133.69247, 134.05156, 134.41103, 134.77235, 135.13736, 135.50645, 135.8795, 136.25783, 136.64098, 137.02933, 137.42108, 137.81506, 138.21031, 138.6064, 139.00264, 139.39864, 139.79347, 140.18669, 140.5778, 140.96692, 141.3536, 141.73607, 142.11363, 142.4859, 142.85222, 143.21268, 143.56712, 143.91518, 144.2572, 144.59436, 144.92719, 145.2555, 145.58058, 145.9044, 146.22737, 146.54959, 146.87091, 147.19122, 147.51085, 147.82996, 148.14868, 148.46783, 148.78761, 149.10857, 149.43124, 149.75533, 150.08118, 150.4091, 150.7392, 151.07094, 151.40532, 151.74324, 152.0845, 152.42871, 152.77493, 153.12257, 153.4717, 153.8214, 154.16893, 154.51227, 154.85153, 155.18715, 155.51936, 155.84863, 156.17549, 156.50032, 156.82468, 157.14928, 157.47502, 157.79926, 158.12263, 158.44583, 158.76866, 159.09123, 159.41125, 159.72737, 160.03738, 160.34085, 160.63832, 160.92969, 161.21524, 161.49547, 161.77098, 162.0406, 162.30403, 162.56102, 162.81308, 163.06026, 163.30179, 163.53682, 163.76416, 163.98448, 164.19734, 164.40302, 164.60352, 164.79962, 164.99248, 165.18314, 165.37401, 165.5658, 165.75833]\n" + ] + } + ], + "source": [ + "print(list(msd_avg)) # in post-processing, simply average over all replicates" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pappuPython", + "language": "python", + "name": "pappupython" + }, + "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.9.16" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}