{ "cells": [ { "cell_type": "markdown", "id": "28f73de5-b88f-465f-ac6a-768f7ad70dba", "metadata": {}, "source": [ "# Compute the alanine dipeptide [$\\Phi,\\Psi$] transition using OpenMM and OPES\n", "\n", "In this brief tutorial, we calculate the Ramachandran [$\\Phi$,$\\Psi$] plot of alanine dipeptide in vacuum using the OPES method." ] }, { "cell_type": "markdown", "id": "e6e82f00-549d-4daf-aceb-0d645de7567e", "metadata": {}, "source": [ "Alanine dipeptide is a popular test system for enhanced sampling algorithms because it is a minimal example of sampling challenges that are posed by many biological systems. The slow motions of the molecule are largely governed by the $\\Psi$ (backbone N-C-C-N) and $\\Phi$ (backbone C-N-C-C) dihedrals." ] }, { "cell_type": "markdown", "id": "ae8fbc05-6630-4a10-83cb-99252cc75f51", "metadata": {}, "source": [ "An interactive visualization of the alanine dipeptide molecule can be shown with nglview. Atoms involved in the [$\\Phi$,$\\Psi$] collective variable (CV) are shown in ball and stick representation, while other atoms are transparent. " ] }, { "cell_type": "code", "execution_count": 2, "id": "6892a570-7bcc-48c9-9ec9-f0e47399d186", "metadata": {}, "outputs": [], "source": [ "import pytraj as pt\n", "import nglview as ngl\n", "\n", "traj = pt.load(\"../data/alanine-dipeptide.pdb\", top='../data/alanine-dipeptide.prmtop')\n", "\n", "view = ngl.show_pytraj(traj)\n", "view.clear_representations()\n", "view.add_ball_and_stick('@6,8,14,16', opacity=1.0) # CV atoms of phi torsion\n", "view.add_ball_and_stick('@4,6,8,14', opacity=1.0) # CV atoms of psi torsion\n", "#view.add_ball_and_stick('@1,4,6,8', opacity=1.0) # CV atoms of theta torsion\n", "\n", "view.add_licorice(opacity=0.5)\n", "\n", "# uncomment for interactive visualization in a notebook\n", "#view" ] }, { "cell_type": "markdown", "id": "23a07d94-39c8-447a-954d-71fec9326cb3", "metadata": {}, "source": [ "While the interactive visualization is not possible on the web page, this is what you should see in a notebook:" ] }, { "cell_type": "markdown", "id": "7735adcf-997e-498d-a3e1-7c765994e9f3", "metadata": {}, "source": [ "![](nglview_snapshot.png)" ] }, { "cell_type": "markdown", "id": "2415f7cc-c021-4442-ae81-4c2187d2eb60", "metadata": {}, "source": [ "## Importance sampling of the [$\\Phi,\\Psi$] plane with OPES\n", "\n", "Below, the OPES sampling algorithm is applied to enhance sampling in the [$\\Phi,\\Psi$] plane. \n", "For this purpose, an OpenMM simulation object is created using AMBER parameter and coordinate files and the OpenMM interface to the adaptive sampling package. " ] }, { "cell_type": "code", "execution_count": 3, "id": "b2ad231f-9d4d-401b-a944-0b454216f044", "metadata": { "scrolled": true }, "outputs": [], "source": [ "from sys import stdout\n", "import numpy as np\n", "\n", "from openmm import *\n", "from openmm.app import *\n", "from openmm.unit import *\n", "\n", "from adaptive_sampling.interface.interface_openmm import AdaptiveSamplingOpenMM\n", "\n", "# ------------------------------------------------------------------------------------\n", "# define collective variables\n", "cv_atoms_psi = [6, 8, 14, 16] # backbone N-C-C-N torsion\n", "cv_atoms_phi = [4, 6, 8, 14] # backbone C-N-C-C torsion\n", "minimum = -180.0 # minimum of the CV\n", "maximum = 180.0 # maximum of the CV\n", "bin_width = 5.0 # bin width along the CV\n", "periodicity = [ # define periodicity of CVs (needs to be given in radians as units are not converted)\n", " [-np.pi, np.pi],\n", " [-np.pi, np.pi],\n", "]\n", "\n", "collective_var = [\n", " [\"torsion\", cv_atoms_psi, minimum, maximum, bin_width],\n", " [\"torsion\", cv_atoms_phi, minimum, maximum, bin_width],\n", "]\n", "\n", "# ------------------------------------------------------------------------------------\n", "# Setup OpenMM\n", "prmtop = AmberPrmtopFile(f\"../data/alanine-dipeptide.prmtop\")\n", "crd = AmberInpcrdFile(f\"../data/alanine-dipeptide.crd\")\n", "system = prmtop.createSystem(\n", " nonbondedMethod=NoCutoff,\n", " constraints=HBonds,\n", ")\n", "\n", "# Initialize the `AdaptiveSamplingOpenMM` interface to couple the OpenMM simulation to a bias potential\n", "# The OpenMM `simulation` object is set up internally, but can still be modified by calling `the_md.simulation` or `the_md.integrator`\n", "the_md = AdaptiveSamplingOpenMM(\n", " crd.positions,\n", " prmtop.topology,\n", " system,\n", " dt=2.0, # timestep in fs\n", " equil_temp=300.0, # temperature of simulation\n", " langevin_damping=1.0, # langevin damping in 1/ps\n", " cv_atoms=np.unique(cv_atoms_phi+cv_atoms_psi) # specifying CV atoms significantly speeds up simulation of large systems, as the bias force will only be calculated for those\n", ")\n", "the_md.integrator.setConstraintTolerance(0.00001)\n", "\n", "# Append OpenMM reporters to simulation for output \n", "the_md.simulation.reporters.append(DCDReporter('alanine-dipeptide.dcd', 100))\n", "the_md.simulation.reporters.append(StateDataReporter(\n", " stdout, \n", " 100,\n", " step=True,\n", " time=True,\n", " potentialEnergy=True,\n", " kineticEnergy=True,\n", " totalEnergy=True,\n", " temperature=True,\n", " speed=False,\n", " separator='\\t')\n", ")" ] }, { "cell_type": "markdown", "id": "5d896eed-2476-4f5b-8b4e-d5e0e30d5a5e", "metadata": {}, "source": [ "OpenMM then is provided with an instance of the `OPES` class as the bias. \n", "Although the OPES implementation let the user toggle between different variants of various features, the defaults for the corresponding settings have turned out to be sufficient for most cases. This minimizes the user input to a small set of parameters. Below is an example of an OPES instance that relies on this implementation's defaults for efficient computation but provides the possibility to take a deeper look into OPES' behaviour for the parameters, that can directly affect the simulations quality.\n", "- Compare the results when providing OPES with a small initial `opes_kernel_std` in comparison to automatically estimating it to fit the starting basin.\n", "- Compare the effects of different barriere heights `opes_barrier` and how limiting the OPES bias to smaller barriers than the system's one." ] }, { "cell_type": "code", "execution_count": 4, "id": "a5a9fcef-7d41-42fb-be06-a66feb2c96da", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", " >>> INFO: Initialize torsion as collective variable:\n", "\t Minimum0:\t-180.0 Degree\n", "\t Maximum0:\t180.0 Degree\n", "\t Bin width0:\t5.0 Degree\n", "\n", " >>> INFO: Initialize torsion as collective variable:\n", "\t Minimum1:\t-180.0 Degree\n", "\t Maximum1:\t180.0 Degree\n", "\t Bin width1:\t5.0 Degree\n", "\t----------------------------------------------\n", "\t Total number of bins:\t\t5184\n", "\n", " >>> INFO: OPES Parameters:\n", "\t ---------------------------------------------\n", "\t Kernel_std:\testimate from 5000 steps\n", "\t Rescaling:\tTrue\n", "\t Adaptive:\tFalse\t(5000 steps)\n", "\t Normalize:\tTrue\t(approximated: True)\n", "\t Explore:\tFalse\n", "\t Barrier:\t11.9503 kcal/mol\n", "\t Bias factor:\t20.045403709021212\n", "\t Read force:\tTrue\n", "\t Kernel merge:\tTrue\t(threshold: 1.0)\n", "\t ---------------------------------------------\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ">>> adaptive-sampling: Module 'ase' not found, will not import 'FENEB'\n", "/home/robert/Bachelor_Thesis/Code/adaptive_sampling/adaptive_sampling/colvars/colvars.py:215: UserWarning: Using torch.cross without specifying the dim arg is deprecated.\n", "Please either pass the dim explicitly or simply use torch.linalg.cross.\n", "The default value of dim will change to agree with that of linalg.cross in a future release. (Triggered internally at /opt/conda/conda-bld/pytorch_1720538438750/work/aten/src/ATen/native/Cross.cpp:62.)\n", " self.cv = torch.atan2(torch.dot(torch.cross(q23_u, n1), n2), torch.dot(n1, n2))\n" ] } ], "source": [ "from adaptive_sampling.sampling_tools.opes import OPES\n", "# --------------------------------------------------------------------------------------\n", "# Setup the sampling algorithm\n", "\n", "opes_kernel_std = None # kernel standard deviation\n", "opes_frequency = 500 # frequency of kernel creation in MD steps\n", "opes_barrier = 50.0 # Barrier parameter in kJ/mol \n", "\n", "the_bias = OPES(\n", " the_md, \n", " collective_var, # collective variable\n", " # OPES parameters\n", " kernel_std = opes_kernel_std,\n", " update_freq = opes_frequency,\n", " energy_barr = opes_barrier,\n", " \n", " # general parameters\n", " output_freq = 1000, # frequency of writing outputs\n", " f_conf = 0.0, # confinement force of CV at boundaries\n", " equil_temp = 300.0, # equilibrium temperature of simulation\n", " periodicity = periodicity, # periodicity of CVs\n", " verbose = True, # print verbose output\n", ")\n", "the_md.set_sampling_algorithm(the_bias) # to take effect the sampling algorithm has to be set in the MD interface" ] }, { "cell_type": "code", "execution_count": 8, "id": "807819d2-caf4-4608-8a9f-c77d789015d1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1100\t2.1999999999999793\t-77.484027451166\t70.72776178943762\t-6.756265661728378\t333.5919402869253\n", "1200\t2.3999999999999573\t-63.594587427905225\t51.409461789394186\t-12.185125638511039\t242.4759623878226\n", "1300\t2.5999999999999353\t-70.42487310290143\t53.21313727800854\t-17.21173582489289\t250.9831113583515\n", "1400\t2.7999999999999132\t-52.20514107819585\t50.758430490201945\t-1.4467105879939055\t239.40533228741566\n", "1500\t2.999999999999891\t-52.521335364644656\t55.41707830529658\t2.895742940651921\t261.3781379358894\n", "1600\t3.199999999999869\t-49.08489702289795\t62.4755462465003\t13.390649223602352\t294.66984626067983\n", "1700\t3.399999999999847\t-46.31697749196262\t59.84464743883469\t13.527669946872074\t282.26104643805064\n", "1800\t3.599999999999825\t-55.343389270210565\t71.54269934628653\t16.19931007607596\t337.43564456829716\n", "1900\t3.799999999999803\t-59.52420207082332\t75.13378924961816\t15.60958717879484\t354.37324613080114\n", "2000\t3.999999999999781\t-60.86357611569707\t48.45052394606737\t-12.413052169629701\t228.51994580578491\n" ] } ], "source": [ "# WARNING: For long simulations, this can become expensive, and it is recommended to perform the computation on an HPC cluster.\n", "if os.path.isfile('CV_traj.dat'):\n", " os.system(\"rm CV_traj.dat opeseabf.out\")\n", "the_md.run(nsteps=1000) # 500000 * 2 fs = 1 ns" ] }, { "cell_type": "markdown", "id": "f87064b7-2090-4a6a-a6b6-2b63370fc60a", "metadata": {}, "source": [ "## Analysis of Results\n", "\n", "### Visualize the trajectory with NGlView\n", "\n", "The following will create an interactive visualization of the trajectory:" ] }, { "cell_type": "code", "execution_count": 10, "id": "fd888874-3b8e-4ea6-8d6d-aa33a22a0805", "metadata": {}, "outputs": [], "source": [ "path = '.' \n", "traj = pt.iterload(f\"{path}/alanine-dipeptide.dcd\", top=\"../data/alanine-dipeptide.pdb\")\n", "\n", "view = ngl.show_pytraj(traj)\n", "view.clear_representations()\n", "view.add_ball_and_stick('@6,8,14,16', opacity=1.0) # CV atoms of phi torsion\n", "view.add_ball_and_stick('@4,6,8,14', opacity=1.0) # CV atoms of psi torsion\n", "#view.add_ball_and_stick('@1,4,6,8', opacity=1.0) # CV atoms of theta torsion\n", "\n", "view.add_licorice(opacity=0.5)\n", "#view" ] }, { "cell_type": "markdown", "id": "286ab551-b83a-44d0-980c-97a1adef21d5", "metadata": {}, "source": [ "### Trajectory of CVs\n", "\n", "All important information for post/processing of the trajectory is stored in `CV_traj.dat`." ] }, { "cell_type": "code", "execution_count": 11, "id": "0b24cd68-2155-4d1f-8bf3-2cca030ce8a1", "metadata": {}, "outputs": [], "source": [ "cv_traj = np.loadtxt(f'{path}/CV_traj.dat', skiprows=1)\n", "cv_phi = cv_traj[:,1] # Phi trajectory (first CV)\n", "cv_psi = cv_traj[:,2] # Psi trajectory (second CV))" ] }, { "cell_type": "code", "execution_count": null, "id": "79e1f61a", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "hiwi", "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.20" } }, "nbformat": 4, "nbformat_minor": 5 }