omm_cookbook_ion_distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "698caaa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk.openmm import app\n",
    "from simtk import openmm\n",
    "from simtk import unit\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c50e2e14",
   "metadata": {},
   "source": [
    "# Prepare a template tip3p water \n",
    "This will be used by `Modeller` to build a water box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "45d7d265",
   "metadata": {},
   "outputs": [],
   "source": [
    "wat = \"\"\"REMARK   1 CREATED WITH OPENMM 7.1.1, 2018-04-06\n",
    "CRYST1   24.668   24.668   24.668  90.00  90.00  90.00 P 1           1 \n",
    "HETATM    1  O   HOH A   1      20.854   7.945   2.695  1.00  0.00           O  \n",
    "HETATM    2  H1  HOH A   1      20.168   7.772   3.339  1.00  0.00           H  \n",
    "HETATM    3  H2  HOH A   1      21.142   8.838   2.885  1.00  0.00           H  \n",
    "END\"\"\"\n",
    "\n",
    "f = open('wat.pdb', 'w')\n",
    "f.write(wat)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e4cc13f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1428 1429\n"
     ]
    }
   ],
   "source": [
    "#load pdb file\n",
    "pdb = app.PDBFile('wat.pdb')\n",
    "# we don't need the protein forcefield here, just water (includes ions):\n",
    "forcefield = app.ForceField('amber14/tip3pfb.xml') \n",
    "\n",
    "#instantiate a Modeller object using the topology and xyz coordinates,\n",
    "#and add salt water:\n",
    "modeller = app.Modeller(pdb.topology, pdb.positions)\n",
    "modeller.addSolvent(forcefield, ionicStrength = 0.15*unit.molar)\n",
    "\n",
    "#take the indices of the first sodium and first chlorine atom\n",
    "sod = [i.index for i in modeller.topology.atoms() if i.residue.name=='NA'][0]\n",
    "cla = [i.index for i in modeller.topology.atoms() if i.residue.name=='CL'][0]\n",
    "\n",
    "print(sod, cla)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6afa2973",
   "metadata": {},
   "source": [
    "# Set up the simulation system\n",
    "With a topology and forcefield, we can create a `System` object. Some other parameters inherent to the calculation of potential energy are also set here (Particle Mesh Ewald (PME) to estimate long-range electrostatics, a cut-off for short-range electrostatics, and bonds to hydrogens being fixed). \n",
    "\n",
    "Following this we can modify the potential energy function by adding external biases (in this case, `CustomBondForce`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "56b2b7bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME,\n",
    "         nonbondedCutoff=0.9*unit.nanometer, constraints=app.HBonds)\n",
    "\n",
    "\n",
    "#Generate a custom bond force biasing the distance between the two ions.\n",
    "#The atoms are not covalently bound, but this is a pairwise restraint \n",
    "#like a covalent bond, so the same terminology is used. In this case, \n",
    "#the actual functional form is the same as a covalent bond, but any\n",
    "#other energy expression could be used.\n",
    "\n",
    "#this string defines a harmonic distance restraint, with distance during the \n",
    "#simulation referred to by openmm as 'r', centred on input parameter 'r0' (nanometers),\n",
    "#and scaled by strength 'k' (kJ/mol). 'r0' and 'k' are variables that can take \n",
    "#any user-defined name and value.\n",
    "nrg_expr = \"k*(r-r0)^2\"\n",
    "ion_force = openmm.CustomBondForce(nrg_expr)\n",
    "#set the user-defined parameters:\n",
    "ion_force.addGlobalParameter('k', 100*unit.kilocalorie_per_mole**2)\n",
    "ion_force.addGlobalParameter('r0', 1*unit.nanometer) \n",
    "\n",
    "#now add a \"bond\"\n",
    "ion_force.addBond(sod, cla, [])\n",
    "\n",
    "ion_force.setUsesPeriodicBoundaryConditions(True) \n",
    "system.addForce(ion_force)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2cb43c3",
   "metadata": {},
   "source": [
    "# Generate a `simulation`\n",
    "The `System` can calculate potential energy, so now we can simulate through time by adding temperature and an integrator, bundled into a `Simulation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a28b8dee",
   "metadata": {},
   "outputs": [],
   "source": [
    "#integrator attributes:\n",
    "TEMPERATURE = 298*unit.kelvin\n",
    "FRICTION = 1/unit.picosecond\n",
    "TIMESTEP = 2*unit.femtosecond\n",
    "integrator = openmm.LangevinIntegrator(TEMPERATURE, FRICTION, TIMESTEP)\n",
    "\n",
    "platform = openmm.Platform.getPlatformByName('CPU')\n",
    "\n",
    "simulation = app.Simulation(modeller.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(modeller.positions)\n",
    "simulation.minimizeEnergy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4faaef02",
   "metadata": {},
   "source": [
    "# Simulate and observe ion distance\n",
    "Note that the distance is about 1 nm, the value to which `r0` was set.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1cf06bf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#all this cell is just to calculate the distance between \n",
    "#the ions. It's not typically necessary in production.\n",
    "\n",
    "startingState = simulation.context.getState()\n",
    "#a list of the x, y, and z lengths of the periodic box\n",
    "BOXLENGTHS = [i[count]/unit.nanometer for count, i in enumerate(startingState.getPeriodicBoxVectors())]\n",
    "\n",
    "def minImgDist(pos1, pos2):\n",
    "    \"\"\"\n",
    "    given known periodic box size, this function takes the minimum distance\n",
    "    between two positions i.e. it performs the minimum image correction\n",
    "    Note: mdtraj is a much better way to do this.\n",
    "    \"\"\"\n",
    "    \n",
    "    #subtract x, y and z coordinates to get displacement\n",
    "    #take the absolute value to simplify the minimum image correction\n",
    "    displacement_xyz = [abs(i) for i in pos1 - pos2]    \n",
    "    \n",
    "    #if the displacement is greater than half a boxlength, then subtract\n",
    "    #one whole boxlength to get the smallest displacement between imaged copies\n",
    "    displacement_xyz = [displ if displ<(blength/2) else displ-blength for displ,blength in zip(displacement_xyz, BOXLENGTHS) ]\n",
    "    distance = math.sqrt(sum([i**2 for i in displacement_xyz]))\n",
    "    \n",
    "    return distance\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b269bdb2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated 0 steps\n",
      "Simulated 100 steps\n",
      "Simulated 200 steps\n",
      "Simulated 300 steps\n",
      "Simulated 400 steps\n",
      "Simulated 500 steps\n",
      "Simulated 600 steps\n",
      "Simulated 700 steps\n",
      "Simulated 800 steps\n",
      "Simulated 900 steps\n"
     ]
    }
   ],
   "source": [
    "nsteps = 100\n",
    "\n",
    "distances = list()\n",
    "for i in range(10):\n",
    "    print(f'Simulated {simulation.currentStep} steps')\n",
    "    #advance atom positions by 'nsteps' timesteps\n",
    "    simulation.step(nsteps)\n",
    "    \n",
    "    #fetch the current state in order to access positions\n",
    "    curr_state = simulation.context.getState(getPositions=True)\n",
    "    positions = curr_state.getPositions()\n",
    "    \n",
    "    #calculate ion distance:\n",
    "    sod_pos = positions[sod]/unit.nanometer\n",
    "    cla_pos = positions[cla]/unit.nanometer\n",
    "    \n",
    "    distance = minImgDist(sod_pos, cla_pos)\n",
    "    distances.append(distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f1192c8b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Timesteps, (200 fs)')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAebUlEQVR4nO3de5RcZZnv8e+vqztJ5w50o5AEExxAo4JAy0UEUVQueozXkch4d3JYitc5Ksw6I8y4Znmdc3AUzEQmg6AHjiMICEicUYRRBNKRSEIQzQGBAJoOmECSTvr2nD/2rqS6u6qrutO7qjv791mrVtV+97t3Pb2T2s9+9+V9FRGYmVl+NTU6ADMzaywnAjOznHMiMDPLOScCM7OccyIwM8u55kYHMFptbW2xcOHCRodhZjaprFmzZktEtJebN+kSwcKFC+ns7Gx0GGZmk4qkRyvN86khM7OccyIwM8s5JwIzs5xzIjAzyzknAjOznHMiMDPLOScCM7OccyIwM8s5JwIzs5xzIjAzyzknAjOznHMiMDPLOScCM7OccyIwM8u5zBKBpJWSNktaP0Kd0yWtlfSApDuyisXMzCrLskVwJXBWpZmS5gKXA2+OiJcA78wwFjMzqyCzRBARdwLPjFDl3cD1EfFYWn9zVrGYmVlljbxGcCRwgKSfS1oj6b2VKkpaJqlTUmdXV1cdQzQz2/81MhE0A8cDbwTOBP5O0pHlKkbEiojoiIiO9vayQ26amdkYNXLM4k3AlojYAeyQdCdwDPC7BsZkZpY7jWwR3AicKqlZ0nTgRODBBsZjZpZLmbUIJF0DnA60SdoEXAy0AETE8oh4UNJtwP3AAHBFRFS81dTMzLKRWSKIiKU11Pkq8NWsYjAzs+r8ZLGZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlXGaJQNJKSZsljTjqmKRXSOqX9I6sYjEzs8qybBFcCZw1UgVJBeDLwKoM4zAzsxFklggi4k7gmSrVPgZcB2zOKg4zMxtZw64RSJoHvBVYXkPdZZI6JXV2dXVlH5yZWY408mLxpcDnIqK/WsWIWBERHRHR0d7enn1kZmY50tzA7+4ArpUE0AacI6kvIm5oYExmZrnTsEQQEYuKnyVdCdzsJGBmVn+ZJQJJ1wCnA22SNgEXAy0AEVH1uoCZmdVHZokgIpaOou77s4rDzMxG5ieLzcxyzonAzCznnAjMzHLOicDMLOecCMzMcs6JwMws55wIzMxyzonAzCznnAjMzHLOicDMLOecCMzMcs6JwMws55wIzMxyzonAzCznnAjMzHLOicDMLOcySwSSVkraLGl9hfnnSbo/fd0l6ZisYjEzs8qybBFcCZw1wvxHgFdHxNHAF4AVGcZiZmYVVB2qUtLJwF8BpwKHAN3AeuAW4LsRsa3cchFxp6SFldYbEXeVTN4NzK89bDMzGy8jtggk/Rj4MLCK5Oj+EGAx8D+BacCNkt48DnF8CPjxCHEsk9QpqbOrq2scvs7MzIoUEZVnSm0RsWXEFYxQJ20R3BwRLx1h+dcAlwOvioinqwXc0dERnZ2d1aqZmVkJSWsioqPcvBFPDQ3dwUuaXbpMRDxTLVFUCexo4Arg7FqSgJmZjb+q1wgAJP134B9Irg8UmxABHD7WL5Z0GHA98J6I+N1Y12NmZvumpkQA/A/gJaM5+pd0DXA60CZpE3Ax0AIQEcuBzwMHAZdLAuir1GwxM7Ps1JoI/h+wczQrjoilVeZ/mORCtJmZNVCtieAi4C5J9wC7i4UR8fFMojIzs7qpNRH8C/AzYB0wkF04ZmZWb7Umgr6I+HSmkZiZWUPU2sXE7elDXYdIOrD4yjQyMzOri1pbBO9O3y8qKdun20fNzGxiqCkRRMSirAMxM7PGqLVFgKRXAgsZ/GTxVRnEZGZmdVTrk8VXAy8E1gL9aXEATgRmZpNcrS2CDmBxjNRDnZmZTUq13jW0Hnh+loGYmVlj1NoiaAM2SLqXwU8Wj8dYBGZm1kC1JoJLsgzCzMwaZ8REIEmRuKNanfEPzaw+brjvCb666iGe3NrNoXNb+cyZR/GWY+flNo6JwNuivqq1CG6XdB1wY0Q8ViyUNAV4FfA+4HaSgerNJp0b7nuCi65fR3dvcjPcE1u7uej6dQB13fFMlDgmAm+L+qs2VOU04IPAecAiYCvJWMUF4CfAZRGxNvMoS3ioyv1H1kd9EcHOnn62dveybWcvW7t7eLa7l607e5Oy7l6+c9cf2NnTP2zZqc1NdCw8ACGS4TIS6dgZCPaUq2Te3s9750gMKldpeTrj9t9uprt3eH+Os6Y18/HXHsG0KQWmtxSYPqVA65QC06c009pS/Jy8prUUmNrctCfGsRiPf5OBgWBnbz87dvelr362Fz/39LF9dx87h5Ul9bfv7uO+x/5Mb//w/dKc1mYufdexLGqbwfwDWmku1Hqvi8HIQ1WOmAiGrKSF5KJxd0RsraH+SuBNwOZyYxYr+d/6deAckrEO3h8Rv6623smcCCZKc3cixDH0qA+gtaXAF9/2smGx9PYPsC3dgW/r7in5vPe9+Nq6s4et3b17dvh9A5X/f7cUVHaHU3T8Cw6g+PsIoPhTCfZO7BmuLyDSqYjBdUt/Y3vLY1CdjZu3V4xjNJpEkiSmFGgdlDgKtLYk5dNbSsv2JpYNT23j+6sfp6dkm7QUxFuOnceRB89iR09furPeu5PfvruPnT17d+I7dvexs7efWk8Wt7YUmDG1mZlTk/cZU5u595Fnqi7XUhALDpzO4W0zObx9Bovaktfh7TNonzl1n5Lh/mpcEsEYvvQ0YDtwVYVEcA7wMZJEcCLw9Yg4sdp6J2siGM2Ob6LGERH0DQS9/QP09gd9/QN7pvv6g76BYnnQO5CUJXWHzB8Y4O9/tIGtO3uHfce0liZevmAu27r72LYz2envKHPEXmrWtGbmtLYwd3oLc1unMKe1hTnTW5Ky1pY98+ak8+am86ZPKfCqL9/OE1u7h61z3txWfnnha2vcqvvulC/9rGwch86dxqpPnkZ3Tz/dvf3s7Elee6f76C6W9faXfO4bUi957UqXKc4bKVGWU7rjnj6lmZlTm5kxtViW7MhnTNm7U69clixfaBq+w660LQ6ZM41vLD2Wh7fs4JEtO3ikK31/egc9fXtbUzOnNg9KDIvaZnB420wWtk1n1rSWUf29+5OGJIL0ixcCN1dIBP8C/DwirkmnHwJOj4inRlrnWBJBVkfAAwOxp6m7fVcfz+3u47ldyeftu3t5blc6vbuPa+59rOwpiCmFZMc3JmM46Fn72FZ6+oefgmhuEocdND3ZYfcP0DuQ7uhLduqj3WmM1SsWHjBsp118Tz5P2bOTnzWteZ9OEewPCXpf9PYP7EkKJ3/xp5T7Fxbwm0vewIwKO+7xNtpt0T8QPLm1O0kKW3bwcNf2Pcniia3dg1onB8+aOihBLEpbFAsOmM6U5uH/jyZC63m8jJQIau5rKAPzgMdLpjelZSMmgtEqf+Hpfnb19vHqow6uugPfnr4ndXr3TG/f1cf2nr6amsAzphTKJgGAnv6BMf24ouxPtrpySQCgbyB48SGzaWkSzYUmWgqiuamJ5oJoKTTRXCwvmd9SSOen9Urnl5YX19VSXG+hiXNX/Io/Pbt7WBzz5rby7+e/ckx/21gUf9SN/rE3Ko6WQhNzWpuY09rCoXNbK7RKWpldxyPp0W6LQlNymmjBgdM57cj2QfN29fbz2DM79yaHtBWx6oE/8cyOnsHrOKA1TRIzWdQ2gye3drPyF4+wK21t7M8XrUdzjeAFwBER8Z+SWoHmiHiuyjILqdwiuAX4YkT8Ip3+KfDZiFhTpu4yYBnAYYcddvyjjz5aU8xQuZlZixlTCsycljRlZ05rYXbx89RmZk5rZtbUZmZNaymp05zWKSmbmhxFVYpjopyCqHccE+VI3PbK27/J1p09Ja2I9H3LDh7Zsp1dZS7cF7W2FHjbcfP2tFBnt+5trc6etvfzzGnj14Iaj5bJPrcIJP01yY74QJLO5+YDy4EzRhXJYJuABSXT84Eny1WMiBXACkhODY3mS54cIQl88W0vq7xTnzq+zeDPnHlU2R/ZZ848aty+YzLFMVGOxG2vvP2bzJ0+hWMPm8Kxhx0wqHxgIPjTc7s4+Ys/K7tcd28/P17/R7Z199I/wulSKbleUZogZrc2D04a0/d+nt06uM7U5gJQn9tpaz019FHgBOAegIj4vaSD9/G7bwIukHQtycXibdWuD4xFpebuvLmtLD3hsPH+uoomyo9sosRRjGV/3clMVv43gaYmccicVuaNsO/45YWv3XN7cvGOtWdL7l57dlffnrLS8ke27NjzeaRWByQ3Tsye1sIzO3qGXZ/r7u3nq6seqnsi2B0RPXvuoZaaYeST1JKuAU4H2iRtAi4GWgAiYjlwK8kdQxtJbh/9wBjir2qiHAHDxPmRTZQ4zCayavsOSXvuhDp0buuo17+7r59nu9OEsau3bOJ4truP/9v5eNnlRzrbMVq1JoI7JP0t0Crp9cBHgB+NtEBELK0yP0haGpmaSEfAZjZ5ZL3vmNpcoH1WgfZZU0es94uNWypexB8vNV0sltQEfAh4A8ndZKuAKxrRx9BkfY7AzGwsxusi/njcPtoKrIyIb6crLKRlO2uOwszMRq0eZzVqTQQ/BV5H8qQwJEngJ0D9bvg2M8uprK/r1fpI5rSI2NMZSvp5ejYhmZlZPdWaCHZIOq44Iel4YPwuWZuZWcPUemrok8C/Syo+8HUI8K5MIjIzs7qqKRFExGpJLwKOIrlr6LcRMbzrSDMzm3RG0+ncK4CF6TLHSiIirsokKjMzq5ta+xq6mqSPobVA8WbWAJwIzMwmuVpbBB3AYg9Sb2a2/6n1rqH1wPOzDMTMzBqj1hZBG7BB0r3AntFEIuLNmURlZmZ1U2siuCTLIMzMrHFqvX30jqwDMTOzxqjpGoGkkyStlrRdUo+kfknPZh2cmZllr9aLxd8ElgK/J+lw7sNpmZmZTXK1JgIiYiNQiIj+iPg3ktHHRiTpLEkPSdoo6cIy8+dI+pGk30h6QFImo5SZmVlltV4s3ilpCrBW0leAp4AZIy2QjllwGfB6koHqV0u6KSI2lFT7KLAhIv6bpHbgIUnfi4ieUf8lZmY2JrW2CN6T1r0A2AEsAN5WZZkTgI0R8XC6Y78WWDKkTgCzlAyGPBN4BuirMSYzMxsHtSaCt0TEroh4NiL+PiI+DbypyjLzgNJRlzelZaW+CbwYeBJYB3wiIgaGrkjSMkmdkjq7urpqDNnMzGpRayJ4X5my91dZRmXKhnZRcSZJ/0WHAi8Hvilp9rCFIlZEREdEdLS3t1cN1szMajfiNQJJS4F3A4sk3VQyazbwdJV1byI5hVQ0n+TIv9QHgC+lfRhtlPQI8CLg3hpiNzOzcVDtYvFdJBeG24B/Kil/Dri/yrKrgSMkLQKeAM4lSSqlHgPOAP5L0vNIxjt4uLbQzcxsPIyYCCLiUeBRSa8DuiNiQNKRJEft66os2yfpAmAVUABWRsQDks5P5y8HvgBcKWkdyamkz0XEln3+q8zMrGaqpWdpSWuAU4EDgLuBTmBnRJyXbXjDdXR0RGdnZ72/1sxsUpO0JiI6ys2r9WKxImInyS2j34iItwKLxytAMzNrnJoTgaSTgfOAW9Ky0QxzaWZmE1StieCTwEXAD9Pz/IcDt2cWlZmZ1c1ouqG+o2T6YeDjWQVlZmb1U+05gksj4pOSfsTwh8E8QpmZ2X6gWovg6vT9a1kHYmZmjVHtOYI16btHKDMz209VOzW0jjKnhIoi4uhxj8jMzOqq2qmhYg+jH03fi6eKzgN2ZhKRmZnVVS1dTCDplIg4pWTWhZJ+CfxDlsGZmVn2an2OYIakVxUnJL2SKiOUmZnZ5FDr08EfAlZKmkNyzWAb8MHMojIzs7qp9YGyNcAx6aAxioht2YZlZmb1Mqr+giLi2awCMTOzxqj1GoGZme2nnAjMzHKu5lND6Z1CC0uXiYirqixzFvB1khHKroiIL5WpczpwKdACbImIV9cak5mZ7buaEoGkq4EXAmuB/rQ4gIqJQFIBuAx4PclA9qsl3RQRG0rqzAUuB86KiMckHTyGv8HMzPZBrS2CDmBx1DKu5V4nABvTLquRdC2wBNhQUufdwPUR8RhARGwexfrNzGwc1HqNYD3w/FGuex7weMn0prSs1JHAAZJ+LmmNpPeWW5GkZZI6JXV2dXWNMgwzMxtJrS2CNmCDpHuB3cXCKuMRqEzZ0BZFM3A8cAbQCvxK0t0R8btBC0WsAFZAMnh9jTGbmVkNak0El4xh3ZuABSXT84Eny9TZEhE7gB2S7gSOAX6HmZnVRU2nhtLxCH4LzEpfD9YwRsFq4AhJiyRNAc4FbhpS50bgVEnNkqYDJwIPjuYPMDOzfVNTIpD0l8C9wDuBvwTukfSOkZaJiD7gAmAVyc79++nA9+dLOj+t8yBwG3B/uv4rImL9WP8YMzMbPdVyI5Ck3wCvL97VI6kd+M+IOCbj+Ibp6OiIzs7Oen+tmdmkJmlNRHSUm1frXUNNQ27tfHoUy5qZ2QRW68Xi2yStAq5Jp98F3JpNSGZmVk+1dkP9GUlvB04huS10RUT8MNPIzMysLmruaygirgOuyzAWMzNrgBETgaTnGP4QGCStgoiI2ZlEZWZmdVNt8PpZ9QrEzMwaw3f+mJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc5lmggknSXpIUkbJV04Qr1XSOqvNuqZmZmNv8wSgaQCcBlwNrAYWCppcYV6XyYZ0tLMzOosyxbBCcDGiHg4InqAa4ElZep9jKR7681l5pmZWcayTATzgMdLpjelZXtImge8FVg+0ookLZPUKamzq6tr3AM1M8uzLBOBypQNHdvgUuBzEdE/0ooiYkVEdERER3t7+3jFZ2ZmjGKEsjHYBCwomZ4PPDmkTgdwrSSANuAcSX0RcUOGcZmZWYksE8Fq4AhJi4AngHOBd5dWiIhFxc+SrgRudhIwM6uvzBJBRPRJuoDkbqACsDIiHpB0fjp/xOsCZmZWH1m2CIiIW4Fbh5SVTQAR8f4sYzEzs/L8ZLGZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlXKaJQNJZkh6StFHShWXmnyfp/vR1l6RjsozHzMyGyywRSCoAlwFnA4uBpZIWD6n2CPDqiDga+AKwIqt4zMysvCxbBCcAGyPi4YjoAa4FlpRWiIi7IuLP6eTdwPwM4zEzszKyTATzgMdLpjelZZV8CPhxuRmSlknqlNTZ1dU1jiGamVmWiUBlyqJsRek1JIngc+XmR8SKiOiIiI729vZxDNHMzJozXPcmYEHJ9HzgyaGVJB0NXAGcHRFPZxiPmZmVkWWLYDVwhKRFkqYA5wI3lVaQdBhwPfCeiPhdhrGYmVkFmbUIIqJP0gXAKqAArIyIBySdn85fDnweOAi4XBJAX0R0ZBWTmZkNp4iyp+0nrI6Ojujs7Gx0GGZmk4qkNZUOtP1ksZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeWcE4GZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORGYmeVcpolA0lmSHpK0UdKFZeZL0j+n8++XdFyW8ZiZ2XCZJQJJBeAy4GxgMbBU0uIh1c4Gjkhfy4BvZRWPmZmVl2WL4ARgY0Q8HBE9wLXAkiF1lgBXReJuYK6kQzKMyczMhshs8HpgHvB4yfQm4MQa6swDniqtJGkZSYsBYLukh8Y31LprA7Y0OogJxNtjMG+PvbwtBtuX7fGCSjOyTAQqUxZjqENErABWjEdQE4GkzkqDSOeRt8dg3h57eVsMltX2yPLU0CZgQcn0fODJMdQxM7MMZZkIVgNHSFokaQpwLnDTkDo3Ae9N7x46CdgWEU8NXZGZmWUns1NDEdEn6QJgFVAAVkbEA5LOT+cvB24FzgE2AjuBD2QVzwSz35zmGifeHoN5e+zlbTFYJttDEcNOyZuZWY74yWIzs5xzIjAzyzkngjqStEDS7ZIelPSApE80OqZGk1SQdJ+kmxsdS6NJmivpB5J+m/4fObnRMTWSpE+lv5P1kq6RNK3RMdWTpJWSNktaX1J2oKT/kPT79P2A8fguJ4L66gP+JiJeDJwEfLRMtxt58wngwUYHMUF8HbgtIl4EHEOOt4ukecDHgY6IeCnJDSfnNjaqursSOGtI2YXATyPiCOCn6fQ+cyKoo4h4KiJ+nX5+juSHPq+xUTWOpPnAG4ErGh1Lo0maDZwG/CtARPRExNaGBtV4zUCrpGZgOjl7xigi7gSeGVK8BPhO+vk7wFvG47ucCBpE0kLgWOCeBofSSJcCnwUGGhzHRHA40AX8W3qq7ApJMxodVKNExBPA14DHSLqc2RYRP2lsVBPC84rPWqXvB4/HSp0IGkDSTOA64JMR8Wyj42kESW8CNkfEmkbHMkE0A8cB34qIY4EdjFOzfzJKz30vARYBhwIzJP1VY6PafzkR1JmkFpIk8L2IuL7R8TTQKcCbJf2BpGfa10r6bmNDaqhNwKaIKLYQf0CSGPLqdcAjEdEVEb3A9cArGxzTRPCnYg/N6fvm8VipE0EdSRLJOeAHI+J/NTqeRoqIiyJifkQsJLkI+LOIyO0RX0T8EXhc0lFp0RnAhgaG1GiPASdJmp7+bs4gxxfPS9wEvC/9/D7gxvFYaZa9j9pwpwDvAdZJWpuW/W1E3Nq4kGwC+RjwvbRvrofJT5crw0TEPZJ+APya5G67+8hZdxOSrgFOB9okbQIuBr4EfF/Sh0iS5TvH5bvcxYSZWb751JCZWc45EZiZ5ZwTgZlZzjkRmJnlnBOBmVnOORHYhCHpIElr09cfJT2Rft4u6fIMv/d0SXV9WEnSpZJOSz9/T9JDaS+bK9OHDkmHcP1nSRsl3S/puJLlz0qX2Sip7BPIktol3ZN2WXFqhTpfk/TaLP5GmzycCGzCiIinI+LlEfFyYDnwv9PpmRHxkQy/+nTq+NSqpAOBk9JOxQC+B7wIeBnQCnw4LT8bOCJ9LQO+lS5fAC5L5y8GllboxfYM4LcRcWxE/FeFcL5BjruysIQTgU146RH7zennSyR9R9JPJP1B0tskfUXSOkm3lRxNHy/pDklrJK0qeSz/45I2pEfY16ad/50PfCptfZyaHklfJ2l1+jql5LuvlvSztD/4v07LD5F0Z7r8+kpH3yXeAdxWnIiIWyMF3AvMT2ctAa5KZ90NzE3/jhOAjRHxcET0kHTRsWTINns58BXgnDSuGZKuTONbJ+lT6Xc/Chwk6flj+bex/YOfLLbJ6IXAa0iOhn8FvD0iPivph8AbJd1CcqS7JCK6JL0L+EfggyRHv4siYrekuRGxVdJyYHtEfA1A0v8haY38QtJhwCrgxel3H00ylsQM4L70u5YCqyLiH9Oj9elV4j+FpC+hQdIk9h6SMRog6aL88ZIqm9KycuUnlq4rItZK+jxJf/4XSDoemJf27Y+kuSXVf53GdF2VuG0/5URgk9GPI6JX0jqSAUuKR9frgIXAUcBLgf9IuqmhQNKVMcD9JN043ADcUGH9rwMWp8sCzJY0K/18Y0R0A92Sbic5Ol8NFM/t3xARa6vEfwhJl9NDXQ7cWXIaR2XqxAjlI3kYOFzSN4BbgNIunTeT9PBpOeVTQzYZ7QaIiAGgN/b2kzJAcnAj4IHi9YaIeFlEvCGt80aS8+vHA2uUDHoyVBNwcsny89KBhGD4DjfSc/2nAU8AV0t6b5X4u4FBwy5KuhhoBz5dUrwJWFAyPZ9kcJZK5RVFxJ9JRj37OfBRBg8GNC2NyXLKicD2Rw8B7UrH/JXUIuklkpqABRFxO8mAOHOBmcBzwKyS5X8CXFCcSM+3Fy2RNE3SQSQXmVdLegHJ2ArfJuld9rh0uasknVAmvgeBvyhZ/4eBM4GlaXIrugl4b3r30Ekkg7M8RdICOULSIiUd1J2b1q1IUhvQFBHXAX/H4C6ujwTWl13QcsGJwPY76QXUdwBflvQbYC3JXUEF4LvpKaX7SK4DbAV+BLy1eLGYdKzc9ILyBpKLyUX3kpxauRv4QkQ8SZIQ1kq6D3g7ydjDkFxPeIrhbkmXKVoOPA/4VRrD59PyW0lO6WwEvg18JP37+kgS1SqSpPL9iHigymaZB/xcSa+3VwIXwZ7rEn8BdFZZ3vZj7n3UrEaSLqHkonKVurOBf42Ist0ES/oF8KZGj0ss6a3AcRHxd42MwxrLLQKzDETEs5WSQOpvgMPqFc8ImoF/anQQ1lhuEZiZ5ZxbBGZmOedEYGaWc04EZmY550RgZpZzTgRmZjn3/wFsacY/lAGAnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(1, 11), distances, '-o')\n",
    "plt.ylim(0,1.75)\n",
    "plt.ylabel('Ion distance (nm)')\n",
    "plt.xlabel(f'Timesteps, ({TIMESTEP*nsteps})')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca2d3f92",
   "metadata": {},
   "source": [
    "# Change the bond parameters mid-simulation\n",
    "The user-defined parameters were `r0`, the centre of the restraining force, and `k`, the strength. To push/pull atoms apart, `r0` can be changed mid-simulation. This might be useful in umbrella sampling to measure the PMF of ion separation.\n",
    "\n",
    "`r0` was added as a `Global` force, accessible outside the force object, so it can be changed quickly on the fly without reinitializing the whole simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "df5134b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#change the value of 'r0'\n",
    "simulation.context.setParameter('r0', 0.3*unit.nanometer)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e750fb2c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated 1000 steps\n",
      "Simulated 1100 steps\n",
      "Simulated 1200 steps\n",
      "Simulated 1300 steps\n",
      "Simulated 1400 steps\n",
      "Simulated 1500 steps\n",
      "Simulated 1600 steps\n",
      "Simulated 1700 steps\n",
      "Simulated 1800 steps\n",
      "Simulated 1900 steps\n"
     ]
    }
   ],
   "source": [
    "for i in range(10):\n",
    "    print(f'Simulated {simulation.currentStep} steps')\n",
    "    #advance atom positions by 'nsteps' timesteps\n",
    "    simulation.step(nsteps)\n",
    "    \n",
    "    #fetch the current state in order to access positions\n",
    "    curr_state = simulation.context.getState(getPositions=True)\n",
    "    positions = curr_state.getPositions()\n",
    "    \n",
    "    #calculate ion distance:\n",
    "    sod_pos = positions[sod]/unit.nanometer\n",
    "    cla_pos = positions[cla]/unit.nanometer\n",
    "    \n",
    "    distance = minImgDist(sod_pos, cla_pos)\n",
    "    distances.append(distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b6ddd5a3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Timesteps, (200 fs)')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmPklEQVR4nO3de5wcZZ3v8c83c0lmkkmGZEImGQgJGBICch1BRF1WVm56BG8r6PG2uiwvRdf1iATdVXfZPSvinuOuyrLRZVnRA17AiIKGVUFWEUyAQBJCINyTAJkkJOQyydx+54+qCc2ke6YnMzU9Pf19v1796u6qp6p+U+nUr+p5qp5HEYGZmVWucaUOwMzMSsuJwMyswjkRmJlVOCcCM7MK50RgZlbhqksdwGA1NTXFnDlzSh2GmVlZue+++zZHxPR888ouEcyZM4fly5eXOgwzs7Ii6elC81w1ZGZW4ZwIzMwqnBOBmVmFcyIwM6twTgRmZhXOicDMrMI5EZiZVTgnAjOzCudEYGZW4ZwIzMwqnBOBmVmFcyIwM6twTgRmZhXOicDMrMJllggkXStpk6RV/ZQ5XdIKSasl/SarWMzMrLAsrwiuA84uNFNSI3A18LaIOBp4d4axmJlZAZklgoi4C9jaT5H3AjdHxDNp+U1ZxWJmZoWVso3gSOAgSXdKuk/SBwoVlHSRpOWSlre1tY1giGZmY18pE0E1cBLwFuAs4G8kHZmvYEQsjojWiGidPj3vkJtmZnaASjlm8Xpgc0TsAnZJugs4Dni0hDGZmVWcUl4R/AR4g6RqSfXAKcCaEsZjZlaRMrsikHQDcDrQJGk98EWgBiAiromINZJ+ATwE9ADfjoiCt5qamVk2MksEEXFhEWWuAq7KKgYzMxuYnyw2M6twTgRmZhXOicDMrMI5EZiZVTgnAjOzCudEYGZW4ZwIzMwqnBOBmVmFcyIwM6twTgRmZhXOicDMrMI5EZiZVTgnAjOzCudEYGZW4ZwIzMwqnBOBmVmFyywRSLpW0iZJ/Y46Juk1krolvSurWMzMrLAsrwiuA87ur4CkKuBKYGmGcZiZWT8ySwQRcRewdYBinwBuAjZlFYeZmfWvZG0EklqAtwPXFFH2IknLJS1va2vLPjgzswpSysbirwGXRUT3QAUjYnFEtEZE6/Tp07OPzMysglSXcNutwI2SAJqAcyV1RcSSEsZkZlZxSpYIImJu72dJ1wE/cxIwMxt5mSUCSTcApwNNktYDXwRqACJiwHYBMzMbGZklgoi4cBBlP5RVHGZm1j8/WWxmVuGcCMzMKpwTgZlZhXMiMDOrcE4EZmYVzonAzKzCORGYmVU4JwIzswrnRGBmVuGcCMzMKpwTgZlZhXMiMDOrcE4EZmYVzonAzKzCORGYmVU4JwIzswqXWSKQdK2kTZJWFZj/PkkPpa+7JR2XVSxmZlZYllcE1wFn9zP/SeCPIuJY4ApgcYaxmJlZAQMOVSnpVOB/Am8AZgLtwCrgVuC7EbE933IRcZekOYXWGxF353y9Bzik+LDNzGy49HtFIOnnwEeBpSRn9zOBhcBfAxOAn0h62zDE8RHg5/3EcZGk5ZKWt7W1DcPmzMyslyKi8EypKSI297uCfsqkVwQ/i4hj+ln+j4GrgddHxJaBAm5tbY3ly5cPVMzMzHJIui8iWvPN67dqqO8BXtLk3GUiYutAiWKAwI4Fvg2cU0wSMDOz4TdgGwGApL8A/o6kfaD3EiKAww90w5JmAzcD74+IRw90PWZmNjRFJQLgM8DRgzn7l3QDcDrQJGk98EWgBiAirgG+AEwDrpYE0FXossXMzLJTbCJ4HNg9mBVHxIUDzP8oSUO0mZmVULGJ4HLgbkn3Ant7J0bEJzOJyszMRkyxieDfgF8DK4Ge7MIxM7ORVmwi6IqIT2caiZmZlUSxXUzckT7UNVPS1N5XppGZmdmIKPaK4L3p++U504Z0+6iZmY0ORSWCiJibdSBmZlYaxV4RIOl1wBxe+WTxdzKIyczMRlCxTxZfDxwBrAC608kBOBGYmZW5Yq8IWoGF0V8PdWZmVpaKvWtoFdCcZSBmZlYaxV4RNAEPS/oDr3yyeDjGIjAzsxIqNhF8KcsgzMysdPpNBJIUid8MVGb4Q7PhsuSBDVy1dC0bt7Uzq7GOS8+az/kntIz4Ooai1Ns3G8sGuiK4Q9JNwE8i4pneiZJqgdcDHwTuIBmo3kahJQ9s4PKbV9LemdzstWFbO5ffvBKg6APpcKxjKEq9fbOxbqChKicAfwa8D5gLbCMZq7gKuB34ZkSsyDzKHJU4VOWBnA3v6ezmqS27eO+37mHrrs795o+vHsfJc6ciCQHjBJIYJwAhpdMQdz66iT2d+/c1OHViLd/+YCtNE8fT1FBLfW3h84qB/oaenmDr7g6e376HF17aw/Mv7eGF7cn7T1ZsZG/X/ttvaazjd4ve1O9+MLNEf0NV9psI+qykhqTRuD0ithVR/lrgrcCmfGMWKxmN5p+Bc0nGOvhQRNw/0HpLkQiGWi0xlOX7ng0D1NVU8Y/veDXnvLqZZ7fu5snNu3lq8y6e3LIred+8i+e27xlw3SfMbqQngAgC6IkgAnoCen8XPRE8+sLOomKtr62iadJ4pk2qpWnS+PRVy4Zt7fz0wY10dr/8W6saJ149azJVVeN4fvseNu3Y84r5kCSipknj2bRjb99N7XPF+cdw5sIZzJg8oagYzSrVsCSCA9joG4GdwHcKJIJzgU+QJIJTgH+OiFMGWu9IJ4L+DsTFHMyLWb67J9jT2c3ujm7aO7rZ3dlFe/r5khseYOuujv3WWzVORERyIE8dVF/DnKaJzJ02kTlNyeuKnz1MW54D6WDOpk/78q/ZsK19v+nTG8bzlXceS9vOvWzeuZfNOzrYsuvlz5t37mXr7g4K/cSqxomT50ylecoEZkyeQPPk8S9/njKB6ZPGU101ruD2q8aJ7nQHnDC7kbOObubMhTM4fPqkov4us0pSkkSQbngO8LMCieDfgDsj4ob0+1rg9Ih4rr91HkgiGMwZeUSwvb2TDdva2bhtD5/54YNsb9+/aqW2ahzHH9q4/wr0yq8rntlGR/f+1RpVgsl1Nezu6M5b7VGMT54xj7lN9cxtmsTcaROZUl+zX5mhJrKhrqO7J3jV524j369MwJNffssBb/9/v/0Yjm6ZwtJVz3P7wy+wcsN2AOYdPClJCkfP4NUtU5DkxmareP0lgqL7GspAC/Bszvf16bR+E8Fg5WtovOymh1i3aQeHT5/Exm3tbNi2J31vZ+O2dnZ3dA+wVujo7qFq3CuP+pHncJcvCQB0B7z12FnU11ZRV1uVvNdUUVdbnfO5ik/c8EDBM/pPv/nIAePsPdgN5SA4lHVUjROzGuvyntHPaqwblu0fOaOBT5wxjw3b2rl99fMsXf08V9+5jm/csY5ZUyZwxMGTuPfJrXSkCdeNzWavNJg2gsOAeRHxS0l1QHVE7BhgmTkUviK4FfjHiPht+v1XwGcj4r48ZS8CLgKYPXv2SU8//XRRMUPhao1cTZNqmdVYx6wpdcl74wRaGpPPf3H9fTz/0v717cVWrRTafrHLD8cZfamV4m/YuquDX615gaWrX+CXa17IW+ag+hr+48Mnc9jUehrra0iarfLzFYWVuyFfEUj6c5ID8VSSzucOAa4BzhhCXOuBQ3O+HwJszFcwIhYDiyGpGhrMRjb2kwTu+MzpzJwygQk1VQXLLDpnQd6D2KVnzS9q+5eeNX9Iyw/HGX2pleJvmDqxlne3Hsq7Ww9l7qJb81ZNvbi7k/O/+TsAGiZUc9i0eg6bOpHZ0+qZM62e2VMncti0eu59fAufW7LKt6/amFVs1dDHgZOBewEi4jFJBw9x27cAl0i6kaSxePtA7QMHolC1REtjHXObJg64/FAPYsNVNVPuB5xS/g2FfgMHN4zn788/hme27ubpLbt5eutuVm/cztLVz9PV0//5RntnN1ctXVv2/y5mUHwi2BsRHb2XzpKqIe9J1j6SbgBOB5okrQe+CNQARMQ1wG0kdwytI7l99MMHEP+AhnpGDkM/iI2FA3k5K/Qb+Ny5R3Hm0fv3pdjV3cPGbXt4eusunt6ym79esirvevu72jQrJ8Umgt9I+hxQJ+nNwMeAn/a3QERcOMD8ILnSyNRYqFqxoRnsb6C6ahyzp9Uze1o9b5gH/3rn40Nq7DYb7YpqLJY0DvgIcCbJXX9LgW+Xoo+hSnyy2EprLDTYmw3H7aN1wLUR8a10hVXptN3DE6LZ6NV7sP/8kpXs2tvNrMYJfPasBU4CNmYUOzDNr0gO/L3qgF8Ofzhmo9P5J7Tw129ZCMCNf36qk4CNKcUmggkRsa/DmfRzfTYhmY1O85sbAHjk+ZdKHInZ8Co2EeySdGLvF0knAb5lwirKkTOSRLD2+X6fozQrO8W2EXwK+KGk3ge+ZgLvySQis1Fq0vhqDp1axyMvOBHY2FJUIoiIZZIWAPNJ7hp6JCL274nNbIxb0DyZR55z1ZCNLYPpdO41wJx0mRMkERHfySQqs1FqQXMDv35kE3s6u/vtmsSsnBTb19D1JH0MrQB6b6YOwInAKsr85ga6e4J1m3ZyTMuUUodjNiyKvSJoBRZ6kHqrdAuaX24wdiKwsaLYu4ZWAft3ymJWYeZMm0ht9TjWusHYxpBirwiagIcl/QHYN0pKRLwtk6jMRqnqqnHMO3gSa9xgbGNIsYngS1kGYVZO5jc38NvHNpc6DLNhU+zto7/JOhCzcrGguYGb79/Ai7s6OGhibanDMRuyotoIJL1W0jJJOyV1SOqW5Gtjq0jzmycD8IifMLYxotjG4m8AFwKPkXQ499F0mlnFefnOIZ8L2dhQbCIgItYBVRHRHRH/QTL6WL8knS1praR1khblmT9F0k8lPShptaRMRikzG04HN4ynsb7GVwQ2ZhTbWLxbUi2wQtJXgOeAfgf8Tccs+CbwZpKB6pdJuiUiHs4p9nHg4Yj4H5KmA2slfS8iOgb9l5iNEEksaG5wIrAxo9grgvenZS8BdgGHAu8YYJmTgXUR8UR6YL8ROK9PmQAalAyGPAnYCnQVGZNZySxonsyjL+ygZ4BB7s3KQbGJ4PyI2BMRL0XE30bEp4G3DrBMC/Bszvf16bRc3wCOAjYCK4G/jIieviuSdJGk5ZKWt7W1FRmyWXbmNzewu6Ob9S+6N3Yrf8Umgg/mmfahAZZRnml9T5/OIum/aBZwPPANSZP3WyhicUS0RkTr9OnTBwzWLGsepMbGkn4TgaQLJf0UmCvplpzXncCWAda9nqQKqdchJGf+uT4M3ByJdcCTwIJB/QVmJeBBamwsGaix+G6ShuEm4J9ypu8AHhpg2WXAPElzgQ3ABcB7+5R5BjgD+G9JM0jGO3iiuNDNSmfS+GpmT613g7GNCf0mgoh4Gnha0p8A7RHRI+lIkrP2lQMs2yXpEmApUAVcGxGrJV2czr8GuAK4TtJKkqqkyyLCz+5bWZjf3OCqIRsTir199C7gDZIOAn4FLCcZqvJ9/S0UEbcBt/WZdk3O543AmYMJ2Gy08CA1NlYU21isiNhNcsvo1yPi7cDC7MIyG/1yB6kxK2dFJwJJp5JcAdyaThvMMJdmY07uIDVm5azYRPAp4HLgx2k9/+HAHZlFZVYGegepcTuBlbvBdEP9m5zvTwCfzCoos3LQO0iN7xyyctdvIpD0tYj4VPoswX7P0nuEMqt0HqTGxoKBrgiuT9+/mnUgZuXIg9TYWDDQcwT3pe8eocwsj9xBak49YlqJozE7MANVDa0kT5VQr4g4dtgjMisjC3L6HHIisHI1UNVQbw+jH0/fe6uK3gfsziQiszJycMN4Dqqv8S2kVtaK6WICSadFxGk5sxZJ+h3wd1kGZzbaSUq7mnAisPJV7HMEEyW9vveLpNcxwAhlZpXCg9RYuSv26eCPANdKmkLSZrAd+LPMojIrI7mD1MyeVl/qcMwGrdgHyu4DjksHjVFEbM82LLPykTtIjROBlaNiq4YASIeqdBIwyzF/Rm8icDuBladBJQIz29/EdJAa3zlk5cqJwGwYeJAaK2dFJwJJr5P0Xkkf6H0VsczZktZKWidpUYEyp0taIWm1JD/BbGVpQXMDT23ZzZ7O7lKHYjZoRTUWS7oeOAJYAfT+0gP4Tj/LVAHfBN5MMpD9Mkm3RMTDOWUagauBsyPiGUkHH8DfYFZyuYPUHNMypdThmA1KsbePtgILI2IwN0qfDKxLu6xG0o3AecDDOWXeC9wcEc8ARMSmQazfbNR4uauJHU4EVnaKrRpaBTQPct0twLM539en03IdCRwk6U5J9xWqbpJ0kaTlkpa3tbUNMgyz7PUOUrPW7QRWhoq9ImgCHpb0B2Bv78QBxiNQnml9ryiqgZOAM4A64PeS7omIR1+xUMRiYDFAa2urH9+0UceD1Fg5KzYRfOkA1r0eODTn+yHAxjxlNkfELmCXpLuA44BHMSszHqTGylVRVUPpeASPAA3pa00RYxQsA+ZJmiupFrgAuKVPmZ8Ab5BULakeOAVYM5g/wGy0WNDcwKYde3lxV0epQzEblKISgaQ/Bf4AvBv4U+BeSe/qb5mI6AIuAZaSHNx/kA58f7Gki9Mya4BfAA+l6/92RKw60D/GrJRyB6kxKyfFVg19HnhN7109kqYDvwR+1N9CEXEbcFufadf0+X4VcFWxAZuNVkd5kBorU8XeNTSuz62dWwaxrFlFmO5BaqxMFXtF8AtJS4Eb0u/voc+Zvlml8yA1Vq6KbSy+lOT2zWNJ7upZHBGXZRmYWTnyIDVWjoq9IiAibgJuyjAWs7LnQWqsHPWbCCTtYP+HwCB5WCwiYnImUZmVqd5BatZ4kBorI/1WDUVEQ0RMzvNqcBIw21/vIDVuMLZy4jt/zIaRB6mxcuREYDbMPEiNlRsnArNhtqC5gSc37/IgNVY2nAjMhtn85gZ6AtZt2lnqUMyK4kRgNswWuM8hKzNOBGbDbM60eg9SY2XFicBsmHmQGis3TgRmGZjf3OBbSK1sOBGYZaB3kJqtHqTGyoATgVkGXh6kxu0ENvplmggknS1praR1khb1U+41kroHGvXMrFz0DlLj6iErB5klAklVwDeBc4CFwIWSFhYodyXJkJZmY4IHqbFykuUVwcnAuoh4IiI6gBuB8/KU+wRJ99ab8swzK0sepMbKSZaJoAV4Nuf7+nTaPpJagLcDrxjHuC9JF0laLml5W1vbsAdqlgUPUmPlIstEoDzT+v6P+BpwWUT02ylLRCyOiNaIaJ0+ffpwxWeWqd5Bap59cXepQzHrV9EjlB2A9cChOd8PATb2KdMK3CgJoAk4V1JXRCzJMC6zEbEgbTB+5PkdHDZtYomjMSssyyuCZcA8SXMl1QIXALfkFoiIuRExJyLmAD8CPuYkYGPFkR6kxspEZlcEEdEl6RKSu4GqgGsjYrWki9P5/bYLmJU7D1Jj5SLLqiEi4jbgtj7T8iaAiPhQlrGYlYIHqbFy4CeLzTLkQWqsHDgRmGVoQfNkD1Jjo54TgVmG5ufcOWQ2WmXaRmBW6TxIzfBY8sAGrlq6lo3b2pnVWMelZ83n/BNaBl7QiuJEYJYhD1IzdEse2MDlN6+kPW1n2bCtnctvXglQdDJwIumfE4FZxuY3N/DbxzYPaR3lfiArNv6u7h42bGvnibZdPN62kyc37+KH962no6vnFeXaO7u59EcPctP965lSV8NB9bU01te84nPyvZZ7ntjM39+6hj2dyToOJJGMdU4EZhlb0NzAzfdvYOuuDqZOrB308sNxRlxK+eJfdPNDPLF5J4c01vP45p082baLJzbv4pktu+nofvmgP6WuZr8k0KuzO9ixp4v1L7azbXcH29s7KbZbp/bObq78xSNlsf8g+xMBJwKzjC3IGaTmdUc0DWrZPZ3d/O1PV+87iPZq7+zmi7esonnKBBY0N9BYP/gEMxJ27Onk7299eL/493T28C+/WgdATZU4bNpEDm+ayBlHHcwRTZM4fPpE5jZNZOrEWl5/5R1s2Na+37pbGutY8vHT9n3v6UkSw7b2Drbt7uTFNDn85Y0r8sb23PY9nP21u2idcxCvmTOV1jlTaWmsG74/fpiMxImAE4FZxhbkDFJTbCJ4eONLfH/ZMyxZsZHt7Z15y2xv7+KCxfcAMHPKBI6aOZkFzQ0smDmZo5obmNs0keqq5MbAoZ5RFrN82469rN64ndUbX+LhjS+xeuN2ntrSf4d7v7n0dFoa6/bFmc+lZ81/xYEQoK6mikvPmv+KcuPGiSn1NUypr+GwaS9P/8ov1uZNJA0Tqjl48gSWPLCR797zDACzpkygdc5UXjPnIFrnTGX+jAbGjdOI7L9e3T3BM1t389gLO3hs006+8et1eU8Erlq61onArFwUO0jN9vZObnlwIz9Y9iwrN2yntnocZx/dzN2Pb2bzzv3HPm6ePIEr33Usjzz3Eo88v4M1z73EXY+20ZXWj9RWj+PIGZOor6nigWe30dmdTN+wrZ1FNz3E7o4u3nXSodRUibTjx7zynZFedtNDPLj+RSaNr2HVhuTgv2nH3n3LHDq1jqNnTuGdJx7CdXc/xZY8Yze3NNYV1Rlf78HuQA/EhRLJFecdw/kntNDV3cMjz+9g+VNbWfb0i9zzxBZueTDpH7NhQjWHNNbx2Kad+/Zrb9VWbmz9KXRG39PTw3GzD+KxF3buO+g/tmknj7ftLFgdlmtjnuR2oBRRXn2lt7a2xvLly0sdhtmgXLD49+zp7HlFVQZARHDPE1v5wfJnuW3lc+zt6uGomZN5T+shnH9CC431tfsdSCA5kP3jO16934Goo6uHx9t2siYnOfxu3eZ+686rxon6miom1FZRV1NFfW0VE9L3upoq7n58M+2d+Q9M4wSvOngSR8+awtGzJnP0rCksnDWZKXU1+8oMJv6sDOaMPCJY/2I7y57ayrKnXuSHy5/dlwT6aqyvYWJtNZPGV1M/vopJ46uZWFvNxPHVTBpfxcTx1Xz3nqd5aU/XgDG2NNYxb8YkjpzRwKsOnsS8gyfxqoMncfbX/rtg1djvFr2p6H0g6b6IaM07z4nALHvv//d79905NKuxjoveOJede7v54fJneWrLbhrGV/O242dxwWtmc0zL5P3O0IdSNTF30a37DQTS6zNnHkl7ZzftHT20d3bR3tHN7o7udFryvnpj/mcgBKy54mwm1FQNGEM53/XU3/57/2sPY1dHF7v2drFrbzc79/Z+7ko+d3TT3U8W/uq7j9t3wJ84Pn8FzXAl0v4SgauGzDK25IEN3PPEln0Hkw3b2vniLQ8DcMrcqXzyjHmcc8xM6moLH1DPP6HlgA+csxrrCp5RXvKmeQMuf9qXf513+VmNdUUlARha/KXW3/674vxj+l02Ijjty79m4/Y9eZd/10mHDLj9oVaNFcNdTJhl7Kqla/fVz+eaMXk83/+LU3nHiYf0mwSG6tKz5lPX54Cdr7E1q+XL3VD+fkl89uwFQ95/55/Qwu8WvYknv/wWfrfoTcOeVH1FYJaxQo16m17am3f6cBvqGeVInJGOZpWw/9xGYJaxQlUrg23sMxuK/toIMq0aknS2pLWS1klalGf++yQ9lL7ulnRclvGYlUKlV63Y6JdZ1ZCkKuCbwJtJBrJfJumWiHg4p9iTwB9FxIuSzgEWA6dkFZNZKZRD1YBVtizbCE4G1kXEEwCSbgTOA/Ylgoi4O6f8PcDATehmZaic75qxsS/LqqEW4Nmc7+vTaYV8BPh5vhmSLpK0XNLytra2YQzRzMyyTAT5nlnP2zIt6Y9JEsFl+eZHxOKIaI2I1unTpw9jiGZmlmXV0Hrg0JzvhwAb+xaSdCzwbeCciNiSYTxmZpZHllcEy4B5kuZKqgUuAG7JLSBpNnAz8P6IeDTDWMzMrIDMrggiokvSJcBSoAq4NiJWS7o4nX8N8AVgGnB12rdKV6H7XM3MLBt+oMzMrAKU7IEyMzMb/ZwIzMwqnBOBmVmFcyIwM6twTgRmZhXOicDMrMI5EZiZVTgnAjOzCudEYGZW4ZwIzMwqnBOBmVmFcyIwM6twTgRmZhXOicDMrMI5EZiZVTgnAjOzCpdpIpB0tqS1ktZJWpRnviT9Szr/IUknZhmPmZntL7NEIKkK+CZwDrAQuFDSwj7FzgHmpa+LgH/NKh4zM8svyyuCk4F1EfFERHQANwLn9SlzHvCdSNwDNEqamWFMZmbWR2aD1wMtwLM539cDpxRRpgV4LreQpItIrhgAdkpaO7yhDpsmYHOpg+jHaI8PRn+Mjm9oHN/QDCW+wwrNyDIRKM+0OIAyRMRiYPFwBJUlScsLDQ49Goz2+GD0x+j4hsbxDU1W8WVZNbQeODTn+yHAxgMoY2ZmGcoyESwD5kmaK6kWuAC4pU+ZW4APpHcPvRbYHhHP9V2RmZllJ7OqoYjoknQJsBSoAq6NiNWSLk7nXwPcBpwLrAN2Ax/OKp4RMtqrr0Z7fDD6Y3R8Q+P4hiaT+BSxX5W8mZlVED9ZbGZW4ZwIzMwqnBPBIEk6VNIdktZIWi3pL/OUOV3Sdkkr0tcXRjjGpyStTLe9PM/8knXtIWl+zn5ZIeklSZ/qU2bE95+kayVtkrQqZ9pUSf8l6bH0/aACy/bblUqG8V0l6ZH03/DHkhoLLNvv7yHD+L4kaUPOv+O5BZYt1f77fk5sT0laUWDZTPdfoWPKiP7+IsKvQbyAmcCJ6ecG4FFgYZ8ypwM/K2GMTwFN/cw/F/g5yXMcrwXuLVGcVcDzwGGl3n/AG4ETgVU5074CLEo/LwKuLPA3PA4cDtQCD/b9PWQY35lAdfr5ynzxFfN7yDC+LwGfKeI3UJL912f+PwFfKMX+K3RMGcnfn68IBikinouI+9PPO4A1JE9Dl5PR0rXHGcDjEfF0Cbb9ChFxF7C1z+TzgP9MP/8ncH6eRYvpSiWT+CLi9ojoSr/eQ/IcTkkU2H/FKNn+6yVJwJ8CNwz3dovRzzFlxH5/TgRDIGkOcAJwb57Zp0p6UNLPJR09spERwO2S7ku75+irUNceI+0CCv/nK+X+6zUj0uda0veD85QZLfvyz0iu8vIZ6PeQpUvSqqtrC1RtjIb99wbghYh4rMD8Edt/fY4pI/b7cyI4QJImATcBn4qIl/rMvp+kuuM44OvAkhEO77SIOJGkd9ePS3pjn/lFde2RJSUPGb4N+GGe2aXef4MxGvbl54Eu4HsFigz0e8jKvwJHAMeT9B/2T3nKlHz/ARfS/9XAiOy/AY4pBRfLM23Q+8+J4ABIqiH5B/teRNzcd35EvBQRO9PPtwE1kppGKr6I2Ji+bwJ+THL5mGs0dO1xDnB/RLzQd0ap91+OF3qrzNL3TXnKlHRfSvog8FbgfZFWGvdVxO8hExHxQkR0R0QP8K0C2y31/qsG3gF8v1CZkdh/BY4pI/b7cyIYpLQ+8d+BNRHxfwqUaU7LIelkkv28ZYTimyipofczSYPiqj7FRkPXHgXPwkq5//q4Bfhg+vmDwE/ylCmmK5VMSDobuAx4W0TsLlCmmN9DVvHltju9vcB2S7b/Un8CPBIR6/PNHIn9188xZeR+f1m1hI/VF/B6kkuvh4AV6etc4GLg4rTMJcBqkhb8e4DXjWB8h6fbfTCN4fPp9Nz4RDJo0OPASqB1hPdhPcmBfUrOtJLuP5Kk9BzQSXKW9RFgGvAr4LH0fWpadhZwW86y55Lc6fF47/4eofjWkdQP9/4Or+kbX6HfwwjFd336+3qI5OA0czTtv3T6db2/u5yyI7r/+jmmjNjvz11MmJlVOFcNmZlVOCcCM7MK50RgZlbhnAjMzCqcE4GZWYVzIrBRQ9K0nN4gn8/puXKnpKsz3O7pkl6X1foLbPNrvU+oSvpe2nvkqrQrhpp0ulSgl9hiepyUNF3SvZIekPSGAmW+KulNWfyNVj6cCGzUiIgtEXF8RBwPXAP83/T7pIj4WIabPh0YsUQgaSrw2kg6QoOka4gFwKuBOuCj6fRzgHnp6yKSLhuQVEXyHMg5JL1UXihpYZ5NnUHysNQJEfHfBcL5OknPllbBnAhs1EvP2H+Wfv6SpP+UdLuSfuLfIekrSvqL/0XO2fRJkn6TdhS2NOdR/U9Kejg9w74x7eTrYuCv0quPN6Rn0jdJWpa+TsvZ9vWSfq2kj/g/T6fPlHRXuvyqQmffOd4F/KL3S0TcFingD7zci2ihXmIH7HFS0vEk3Rifm8Y1UdJ1aXwrJf1Vuu2ngWmSmg/k38bGhswGrzfL0BHAH5OcDf8eeGdEfFbSj4G3SLqV5Ez3vIhok/Qe4B9IeuhcBMyNiL2SGiNim6RrgJ0R8VUASf+P5Grkt5JmA0uBo9JtH0syhsNE4IF0WxcCSyPiH9Kz9foB4j8N+FHfiWkSez/QO9hRoZ4l800/JXddEdE7oE9rRFwi6SSgJSKOSbfVmFP8/jSmmwaI28YoJwIrRz+PiE5JK0kG5ug9u14JzAHmA8cA/5V2WVRF0r0AJI/xf0/SEgr3avonwMJ0WYDJvf3NAD+JiHagXdIdJGfny4Deuv0lEbFigPhnAm15pl8N3JVTjVOoZ8kD6XHyCeBwSV8HbgVuz5m3iaTbAqtQrhqycrQXIJJeLTvj5X5SekhObgSs7m1viIhXR8SZaZm3kNSvnwTcp6T3yb7GAafmLN8SyYAhsP8BN9K6/jcCG4DrJX1ggPjbgQm5EyR9EZgOfDpncqGeJQfd42REvAgcB9wJfBz4ds7sCWlMVqGcCGwsWgtMl3QqJFUuko6WNA44NCLuAD4LNAKTgB0kQwT2up2k4zvS5Y/PmXeepAmSppE0Mi+TdBiwKSK+RdKL5Inpct9R0ntqX2uAV+Ws/6PAWcCFaXLrVaiX2EH3OKmkG+9xEXET8De9MaaOZIR6JLXRyYnAxpy0AfVdwJWSHiTpzfF1JFVE302rlB4gaQfYBvwUeHtvYzHwSaA1bVB+mKQxudcfSKpW7gGuiKSv+tOBFZIeAN4J/HNa9lherpLKdWu6TK9rgBnA79MYvpBOv42kSmcdSX/+H0v/vi6SRLWUJKn8ICJWD7BbWoA7lQzQfh1wOexrl3gVkMmg9lYe3PuoWZEkfYmcRuUByk4G/j0i3l1g/m+Bt6aJqGQkvZ1k4PS/KWUcVlq+IjDLQCSjrOVNAqn/BcweqXj6UU3+ISStgviKwMyswvmKwMyswjkRmJlVOCcCM7MK50RgZlbhnAjMzCrc/wfh5+yUh1dSFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(1, 21), distances, '-o')\n",
    "plt.ylim(0,1.75)\n",
    "plt.ylabel('Ion distance (nm)')\n",
    "plt.xlabel(f'Timesteps, ({TIMESTEP*nsteps})')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c754cfd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}