Molecular dynamics protocols in HTMD#

In this tutorial, we should how to equilibrate and prepare a system for productive simulations in HTMD.

First, let’s import HTMD:

from htmd.ui import *
from moleculekit.config import config

config(viewer='webgl')

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.
2024-12-05 13:27:36,940 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.
/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).
  from pandas.core.computation.check import NUMEXPR_INSTALLED

Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049
HTMD Documentation at: https://software.acellera.com/htmd/

You are on the latest HTMD version (2.4.2+20.g74f7135a7.dirty).

Build a sample system#

We need a built system to perform simulations. We can easily build one with:

tryp = Molecule("3PTB")
tryp.filter("protein")
tryp_op = systemPrepare(tryp)
tryp_seg = autoSegment(tryp_op)
tryp_solv = solvate(tryp_seg, pad=10)

2024-12-05 13:27:39,407 - moleculekit.molecule - INFO - Removed 72 atoms. 1629 atoms remaining in the molecule.

---- Molecule chain report ----
Chain A:
    First residue: ILE    16
    Final residue: ASN   245
---- End of chain report ----

2024-12-05 13:27:40,621 - moleculekit.tools.preparation - INFO - Modified residue CYS    22 A to CYX
2024-12-05 13:27:40,621 - moleculekit.tools.preparation - INFO - Modified residue HIS    40 A to HIE
2024-12-05 13:27:40,621 - moleculekit.tools.preparation - INFO - Modified residue CYS    42 A to CYX
2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue HIS    57 A to HIP
2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue CYS    58 A to CYX
2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue HIS    91 A to HID
2024-12-05 13:27:40,622 - moleculekit.tools.preparation - INFO - Modified residue CYS   128 A to CYX
2024-12-05 13:27:40,623 - moleculekit.tools.preparation - INFO - Modified residue CYS   136 A to CYX
2024-12-05 13:27:40,623 - moleculekit.tools.preparation - INFO - Modified residue CYS   157 A to CYX
2024-12-05 13:27:40,624 - moleculekit.tools.preparation - INFO - Modified residue CYS   168 A to CYX
2024-12-05 13:27:40,625 - moleculekit.tools.preparation - INFO - Modified residue CYS   182 A to CYX
2024-12-05 13:27:40,625 - moleculekit.tools.preparation - INFO - Modified residue CYS   191 A to CYX
2024-12-05 13:27:40,625 - moleculekit.tools.preparation - INFO - Modified residue CYS   201 A to CYX
2024-12-05 13:27:40,626 - moleculekit.tools.preparation - INFO - Modified residue CYS   220 A to CYX
2024-12-05 13:27:40,627 - moleculekit.tools.preparation - INFO - Modified residue CYS   232 A to CYX
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.
2024-12-05 13:27:40,628 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    TYR    39 A (pKa= 8.24)
2024-12-05 13:27:40,629 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    HIS    57 A (pKa= 7.46)
2024-12-05 13:27:40,629 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    GLU    70 A (pKa= 7.06)
2024-12-05 13:27:40,630 - moleculekit.tools.preparation - WARNING - Dubious protonation state:    ASP   189 A (pKa= 6.49)
2024-12-05 13:27:40,640 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 16 and 245.
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
2024-12-05 13:27:40,912 - htmd.builder.solvate - INFO - Replicating 4 water segments, 2 by 1 by 2

2024-12-05 13:27:41,299 - htmd.builder.solvate - INFO - 6943 water molecules were added to the system.
██████████████████████████████████████████████████████████████████████████| 4/4 [00:00<00:00, 11.70it/s]

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):

tryp_amber = amber.build(tryp_solv, outdir='build-amber')

2024-12-05 13:27:42,181 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-12-05 13:27:42,195 - htmd.builder.builder - INFO - 6 disulfide bonds were added

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 26, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 42, insertion: '', segid: 'P0'>

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 117, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 184, insertion: '', segid: 'P0'>

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 174, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 198, insertion: '', segid: 'P0'>

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 149, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 163, insertion: '', segid: 'P0'>

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 8, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 138, insertion: '', segid: 'P0'>

Disulfide Bond between: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 110, insertion: '', segid: 'P0'>
                   and: UniqueResidueID<resname: 'CYX', chain: 'A', resid: 211, insertion: '', segid: 'P0'>

2024-12-05 13:27:42,717 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 13:27:43,634 - htmd.builder.amber - INFO - Finished building.
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.
2024-12-05 13:27:46,983 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 13:27:47,947 - htmd.builder.amber - INFO - Finished building.
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.
  warnings.warn(

2024-12-05 13:27:51,311 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224

Equilibration protocol#

First, let’s import the setup_equilibration function, which is an MD protocol:

from acemd.protocols import setup_equilibration

2024-12-05 13:27:51,810 - rdkit - INFO - Enabling RDKit 2024.03.5 jupyter extensions
2024-12-05 13:27:51,894 - acemd - INFO - # You are on the latest ACEMD version (4.0.1).

The setup_equilibration function has sensible defaults for an equilibration MD simulation but it also allows the user to override any options:

setup_equilibration('./build-amber/', './equil', run="1ns", temperature=300)

2024-12-05 13:27:51,898 - acemd - INFO - Copying ./build-amber/structure.pdb to ./equil/structure.pdb
2024-12-05 13:27:51,900 - acemd - INFO - Copying ./build-amber/structure.prmtop to ./equil/structure.prmtop

One can inspect what the Equilibration object has created with the write method:

%ls equil/

input.yaml  [0m[01;32mrun.sh[0m*  structure.pdb  structure.prmtop

Run the equilibration#

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:

local = LocalGPUQueue()
local.submit('./equil/')
local.wait()

2024-12-05 13:28:02,499 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-12-05 13:28:02,548 - jobqueues.localqueue - INFO - Using GPU devices 0
2024-12-05 13:28:02,549 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-12-05 13:28:02,594 - jobqueues.localqueue - INFO - Queueing /home/sdoerr/Work/htmd/tutorials/equil
2024-12-05 13:28:02,596 - jobqueues.localqueue - INFO - Running /home/sdoerr/Work/htmd/tutorials/equil on device 0
2024-12-05 13:31:44,008 - jobqueues.localqueue - INFO - Completed /home/sdoerr/Work/htmd/tutorials/equil

Let’s check which files have been created by the simulation

%ls equil/

input.yaml      minimized.coor  output.xsc       restraint_1.sel
jobqueues.done  output.coor     output.xtc       [0m[01;32mrun.sh[0m*
[01;32mjob.sh[0m*         output.csv      restart.chk      structure.pdb
log.txt         output.vel      restraint_0.sel  structure.prmtop

Production protocol#

Now let’s use the setup_production function and do the same as before to perform a short production MD simulation:

from acemd.protocols import setup_production

setup_production('equil', 'prod', run="1ns", temperature=300)

2024-12-05 13:32:02,746 - acemd - INFO - Copying equil/output.coor to prod/input.coor
2024-12-05 13:32:02,747 - acemd - INFO - Copying equil/output.xsc to prod/input.xsc
2024-12-05 13:32:02,748 - acemd - INFO - Copying equil/structure.prmtop to prod/structure.prmtop

%ls prod/

input.coor  input.xsc  input.yaml  [0m[01;32mrun.sh[0m*  structure.prmtop

Run the production#

local = LocalGPUQueue()
local.submit('./prod/')
local.wait()

2024-12-05 13:32:26,221 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-12-05 13:32:26,275 - jobqueues.localqueue - INFO - Using GPU devices 0
2024-12-05 13:32:26,276 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-12-05 13:32:26,321 - jobqueues.localqueue - INFO - Queueing /home/sdoerr/Work/htmd/tutorials/prod
2024-12-05 13:32:26,323 - jobqueues.localqueue - INFO - Running /home/sdoerr/Work/htmd/tutorials/prod on device 0
2024-12-05 13:35:27,738 - jobqueues.localqueue - INFO - Completed /home/sdoerr/Work/htmd/tutorials/prod

When the above process finished, one can check that a trajectory was produced:

%ls prod/output.xtc

prod/output.xtc

Quickly visualize the trajectory#

traj = Molecule('prod/structure.prmtop')
traj.read('prod/output.xtc')
traj.wrap()
traj.align('protein and name CA')
traj.view()

NGLWidget(max_frame=9)

Note: Waters are not shown for clarity