long range correction experiments

9_debug_long_correction_slowdown.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8c863835",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.\n"
     ]
    }
   ],
   "source": [
    "from simtk import openmm\n",
    "from simtk.openmm import app, unit\n",
    "\n",
    "import time\n",
    "import random\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2d34af07",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_2_particle_system_and_positions(expression):\n",
    "    \n",
    "    # Create system and add particles\n",
    "    system = openmm.System()\n",
    "    system.addParticle(22.99)\n",
    "    system.addParticle(35.45)\n",
    "    system.setDefaultPeriodicBoxVectors(openmm.Vec3(9, 0, 0), openmm.Vec3(0, 9, 0), openmm.Vec3(0, 0, 9))\n",
    "    \n",
    "    # Create CustomNonbondedForce (will add particles to the force later)\n",
    "    nb_force = openmm.CustomNonbondedForce(expression)\n",
    "    system.addForce(nb_force)\n",
    "    nb_force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)\n",
    "    nb_force.setCutoffDistance(0.1 * unit.nanometer)\n",
    "    \n",
    "    # Create positions\n",
    "    positions = [openmm.Vec3(random.uniform(0, 9)*0.1, 0, 0) for _ in range(2)]\n",
    "    \n",
    "    return system, positions\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "59b1bf0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_experiment(system, positions):\n",
    "    # Set up context\n",
    "    integrator = openmm.LangevinMiddleIntegrator(300, 1.0, 0.004)\n",
    "    platform = openmm.Platform.getPlatformByName(\"CUDA\")\n",
    "    context = openmm.Context(system, integrator, platform)\n",
    "    context.setPositions(positions)\n",
    "\n",
    "    # Time getState()\n",
    "    force_group_idx = 0\n",
    "    times = []\n",
    "    for i in range(5):\n",
    "        context.setParameter('lambda', 1 - i*0.1)\n",
    "        initial_time = time.time()\n",
    "        context.getState(getEnergy=True, groups=2**force_group_idx)\n",
    "        elapsed_time = (time.time() - initial_time)\n",
    "        print(elapsed_time)\n",
    "        times.append(elapsed_time)\n",
    "    print(\"mean: \", np.mean(times[1:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56d617d7",
   "metadata": {},
   "source": [
    "# With simple expression and no per particle params, long range correction off"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b1dd56e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_expression = \"U_sterics; U_sterics = lambda * 4 * epsilon * x * (x - 1.0); x = (sigma / r_eff_sterics)^6; sigma = 1; epsilon = 1; r_eff_sterics = sqrt(r^2 + w_sterics^2); w_sterics = 0.3;\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2cb4f275",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0013477802276611328\n",
      "0.00010633468627929688\n",
      "9.799003601074219e-05\n",
      "8.940696716308594e-05\n",
      "8.630752563476562e-05\n",
      "mean:  9.500980377197266e-05\n"
     ]
    }
   ],
   "source": [
    "# Create system and customize it for this experiment\n",
    "system, positions = get_2_particle_system_and_positions(simple_expression)\n",
    "nb_force = system.getForce(0)\n",
    "nb_force.setUseLongRangeCorrection(False)\n",
    "nb_force.addGlobalParameter('lambda', 1.0)\n",
    "nb_force.addParticle([])\n",
    "nb_force.addParticle([])\n",
    "\n",
    "# Run experiment\n",
    "run_experiment(system, positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3766834d",
   "metadata": {},
   "source": [
    "# With simple expression and no per particle params, long range correction on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "af938912",
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_expression = \"U_sterics; U_sterics = lambda * 4 * epsilon * x * (x - 1.0); x = (sigma / r_eff_sterics)^6; sigma = 1; epsilon = 1; r_eff_sterics = sqrt(r^2 + w_sterics^2); w_sterics = 0.3;\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cb31f950",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.005000591278076172\n",
      "0.0017299652099609375\n",
      "0.0016379356384277344\n",
      "0.001554727554321289\n",
      "0.0016150474548339844\n",
      "mean:  0.0016344189643859863\n"
     ]
    }
   ],
   "source": [
    "# Create system and customize it for this experiment\n",
    "system, positions = get_2_particle_system_and_positions(simple_expression)\n",
    "nb_force = system.getForce(0)\n",
    "nb_force.setUseLongRangeCorrection(True)\n",
    "nb_force.addGlobalParameter('lambda', 1.0)\n",
    "nb_force.addParticle([])\n",
    "nb_force.addParticle([])\n",
    "\n",
    "# Run experiment\n",
    "run_experiment(system, positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a49d9b13",
   "metadata": {},
   "source": [
    "# With simple expression and many per particle params, long range correction off"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ea41e69c",
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_expression = \"U_sterics; U_sterics = lambda * 4 * epsilon * x * (x - 1.0); x = (sigma / r_eff_sterics)^6; sigma = 1; epsilon = 1; r_eff_sterics = sqrt(r^2 + w_sterics^2); w_sterics = 0.3;\"\n",
    "# simple_expression = \"U_sterics; U_sterics = lambda * 4 * epsilon * x * (x - 1.0); x = (sigma / r_eff_sterics)^6; sigma = sigma_old1; epsilon = epsilon_old1 * is_rest * (is_unique_old + is_core + is_unique_new); r_eff_sterics = sqrt(r^2 + w_sterics^2); w_sterics = 0.3;\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9355261c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.001435995101928711\n",
      "0.00011420249938964844\n",
      "9.5367431640625e-05\n",
      "9.560585021972656e-05\n",
      "9.608268737792969e-05\n",
      "mean:  0.00010031461715698242\n"
     ]
    }
   ],
   "source": [
    "# Create system and customize it for this experiment\n",
    "system, positions = get_2_particle_system_and_positions(simple_expression)\n",
    "nb_force = system.getForce(0)\n",
    "nb_force.setUseLongRangeCorrection(False)\n",
    "nb_force.addGlobalParameter('lambda', 1.0)\n",
    "\n",
    "# Add per-particle parameters for rest scaling -- these three sets are disjoint\n",
    "nb_force.addPerParticleParameter(\"is_rest\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solute\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solvent\")\n",
    "\n",
    "# Add per-particle parameters for alchemical scaling -- these sets are also disjoint\n",
    "nb_force.addPerParticleParameter('is_environment')\n",
    "nb_force.addPerParticleParameter('is_core')\n",
    "nb_force.addPerParticleParameter('is_unique_old')\n",
    "nb_force.addPerParticleParameter('is_unique_new')\n",
    "\n",
    "# Add per-particle parameters for defining energy\n",
    "nb_force.addPerParticleParameter(\"sigma_old\")\n",
    "nb_force.addPerParticleParameter(\"sigma_new\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_old\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_new\")\n",
    "\n",
    "nb_force.addParticle([1, 0, 0, 0, 0, 1, 0, 0.24, 0.24, 0.37, 0.37])\n",
    "nb_force.addParticle([1, 0, 0, 0, 1, 0, 0, 0.45, 0.45, 0.15, 0.15])\n",
    "\n",
    "\n",
    "# Run experiment\n",
    "run_experiment(system, positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddb3d407",
   "metadata": {},
   "source": [
    "# With simple expression and many per particle params, long range correction on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ce2726ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_expression = \"U_sterics; U_sterics = lambda * 4 * epsilon * x * (x - 1.0); x = (sigma / r_eff_sterics)^6; sigma = 1; epsilon = 1; r_eff_sterics = sqrt(r^2 + w_sterics^2); w_sterics = 0.3;\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "46ed7452",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0053102970123291016\n",
      "0.0018839836120605469\n",
      "0.001558065414428711\n",
      "0.0017595291137695312\n",
      "0.0015594959259033203\n",
      "mean:  0.0016902685165405273\n"
     ]
    }
   ],
   "source": [
    "# Create system and customize it for this experiment\n",
    "system, positions = get_2_particle_system_and_positions(simple_expression)\n",
    "nb_force = system.getForce(0)\n",
    "nb_force.setUseLongRangeCorrection(True)\n",
    "nb_force.addGlobalParameter('lambda', 1.0)\n",
    "\n",
    "# Add per-particle parameters for rest scaling -- these three sets are disjoint\n",
    "nb_force.addPerParticleParameter(\"is_rest\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solute\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solvent\")\n",
    "\n",
    "# Add per-particle parameters for alchemical scaling -- these sets are also disjoint\n",
    "nb_force.addPerParticleParameter('is_environment')\n",
    "nb_force.addPerParticleParameter('is_core')\n",
    "nb_force.addPerParticleParameter('is_unique_old')\n",
    "nb_force.addPerParticleParameter('is_unique_new')\n",
    "\n",
    "# Add per-particle parameters for defining energy\n",
    "nb_force.addPerParticleParameter(\"sigma_old\")\n",
    "nb_force.addPerParticleParameter(\"sigma_new\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_old\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_new\")\n",
    "\n",
    "nb_force.addParticle([1, 0, 0, 0, 0, 1, 0, 0.24, 0.24, 0.37, 0.37])\n",
    "nb_force.addParticle([1, 0, 0, 0, 1, 0, 0, 0.45, 0.45, 0.15, 0.15])\n",
    "\n",
    "\n",
    "# Run experiment\n",
    "run_experiment(system, positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e7b9f89",
   "metadata": {},
   "source": [
    "# With complex expression and many per particle params, long range correction off"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cbe4256b",
   "metadata": {},
   "outputs": [],
   "source": [
    "_default_sterics_expression_list = [\n",
    "\n",
    "        \"U_sterics;\",\n",
    "\n",
    "        # Define sterics functional form\n",
    "        \"U_sterics = 4 * epsilon * x * (x - 1.0);\"\n",
    "        \"x = (sigma / r_eff_sterics)^6;\"\n",
    "\n",
    "        # Define sigma\n",
    "        \"sigma = (sigma1 + sigma2) / 2;\",\n",
    "\n",
    "        # Define sigma1 and sigma2 (with alchemical scaling)\n",
    "        \"sigma1 = is_unique_old1 * sigma_old1 + is_unique_new1 * sigma_new1 + is_core1 * (max(0.05, lambda * sigma_old1 + lambda * sigma_new1)) + is_environment1 * sigma_old1;\",\n",
    "        \"sigma2 = is_unique_old2 * sigma_old2 + is_unique_new2 * sigma_new2 + is_core2 * (max(0.05, lambda * sigma_old2 + lambda * sigma_new2)) + is_environment2 * sigma_old2;\",\n",
    "\n",
    "        # Define epsilon (with rest scaling)\n",
    "        \"epsilon = p1_sterics_rest_scale * p2_sterics_rest_scale * sqrt(epsilon1 * epsilon2);\",\n",
    "\n",
    "        # Define epsilon1 and epsilon2 (with alchemical scaling)\n",
    "        \"epsilon1 = is_unique_old1 * old_epsilon_scaled1 + is_unique_new1 * new_epsilon_scaled1 + is_core1 * (old_epsilon_scaled1 + new_epsilon_scaled1) + is_environment1 * epsilon_old1;\",\n",
    "        \"old_epsilon_scaled1 = lambda * epsilon_old1;\",\n",
    "        \"new_epsilon_scaled1 = lambda * epsilon_new1;\",\n",
    "\n",
    "        \"epsilon2 = is_unique_old2 * old_epsilon_scaled2 + is_unique_new2 * new_epsilon_scaled2 + is_core2 * (old_epsilon_scaled2 + new_epsilon_scaled2) + is_environment2 * epsilon_old2;\",\n",
    "        \"old_epsilon_scaled2 = lambda * epsilon_old2;\",\n",
    "        \"new_epsilon_scaled2 = lambda * epsilon_new2;\",\n",
    "\n",
    "        # Define rest scale factors (normal rest)\n",
    "        \"p1_sterics_rest_scale = is_rest1 * lambda + is_nonrest_solute1 + is_nonrest_solvent1;\",\n",
    "        \"p2_sterics_rest_scale = is_rest2 * lambda + is_nonrest_solute2 + is_nonrest_solvent2;\",\n",
    "\n",
    "        # # Define rest scale factors (scaled water rest)\n",
    "        # \"p1_sterics_rest_scale = lambda_rest_sterics * (is_rest1 + is_nonrest_solvent1 * is_rest2) + 1 * (is_nonrest_solute1 + is_nonrest_solvent1 * is_nonrest_solute2 + is_nonrest_solvent1 * is_nonrest_solvent2);\",\n",
    "        # \"p2_sterics_rest_scale = lambda_rest_sterics * (is_rest2 + is_nonrest_solvent2 * is_rest1) + 1 * (is_nonrest_solute2 + is_nonrest_solvent2 * is_nonrest_solute1 + is_nonrest_solvent2 * is_nonrest_solvent1);\",\n",
    "\n",
    "        # Define r_eff\n",
    "        \"r_eff_sterics = sqrt(r^2 + w_sterics^2);\",\n",
    "\n",
    "        # Define 4th dimension terms:\n",
    "        \"w_sterics = w_scale * r_cutoff * (is_unique_old * lambda + is_unique_new * lambda);\", # because we want w for unique old atoms to go from 0 to 1 and the opposite for unique new atoms\n",
    "        \"is_unique_old = step(is_unique_old1 + is_unique_old2 - 0.1);\", # if at least one of the particles in the interaction is unique old\n",
    "        \"is_unique_new = step(is_unique_new1 + is_unique_new2 - 0.1);\", # if at least one of the particles in the interaction is unique new\n",
    "        \"w_scale = 0.3;\",\n",
    "        \"r_cutoff = 1.0;\"\n",
    "\n",
    "    ]\n",
    "\n",
    "_default_sterics_expression = ' '.join(_default_sterics_expression_list)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "67b0ec2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.018584251403808594\n",
      "0.0002193450927734375\n",
      "0.00018906593322753906\n",
      "0.00019788742065429688\n",
      "0.0001819133758544922\n",
      "mean:  0.0001970529556274414\n"
     ]
    }
   ],
   "source": [
    "# Create system and customize it for this experiment\n",
    "system, positions = get_2_particle_system_and_positions(_default_sterics_expression)\n",
    "nb_force = system.getForce(0)\n",
    "nb_force.setUseLongRangeCorrection(False)\n",
    "\n",
    "# Add global parameters\n",
    "nb_force.addGlobalParameter(f\"lambda\", 1.0)\n",
    "\n",
    "# Add per-particle parameters for rest scaling -- these three sets are disjoint\n",
    "nb_force.addPerParticleParameter(\"is_rest\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solute\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solvent\")\n",
    "\n",
    "# Add per-particle parameters for alchemical scaling -- these sets are also disjoint\n",
    "nb_force.addPerParticleParameter('is_environment')\n",
    "nb_force.addPerParticleParameter('is_core')\n",
    "nb_force.addPerParticleParameter('is_unique_old')\n",
    "nb_force.addPerParticleParameter('is_unique_new')\n",
    "\n",
    "# Add per-particle parameters for defining energy\n",
    "nb_force.addPerParticleParameter(\"sigma_old\")\n",
    "nb_force.addPerParticleParameter(\"sigma_new\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_old\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_new\")\n",
    "\n",
    "nb_force.addParticle([1, 0, 0, 0, 0, 1, 0, 0.24, 0.24, 0.37, 0.37])\n",
    "nb_force.addParticle([1, 0, 0, 0, 1, 0, 0, 0.45, 0.45, 0.15, 0.15])\n",
    "\n",
    "# Run experiment\n",
    "run_experiment(system, positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7173a2f",
   "metadata": {},
   "source": [
    "# With complex expression and many per particle params, long range correction on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4a6fb7f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "_default_sterics_expression_list = [\n",
    "\n",
    "        \"U_sterics;\",\n",
    "\n",
    "        # Define sterics functional form\n",
    "        \"U_sterics = 4 * epsilon * x * (x - 1.0);\"\n",
    "        \"x = (sigma / r_eff_sterics)^6;\"\n",
    "\n",
    "        # Define sigma\n",
    "        \"sigma = (sigma1 + sigma2) / 2;\",\n",
    "\n",
    "        # Define sigma1 and sigma2 (with alchemical scaling)\n",
    "        \"sigma1 = is_unique_old1 * sigma_old1 + is_unique_new1 * sigma_new1 + is_core1 * (max(0.05, lambda * sigma_old1 + lambda * sigma_new1)) + is_environment1 * sigma_old1;\",\n",
    "        \"sigma2 = is_unique_old2 * sigma_old2 + is_unique_new2 * sigma_new2 + is_core2 * (max(0.05, lambda * sigma_old2 + lambda * sigma_new2)) + is_environment2 * sigma_old2;\",\n",
    "\n",
    "        # Define epsilon (with rest scaling)\n",
    "        \"epsilon = p1_sterics_rest_scale * p2_sterics_rest_scale * sqrt(epsilon1 * epsilon2);\",\n",
    "\n",
    "        # Define epsilon1 and epsilon2 (with alchemical scaling)\n",
    "        \"epsilon1 = is_unique_old1 * old_epsilon_scaled1 + is_unique_new1 * new_epsilon_scaled1 + is_core1 * (old_epsilon_scaled1 + new_epsilon_scaled1) + is_environment1 * epsilon_old1;\",\n",
    "        \"old_epsilon_scaled1 = lambda * epsilon_old1;\",\n",
    "        \"new_epsilon_scaled1 = lambda * epsilon_new1;\",\n",
    "\n",
    "        \"epsilon2 = is_unique_old2 * old_epsilon_scaled2 + is_unique_new2 * new_epsilon_scaled2 + is_core2 * (old_epsilon_scaled2 + new_epsilon_scaled2) + is_environment2 * epsilon_old2;\",\n",
    "        \"old_epsilon_scaled2 = lambda * epsilon_old2;\",\n",
    "        \"new_epsilon_scaled2 = lambda * epsilon_new2;\",\n",
    "\n",
    "        # Define rest scale factors (normal rest)\n",
    "        \"p1_sterics_rest_scale = is_rest1 * lambda + is_nonrest_solute1 + is_nonrest_solvent1;\",\n",
    "        \"p2_sterics_rest_scale = is_rest2 * lambda + is_nonrest_solute2 + is_nonrest_solvent2;\",\n",
    "\n",
    "        # # Define rest scale factors (scaled water rest)\n",
    "        # \"p1_sterics_rest_scale = lambda_rest_sterics * (is_rest1 + is_nonrest_solvent1 * is_rest2) + 1 * (is_nonrest_solute1 + is_nonrest_solvent1 * is_nonrest_solute2 + is_nonrest_solvent1 * is_nonrest_solvent2);\",\n",
    "        # \"p2_sterics_rest_scale = lambda_rest_sterics * (is_rest2 + is_nonrest_solvent2 * is_rest1) + 1 * (is_nonrest_solute2 + is_nonrest_solvent2 * is_nonrest_solute1 + is_nonrest_solvent2 * is_nonrest_solvent1);\",\n",
    "\n",
    "        # Define r_eff\n",
    "        \"r_eff_sterics = sqrt(r^2 + w_sterics^2);\",\n",
    "\n",
    "        # Define 4th dimension terms:\n",
    "        \"w_sterics = w_scale * r_cutoff * (is_unique_old * lambda + is_unique_new * lambda);\", # because we want w for unique old atoms to go from 0 to 1 and the opposite for unique new atoms\n",
    "        \"is_unique_old = step(is_unique_old1 + is_unique_old2 - 0.1);\", # if at least one of the particles in the interaction is unique old\n",
    "        \"is_unique_new = step(is_unique_new1 + is_unique_new2 - 0.1);\", # if at least one of the particles in the interaction is unique new\n",
    "        \"w_scale = 0.3;\",\n",
    "        \"r_cutoff = 1.0;\"\n",
    "\n",
    "    ]\n",
    "\n",
    "_default_sterics_expression = ' '.join(_default_sterics_expression_list)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c6d67485",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.012380361557006836\n",
      "0.006872892379760742\n",
      "0.006464958190917969\n",
      "0.006218433380126953\n",
      "0.006335735321044922\n",
      "mean:  0.0064730048179626465\n"
     ]
    }
   ],
   "source": [
    "# Create system and customize it for this experiment\n",
    "system, positions = get_2_particle_system_and_positions(_default_sterics_expression)\n",
    "nb_force = system.getForce(0)\n",
    "nb_force.setUseLongRangeCorrection(True)\n",
    "\n",
    "# Add global parameters\n",
    "nb_force.addGlobalParameter(f\"lambda\", 1.0)\n",
    "\n",
    "# Add per-particle parameters for rest scaling -- these three sets are disjoint\n",
    "nb_force.addPerParticleParameter(\"is_rest\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solute\")\n",
    "nb_force.addPerParticleParameter(\"is_nonrest_solvent\")\n",
    "\n",
    "# Add per-particle parameters for alchemical scaling -- these sets are also disjoint\n",
    "nb_force.addPerParticleParameter('is_environment')\n",
    "nb_force.addPerParticleParameter('is_core')\n",
    "nb_force.addPerParticleParameter('is_unique_old')\n",
    "nb_force.addPerParticleParameter('is_unique_new')\n",
    "\n",
    "# Add per-particle parameters for defining energy\n",
    "nb_force.addPerParticleParameter(\"sigma_old\")\n",
    "nb_force.addPerParticleParameter(\"sigma_new\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_old\")\n",
    "nb_force.addPerParticleParameter(\"epsilon_new\")\n",
    "\n",
    "nb_force.addParticle([1, 0, 0, 0, 0, 1, 0, 0.24, 0.24, 0.37, 0.37])\n",
    "nb_force.addParticle([1, 0, 0, 0, 1, 0, 0, 0.45, 0.45, 0.15, 0.15])\n",
    "\n",
    "# Run experiment\n",
    "run_experiment(system, positions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3e2e947",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62754c6e",
   "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}