{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "6b971986", "metadata": {}, "source": [ "## Alchemical free energy calculations\n", "\n", "*Computing the free energy of inserting a Lennard-Jones particle in a Lennard-Jones fluid.*\n", "\n", "This tutorial is described in [OpenMM 7 publication](http://dx.doi.org/10.1371/journal.pcbi.1005659).\n", "\n", "### Basic concept\n", "\n", "OpenMM's custom forces---which allow the programmer to express a potential algebraically, potentially with multiple parameters that can be adjusted on the fly---allow a great deal of flexibility and simplicity in encoding potentials while still achieving high performance on GPUs. One common use of this facility is to convert standard interactions (such as Lennard-Jones potentials) into alchemically-modified potentials for the purposes of computing free energy differences. The alchemical free energy code [YANK](http://github.com/choderalab/yank), for example, uses a variety of custom forces to represent alchemically-modified potentials for [protein-ligand alchemical binding free energy calculations](http://dx.doi.org/10.1007/s10822-013-9689-8).\n", "\n", "### Defining alchemically-modified potentials\n", "\n", "As a simple example of how this is facilitated by custom forces, consider computing the chemical potential of liquid argon by estimating the free energy of alchemically annihilating a Lennard-Jones particle. First, we create a simple Lennard-Jones fluid to represent liquid argon at 120 K and 80 atm, which can be conveniently done using the `testsystems` module of the conda-installable [openmmtools](http://github.com/choderalab/openmmtools) package:" ] }, { "cell_type": "code", "execution_count": null, "id": "760b92bf", "metadata": {}, "outputs": [], "source": [ "# openmmtools is also needed for this notebook\n", "!mamba install -y -c conda-forge openmmtools" ] }, { "cell_type": "code", "execution_count": null, "id": "dea68130", "metadata": {}, "outputs": [], "source": [ "\n", "from openmm.app import *\n", "from openmm import *\n", "from openmm.unit import *\n", "from openmmtools.testsystems import LennardJonesFluid\n", "\n", "\n", "# Simulation settings\n", "pressure = 80*atmospheres\n", "temperature = 120*kelvin\n", "collision_rate = 5/picoseconds\n", "timestep = 2.5*femtoseconds\n", "\n", "\n", "# Create a Lennard Jones test fluid\n", "sigma = 3.4*angstrom\n", "epsilon = 0.238 * kilocalories_per_mole\n", "fluid = LennardJonesFluid(sigma=sigma, epsilon=epsilon)\n", "[topology, system, positions] = [fluid.topology, fluid.system, fluid.positions]\n", "\n", "# Add a barostat\n", "barostat = MonteCarloBarostat(pressure, temperature)\n", "system.addForce(barostat)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "695bedcf", "metadata": {}, "source": [ "To allow one of the Lennard-Jones particles to be alchemically eliminated, we create a [CustomNonbondedForce](http://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomNonbondedForce.html) that will compute the interactions between the alchemical particle and the remaining chemical particles using a softcore potential. The alchemically-modified particle has its Lennard-Jones well depth (epsilon parameter) set to zero in the original NonbondedForce, while the `CustomNonbondedForce` is set to evaluate only the interactions between the alchemically-modified particle and the remaining particles using [addInteractionGroup()](http://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomNonbondedForce.html#openmm.openmm.CustomNonbondedForce.addInteractionGroup) to specify only interactions between these groups are to be computed. A global context parameter lambda is created to control the coupling of the alchemically-modified particle with the rest of the system during the simulation. The Lennard-Jones parameters sigma and epsilon are implemented as per-particle parameters, though this is not strictly necessary in this case since all particles are equivalent." ] }, { "cell_type": "code", "execution_count": null, "id": "e5ebb1d7", "metadata": {}, "outputs": [], "source": [ "# Retrieve the NonbondedForce\n", "forces = { force.__class__.__name__ : force for force in system.getForces() }\n", "nbforce = forces['NonbondedForce']\n", "\n", "# Add a CustomNonbondedForce to handle only alchemically-modified interactions\n", "\n", "# Make two sets of particles, one that contains just the particle we will alchemically annihilate\n", "# and the other which contains all the other particles.\n", "alchemical_particles = set([0])\n", "chemical_particles = set(range(system.getNumParticles())) - alchemical_particles\n", "\n", "\n", "# Define the energy function for the CustomNonbondedForce\n", "# when lambda is 1.0 it is a normal LJ potential, when lambda is 0.0 the interaction vanishes \n", "energy_function = 'lambda*4*epsilon*x*(x-1.0); x = (sigma/reff_sterics)^6;'\n", "energy_function += 'reff_sterics = sigma*(0.5*(1.0-lambda) + (r/sigma)^6)^(1/6);'\n", "energy_function += 'sigma = 0.5*(sigma1+sigma2); epsilon = sqrt(epsilon1*epsilon2);'\n", "custom_force = CustomNonbondedForce(energy_function)\n", "\n", "# Add lambda as a parameter we can change during the simulation\n", "custom_force.addGlobalParameter('lambda', 1.0)\n", "\n", "# set the values of sigma and epsilon by copying them from the existing NonBondedForce\n", "custom_force.addPerParticleParameter('sigma')\n", "custom_force.addPerParticleParameter('epsilon')\n", "for index in range(system.getNumParticles()):\n", " [charge, sigma, epsilon] = nbforce.getParticleParameters(index)\n", " custom_force.addParticle([sigma, epsilon])\n", " if index in alchemical_particles:\n", " # remove the alchemical particle from the existing NonBondedForce\n", " nbforce.setParticleParameters(index, charge*0, sigma, epsilon*0)\n", "\n", "# Set the custom force to occur between just the alchemical particle and the other particles\n", "custom_force.addInteractionGroup(alchemical_particles, chemical_particles)\n", "system.addForce(custom_force)\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "e1945939", "metadata": {}, "source": [ "### Simulating alchemically-modified systems\n", "\n", "We then create a `LangevinIntegrator` and `Simulation` to run the simulation, and run a series of simulations at different values of lambda by using `simulation.context.setParameter()` to update the alchemical parameter on the fly. For each configuration sample that is collected, we can easily scan through the energy at different lambda values by simply alternating between `simulation.context.setParameter()` to update lambda and `simulation.context.getState()` to retrieve potential energies at the new alchemical state." ] }, { "cell_type": "code", "execution_count": null, "id": "6010430f", "metadata": {}, "outputs": [], "source": [ "# Create an integrator\n", "integrator = LangevinIntegrator(temperature, collision_rate, timestep)\n", "\n", "# Create a simulation\n", "simulation = Simulation(topology, system, integrator)\n", "simulation.context.setPositions(positions)\n", "\n", "# Minimize energy\n", "print('Minimizing energy...')\n", "#LocalEnergyMinimizer.minimize(context)\n", "simulation.minimizeEnergy()\n", "\n", "\n", "# Collect data\n", "\n", "# number of steps per sample\n", "nsteps = 2500\n", "\n", "# number of samples to collect per alchemical state\n", "niterations = 500\n", "\n", "import numpy as np\n", "lambdas = np.linspace(1.0, 0.0, 10) # alchemical lambda schedule\n", "nstates = len(lambdas)\n", "u_kln = np.zeros([nstates,nstates,niterations], np.float64)\n", "kT = AVOGADRO_CONSTANT_NA * BOLTZMANN_CONSTANT_kB * integrator.getTemperature()\n", "for k in range(nstates):\n", " for iteration in range(niterations):\n", " print('state %5d iteration %5d / %5d' % (k, iteration, niterations))\n", " # Set alchemical state\n", " simulation.context.setParameter('lambda', lambdas[k])\n", " # Run some dynamics\n", " simulation.step(nsteps)\n", " # Compute energies at all alchemical states\n", " for l in range(nstates):\n", " simulation.context.setParameter('lambda', lambdas[l])\n", " u_kln[k,l,iteration] = simulation.context.getState(getEnergy=True).getPotentialEnergy() / kT\n", "\n" ] }, { "cell_type": "markdown", "id": "05abedbb", "metadata": {}, "source": [ "### Analyzing the data with MBAR\n", "\n", "Finally, the [multistate Bennett acceptance ratio (MBAR)](https://dx.doi.org/10.1063%2F1.2978177) is used to estimate the free energy of annihilating the particle using the conda-installable [pymbar](http://pymbar.org/) package. In order to estimate how much data must be discarded to equilibration, we use a scheme for [automated equilibration detection](http://dx.doi.org/10.1021/acs.jctc.5b00784) and subsequent extraction of decorrelated samples found in the [pymbar.timeseries](http://github.com/choderalab/pymbar) module." ] }, { "cell_type": "code", "execution_count": null, "id": "e61ed834", "metadata": {}, "outputs": [], "source": [ "\n", "# Estimate free energy of Lennard-Jones particle insertion\n", "from pymbar import MBAR, timeseries\n", "\n", "## Subsample data to extract uncorrelated equilibrium timeseries\n", "N_k = np.zeros([nstates], np.int32) # number of uncorrelated samples\n", "for k in range(nstates):\n", " [nequil, g, Neff_max] = timeseries.detectEquilibration(u_kln[k,k,:])\n", " indices = timeseries.subsampleCorrelatedData(u_kln[k,k,:], g=g)\n", " N_k[k] = len(indices)\n", " u_kln[k,:,0:N_k[k]] = u_kln[k,:,indices].T\n", "\n", "# Compute free energy differences\n", "mbar = MBAR(u_kln, N_k)\n", "\n", "# dont compute uncertainties here, if you do it may fail with an error for\n", "# pymbar versions > 3.0.3. See this issue: https://github.com/choderalab/pymbar/issues/419\n", "[DeltaF_ij] = mbar.getFreeEnergyDifferences(compute_uncertainty=False)\n", "\n", "print(\"Free energy change to insert a particle = \",DeltaF_ij[nstates-1][0])" ] } ], "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.9" }, "nbsphinx": { "execute": "never" }, "tags": [ "tutorial" ], "vscode": { "interpreter": { "hash": "0196bde48eca6e98b9d405f49c02d8b38bda65759256048592a275ce0440cb5b" } } }, "nbformat": 4, "nbformat_minor": 5 }