{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Molecular dynamics protocols in HTMD" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "In this tutorial, we should how to equilibrate and prepare a system for productive simulations in HTMD.\n", "\n", "First, let's import HTMD:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:36,940 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but \"NUMEXPR_MAX_THREADS\" not set, so enforcing safe limit of 8.\n", "2024-12-05 13:27:36,940 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.\n", "/home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/pandas/core/computation/expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).\n", " from pandas.core.computation.check import NUMEXPR_INSTALLED\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049\n", "HTMD Documentation at: https://software.acellera.com/htmd/\n", "\n", "You are on the latest HTMD version (2.4.2+20.g74f7135a7.dirty).\n", "\n" ] } ], "source": [ "from htmd.ui import *\n", "from moleculekit.config import config\n", "\n", "config(viewer='webgl')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Build a sample system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need a built system to perform simulations. We can easily build one with:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:39,407 - moleculekit.molecule - INFO - Removed 72 atoms. 1629 atoms remaining in the molecule.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "---- Molecule chain report ----\n", "Chain A:\n", " First residue: ILE 16 \n", " Final residue: ASN 245 \n", "---- End of chain report ----\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:40,621 - moleculekit.tools.preparation - INFO - Modified residue CYS 22 A to CYX\n", "2024-12-05 13:27:40,621 - moleculekit.tools.preparation - INFO - Modified residue HIS 40 A to HIE\n", "2024-12-05 13:27:40,621 - moleculekit.tools.preparation - INFO - Modified residue CYS 42 A to CYX\n", "2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue HIS 57 A to HIP\n", "2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue CYS 58 A to CYX\n", "2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue HIS 91 A to HID\n", "2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue CYS 128 A to CYX\n", "2024-12-05 13:27:40,623 - moleculekit.tools.preparation - INFO - Modified residue CYS 136 A to CYX\n", "2024-12-05 13:27:40,623 - moleculekit.tools.preparation - INFO - Modified residue CYS 157 A to CYX\n", "2024-12-05 13:27:40,624 - moleculekit.tools.preparation - INFO - Modified residue CYS 168 A to CYX\n", "2024-12-05 13:27:40,625 - moleculekit.tools.preparation - INFO - Modified residue CYS 182 A to CYX\n", "2024-12-05 13:27:40,625 - moleculekit.tools.preparation - INFO - Modified residue CYS 191 A to CYX\n", "2024-12-05 13:27:40,625 - moleculekit.tools.preparation - INFO - Modified residue CYS 201 A to CYX\n", "2024-12-05 13:27:40,626 - moleculekit.tools.preparation - INFO - Modified residue CYS 220 A to CYX\n", "2024-12-05 13:27:40,627 - moleculekit.tools.preparation - INFO - Modified residue CYS 232 A to CYX\n", "2024-12-05 13:27:40,628 - moleculekit.tools.preparation - WARNING - Dubious protonation state: the pKa of 4 residues is within 1.0 units of pH 7.4.\n", "2024-12-05 13:27:40,628 - moleculekit.tools.preparation - WARNING - Dubious protonation state: TYR 39 A (pKa= 8.24)\n", "2024-12-05 13:27:40,629 - moleculekit.tools.preparation - WARNING - Dubious protonation state: HIS 57 A (pKa= 7.46)\n", "2024-12-05 13:27:40,629 - moleculekit.tools.preparation - WARNING - Dubious protonation state: GLU 70 A (pKa= 7.06)\n", "2024-12-05 13:27:40,630 - moleculekit.tools.preparation - WARNING - Dubious protonation state: ASP 189 A (pKa= 6.49)\n", "2024-12-05 13:27:40,640 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 16 and 245.\n", "2024-12-05 13:27:40,646 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb\n", "2024-12-05 13:27:40,912 - htmd.builder.solvate - INFO - Replicating 4 water segments, 2 by 1 by 2\n", "\n", "2024-12-05 13:27:41,299 - htmd.builder.solvate - INFO - 6943 water molecules were added to the system.\n", "██████████████████████████████████████████████████████████████████████████| 4/4 [00:00<00:00, 11.70it/s]" ] } ], "source": [ "tryp = Molecule(\"3PTB\")\n", "tryp.filter(\"protein\")\n", "tryp_op = systemPrepare(tryp)\n", "tryp_seg = autoSegment(tryp_op)\n", "tryp_solv = solvate(tryp_seg, pad=10)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "With this solvated system and since HTMD is force-field agnostic, one can either build in CHARMM or Amber (for details on system building, see the corresponding tutorials):" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:42,181 - htmd.builder.amber - INFO - Detecting disulfide bonds.\n", "2024-12-05 13:27:42,195 - htmd.builder.builder - INFO - 6 disulfide bonds were added\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Disulfide Bond between: UniqueResidueID\n", " and: UniqueResidueID\n", "\n", "Disulfide Bond between: UniqueResidueID\n", " and: UniqueResidueID\n", "\n", "Disulfide Bond between: UniqueResidueID\n", " and: UniqueResidueID\n", "\n", "Disulfide Bond between: UniqueResidueID\n", " and: UniqueResidueID\n", "\n", "Disulfide Bond between: UniqueResidueID\n", " and: UniqueResidueID\n", "\n", "Disulfide Bond between: UniqueResidueID\n", " and: UniqueResidueID\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:42,717 - htmd.builder.amber - INFO - Starting the build.\n", "2024-12-05 13:27:43,634 - htmd.builder.amber - INFO - Finished building.\n", "2024-12-05 13:27:46,066 - htmd.builder.ionize - INFO - Adding 7 anions + 0 cations for neutralizing and 0 ions for the given salt concentration 0 M.\n", "2024-12-05 13:27:46,983 - htmd.builder.amber - INFO - Starting the build.\n", "2024-12-05 13:27:47,947 - htmd.builder.amber - INFO - Finished building.\n", "2024-12-05 13:27:50,278 - py.warnings - WARNING - /home/sdoerr/miniforge3/envs/htmd/lib/python3.10/site-packages/Bio/pairwise2.py:278: BiopythonDeprecationWarning: Bio.pairwise2 has been deprecated, and we intend to remove it in a future release of Biopython. As an alternative, please consider using Bio.Align.PairwiseAligner as a replacement, and contact the Biopython developers if you still need the Bio.pairwise2 module.\n", " warnings.warn(\n", "\n", "2024-12-05 13:27:51,311 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224\n" ] } ], "source": [ "tryp_amber = amber.build(tryp_solv, outdir='build-amber')" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Equilibration protocol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's import the `setup_equilibration` function, which is an MD protocol:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:51,810 - rdkit - INFO - Enabling RDKit 2024.03.5 jupyter extensions\n", "2024-12-05 13:27:51,894 - acemd - INFO - # You are on the latest ACEMD version (4.0.1).\n" ] } ], "source": [ "from acemd.protocols import setup_equilibration" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The `setup_equilibration` function has sensible defaults for an equilibration MD simulation but it also allows the user to override any options:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:27:51,898 - acemd - INFO - Copying ./build-amber/structure.pdb to ./equil/structure.pdb\n", "2024-12-05 13:27:51,900 - acemd - INFO - Copying ./build-amber/structure.prmtop to ./equil/structure.prmtop\n" ] } ], "source": [ "setup_equilibration('./build-amber/', './equil', run=\"1ns\", temperature=300)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One can inspect what the `Equilibration` object has created with the `write` method:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "input.yaml \u001b[0m\u001b[01;32mrun.sh\u001b[0m* structure.pdb structure.prmtop\n" ] } ], "source": [ "%ls equil/" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Run the equilibration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's use the queue resources of HTMD to run the simulation. One can use the local computer and the `LocalGPUQueue` class to submit a job. The equilibration time is short (see above) for demonstration purposes:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:28:02,499 - jobqueues.util - INFO - Trying to determine all GPU devices\n", "2024-12-05 13:28:02,548 - jobqueues.localqueue - INFO - Using GPU devices 0\n", "2024-12-05 13:28:02,549 - jobqueues.util - INFO - Trying to determine all GPU devices\n", "2024-12-05 13:28:02,594 - jobqueues.localqueue - INFO - Queueing /home/sdoerr/Work/htmd/tutorials/equil\n", "2024-12-05 13:28:02,596 - jobqueues.localqueue - INFO - Running /home/sdoerr/Work/htmd/tutorials/equil on device 0\n", "2024-12-05 13:31:44,008 - jobqueues.localqueue - INFO - Completed /home/sdoerr/Work/htmd/tutorials/equil\n" ] } ], "source": [ "local = LocalGPUQueue()\n", "local.submit('./equil/')\n", "local.wait()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check which files have been created by the simulation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "input.yaml minimized.coor output.xsc restraint_1.sel\n", "jobqueues.done output.coor output.xtc \u001b[0m\u001b[01;32mrun.sh\u001b[0m*\n", "\u001b[01;32mjob.sh\u001b[0m* output.csv restart.chk structure.pdb\n", "log.txt output.vel restraint_0.sel structure.prmtop\n" ] } ], "source": [ "%ls equil/" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Production protocol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's use the `setup_production` function and do the same as before to perform a short production MD simulation:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "from acemd.protocols import setup_production" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:32:02,746 - acemd - INFO - Copying equil/output.coor to prod/input.coor\n", "2024-12-05 13:32:02,747 - acemd - INFO - Copying equil/output.xsc to prod/input.xsc\n", "2024-12-05 13:32:02,748 - acemd - INFO - Copying equil/structure.prmtop to prod/structure.prmtop\n" ] } ], "source": [ "setup_production('equil', 'prod', run=\"1ns\", temperature=300)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "input.coor input.xsc input.yaml \u001b[0m\u001b[01;32mrun.sh\u001b[0m* structure.prmtop\n" ] } ], "source": [ "%ls prod/" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Run the production" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2024-12-05 13:32:26,221 - jobqueues.util - INFO - Trying to determine all GPU devices\n", "2024-12-05 13:32:26,275 - jobqueues.localqueue - INFO - Using GPU devices 0\n", "2024-12-05 13:32:26,276 - jobqueues.util - INFO - Trying to determine all GPU devices\n", "2024-12-05 13:32:26,321 - jobqueues.localqueue - INFO - Queueing /home/sdoerr/Work/htmd/tutorials/prod\n", "2024-12-05 13:32:26,323 - jobqueues.localqueue - INFO - Running /home/sdoerr/Work/htmd/tutorials/prod on device 0\n", "2024-12-05 13:35:27,738 - jobqueues.localqueue - INFO - Completed /home/sdoerr/Work/htmd/tutorials/prod\n" ] } ], "source": [ "local = LocalGPUQueue()\n", "local.submit('./prod/')\n", "local.wait()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When the above process finished, one can check that a trajectory was produced:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "prod/output.xtc\n" ] } ], "source": [ "%ls prod/output.xtc" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Quickly visualize the trajectory" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a810a0f07599487a880a6816846bbd58", "version_major": 2, "version_minor": 0 }, "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "06b47eab2def4fcbb9a34468d5585b1e", "version_major": 2, "version_minor": 0 }, "text/plain": [ "NGLWidget(max_frame=9)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "traj = Molecule('prod/structure.prmtop')\n", "traj.read('prod/output.xtc')\n", "traj.wrap()\n", "traj.align('protein and name CA')\n", "traj.view()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** Waters are not shown for clarity" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Slideshow", "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.15" } }, "nbformat": 4, "nbformat_minor": 4 }