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-06-11 15:59:07,979 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
2024-06-11 15:59:07,980 - numexpr.utils - INFO - NumExpr defaulting to 8 threads.
2024-06-11 15:59:08,093 - rdkit - INFO - Enabling RDKit 2022.09.1 jupyter extensions
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.3.28+3.ga08b1aa18.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-06-11 15:59:10,521 - 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-06-11 15:59:11,807 - moleculekit.tools.preparation - INFO - Modified residue CYS 22 A to CYX
2024-06-11 15:59:11,808 - moleculekit.tools.preparation - INFO - Modified residue HIS 40 A to HIE
2024-06-11 15:59:11,809 - moleculekit.tools.preparation - INFO - Modified residue CYS 42 A to CYX
2024-06-11 15:59:11,809 - moleculekit.tools.preparation - INFO - Modified residue HIS 57 A to HIP
2024-06-11 15:59:11,810 - moleculekit.tools.preparation - INFO - Modified residue CYS 58 A to CYX
2024-06-11 15:59:11,811 - moleculekit.tools.preparation - INFO - Modified residue HIS 91 A to HID
2024-06-11 15:59:11,811 - moleculekit.tools.preparation - INFO - Modified residue CYS 128 A to CYX
2024-06-11 15:59:11,812 - moleculekit.tools.preparation - INFO - Modified residue CYS 136 A to CYX
2024-06-11 15:59:11,812 - moleculekit.tools.preparation - INFO - Modified residue CYS 157 A to CYX
2024-06-11 15:59:11,813 - moleculekit.tools.preparation - INFO - Modified residue CYS 168 A to CYX
2024-06-11 15:59:11,813 - moleculekit.tools.preparation - INFO - Modified residue CYS 182 A to CYX
2024-06-11 15:59:11,814 - moleculekit.tools.preparation - INFO - Modified residue CYS 191 A to CYX
2024-06-11 15:59:11,814 - moleculekit.tools.preparation - INFO - Modified residue CYS 201 A to CYX
2024-06-11 15:59:11,814 - moleculekit.tools.preparation - INFO - Modified residue CYS 220 A to CYX
2024-06-11 15:59:11,815 - moleculekit.tools.preparation - INFO - Modified residue CYS 232 A to CYX
2024-06-11 15:59:11,815 - moleculekit.tools.preparation - WARNING - Dubious protonation state: the pKa of 4 residues is within 1.0 units of pH 7.4.
2024-06-11 15:59:11,816 - moleculekit.tools.preparation - WARNING - Dubious protonation state: TYR 39 A (pKa= 8.24)
2024-06-11 15:59:11,816 - moleculekit.tools.preparation - WARNING - Dubious protonation state: HIS 57 A (pKa= 7.46)
2024-06-11 15:59:11,817 - moleculekit.tools.preparation - WARNING - Dubious protonation state: GLU 70 A (pKa= 7.06)
2024-06-11 15:59:11,817 - moleculekit.tools.preparation - WARNING - Dubious protonation state: ASP 189 A (pKa= 6.49)
2024-06-11 15:59:11,828 - moleculekit.tools.autosegment - INFO - Created segment P0 between resid 16 and 245.
2024-06-11 15:59:11,834 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-06-11 15:59:12,080 - htmd.builder.solvate - INFO - Replicating 4 water segments, 2 by 1 by 2
Solvating: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4/4 [00:01<00:00, 3.99it/s]
2024-06-11 15:59:13,127 - htmd.builder.solvate - INFO - 6943 water molecules were added to the system.
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-06-11 15:59:13,940 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-06-11 15:59:13,944 - 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-06-11 15:59:14,415 - htmd.builder.amber - INFO - Starting the build.
2024-06-11 15:59:15,264 - htmd.builder.amber - INFO - Finished building.
2024-06-11 15:59:17,575 - htmd.builder.ionize - INFO - Adding 7 anions + 0 cations for neutralizing and 0 ions for the given salt concentration 0 M.
2024-06-11 15:59:18,458 - htmd.builder.amber - INFO - Starting the build.
2024-06-11 15:59:19,344 - htmd.builder.amber - INFO - Finished building.
2024-06-11 15:59:21,460 - 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-06-11 15:59:22,474 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224
Equilibration protocol#
First, let’s import the Equilibration class, which is an MD
protocol:
from htmd.protocols.equilibration_v3 import Equilibration
The MD protocols are not imported automatically, the user must choose which one to use.
Now let’s start an Equilibration object, which already has sensible
defaults for an equilibration MD simulation and let’s define the
remaining ones:
md = Equilibration()
md.runtime = 1000
md.timeunits = 'fs'
md.temperature = 300
md.useconstantratio = False # only for membrane sims
# # Add a 10A flat bottom potential to prevent the ligand from diffusing from original position during equilibration
# width = np.array([10, 10, 10])
# flatbot = GroupRestraint('segname L and noh', width, [(5, '0ns')])
# md.restraints = [flatbot] + md.defaultEquilRestraints('20ns')
md.write('./build-amber/', './equil')
One can inspect what the Equilibration object has created with the
write method:
%ls equil/
input parameters [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-06-11 15:59:25,712 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-06-11 15:59:25,752 - jobqueues.localqueue - INFO - Using GPU devices 0
2024-06-11 15:59:25,753 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-06-11 15:59:25,783 - jobqueues.localqueue - INFO - Queueing /home/sdoerr/Work/htmd/tutorials/equil
2024-06-11 15:59:25,784 - jobqueues.localqueue - INFO - Running /home/sdoerr/Work/htmd/tutorials/equil on device 0
2024-06-11 15:59:42,979 - jobqueues.localqueue - INFO - Completed /home/sdoerr/Work/htmd/tutorials/equil
Production protocol#
Now let’s use the Production class and do the same as before to
perform a short production MD simulation:
from htmd.protocols.production_v6 import Production
md = Production()
md.runtime = 1
md.timeunits = 'ns'
md.temperature = 300
md.acemd.bincoordinates = 'output.coor'
md.acemd.extendedsystem = 'output.xsc'
md.write('equil','prod')
%ls prod/
input input.xsc [0m[01;32mrun.sh[0m* structure.prmtop
input.coor parameters structure.pdb
Run the production#
local = LocalGPUQueue()
local.submit('./prod/')
local.wait()
2024-06-11 16:00:02,235 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-06-11 16:00:02,280 - jobqueues.localqueue - INFO - Using GPU devices 0
2024-06-11 16:00:02,281 - jobqueues.util - INFO - Trying to determine all GPU devices
2024-06-11 16:00:02,309 - jobqueues.localqueue - INFO - Queueing /home/sdoerr/Work/htmd/tutorials/prod
2024-06-11 16:00:02,310 - jobqueues.localqueue - INFO - Running /home/sdoerr/Work/htmd/tutorials/prod on device 0
2024-06-11 16:02:42,625 - 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.pdb')
traj.read('prod/output.xtc')
traj.wrap()
traj.align('protein and name CA')
traj.view()
2024-06-11 16:02:47,993 - moleculekit.molecule - WARNING - Wrapping detected 0 bonds and 24045 atoms. Ignore this message if you believe this is accurate, otherwise make sure you have loaded a topology containing all the bonds of the system before wrapping. The results may be inaccurate. If you want to use guessed bonds use the guessBonds argument.
NGLWidget(max_frame=9)
Note: Waters are not shown for clarity