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 *
config(viewer='webgl')

Please cite HTMD: Doerr et al.(2016)JCTC,12,1845.
https://dx.doi.org/10.1021/acs.jctc.6b00049
Documentation: http://software.acellera.com/
To update: conda update htmd -c acellera -c psi4

You are on the latest HTMD version (unpackaged : /home/joao/maindisk/software/repos/Acellera/htmd/htmd).

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)

2018-03-17 02:36:51,151 - htmd.molecule.readers - INFO - Using local copy for 3PTB: /home/joao/maindisk/software/repos/Acellera/htmd/htmd/data/pdb/3ptb.pdb
2018-03-17 02:36:51,421 - htmd.molecule.molecule - WARNING - Residue insertions were detected in the Molecule. It is recommended to renumber the residues using the Molecule.renumberResidues() method.
2018-03-17 02:36:51,457 - htmd.molecule.molecule - INFO - Removed 72 atoms. 1629 atoms remaining in the molecule.
2018-03-17 02:36:51,555 - propka - INFO - No pdbfile provided
2018-03-17 02:37:00,607 - htmd.builder.preparationdata - INFO - The following residues are in a non-standard state: CYS     6  A (CYX), HIS    22  A (HIE), CYS    24  A (CYX), HIS    39  A (HIP), CYS    40  A (CYX), GLU    51  A (GLH), HIS    72  A (HID), CYS   108  A (CYX), CYS   115  A (CYX), CYS   136  A (CYX), CYS   147  A (CYX), CYS   161  A (CYX), CYS   172  A (CYX), CYS   182  A (CYX), CYS   196  A (CYX), CYS   209  A (CYX)
2018-03-17 02:37:00,661 - htmd.builder.preparationdata - WARNING - Dubious protonation state: the pKa of 4 residues is within 1.0 units of pH 7.0.
2018-03-17 02:37:00,664 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    HIS    39  A (pKa= 7.46)
2018-03-17 02:37:00,666 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    GLU    51  A (pKa= 7.06)
2018-03-17 02:37:00,667 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    ASP   170  A (pKa= 6.49)
2018-03-17 02:37:00,669 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    N+      0T A (pKa= 7.49)
2018-03-17 02:37:00,748 - htmd.builder.preparationdata - WARNING - Found N-terminus 80.7% buried (> 50.0% threshold)
2018-03-17 02:37:00,749 - htmd.builder.preparationdata - WARNING - Found C-terminus involved in H bonds
2018-03-17 02:37:00,810 - htmd.builder.solvate - INFO - Using water pdb file at: /data/joao/maindisk/software/repos/Acellera/htmd/htmd/builder/wat.pdb
2018-03-17 02:37:02,115 - htmd.builder.solvate - INFO - Replicating 4 water segments, 2 by 1 by 2
Solvating: 100%|██████████| 4/4 [00:03<00:00,  1.10it/s]
2018-03-17 02:37:05,791 - 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')
tryp_charmm = charmm.build(tryp_solv, outdir='build-charmm')

2018-03-17 02:37:32,679 - htmd.builder.charmm - INFO - Writing out segments.
2018-03-17 02:37:34,330 - htmd.builder.builder - INFO - 6 disulfide bonds were added

Bond between A: [serial 351 resid 24 resname CYX chain A segid P0]
             B: [serial 571 resid 40 resname CYX chain A segid P0]

Bond between A: [serial 1684 resid 115 resname CYX chain A segid P0]
             B: [serial 2612 resid 182 resname CYX chain A segid P0]

Bond between A: [serial 2147 resid 147 resname CYX chain A segid P0]
             B: [serial 2354 resid 161 resname CYX chain A segid P0]

Bond between A: [serial 1605 resid 108 resname CYX chain A segid P0]
             B: [serial 3005 resid 209 resname CYX chain A segid P0]

Bond between A: [serial 2495 resid 172 resname CYX chain A segid P0]
             B: [serial 2800 resid 196 resname CYX chain A segid P0]

Bond between A: [serial 92 resid 6 resname CYX chain A segid P0]
             B: [serial 1989 resid 136 resname CYX chain A segid P0]

2018-03-17 02:37:34,835 - htmd.builder.charmm - INFO - Starting the build.
2018-03-17 02:37:35,225 - htmd.builder.charmm - WARNING - Failed to set coordinates for 9 atoms.
2018-03-17 02:37:35,227 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 3 atoms due to bad angles.
2018-03-17 02:37:35,228 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 16 atoms.
2018-03-17 02:37:35,229 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/build-charmm/log.txt for further information.
2018-03-17 02:37:35,230 - htmd.builder.charmm - INFO - Finished building.
2018-03-17 02:37:36,838 - htmd.builder.ionize - INFO - Adding 8 anions + 0 cations for neutralizing and 0 ions for the given salt concentration.
2018-03-17 02:37:37,042 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A
2018-03-17 02:37:37,043 - htmd.builder.ionize - INFO - Min distance between ions: 5A
2018-03-17 02:37:37,044 - htmd.builder.ionize - INFO - Placing 8 ions.
2018-03-17 02:37:39,577 - htmd.builder.charmm - INFO - Writing out segments.
2018-03-17 02:37:41,937 - htmd.builder.charmm - INFO - Starting the build.
2018-03-17 02:37:42,333 - htmd.builder.charmm - WARNING - Failed to set coordinates for 9 atoms.
2018-03-17 02:37:42,335 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 3 atoms due to bad angles.
2018-03-17 02:37:42,337 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 16 atoms.
2018-03-17 02:37:42,338 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/build-charmm/log.txt for further information.
2018-03-17 02:37:42,339 - htmd.builder.charmm - INFO - Finished building.

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-charmm/', './equil')

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

%ls equil/

input  parameters  [0m[38;5;34mrun.sh[0m*  structure.pdb  structure.psf

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

2018-03-17 02:39:07,892 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices
2018-03-17 02:39:07,953 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3
2018-03-17 02:39:08,112 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices
2018-03-17 02:39:08,175 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3
2018-03-17 02:39:08,179 - htmd.queues.localqueue - INFO - Queueing /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/equil
2018-03-17 02:39:08,181 - htmd.queues.localqueue - INFO - Running /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/equil on device 0
2018-03-17 02:39:33,879 - htmd.queues.localqueue - INFO - Completed /data/joao/maindisk/software/repos/Acellera/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.coor  input.xsc  parameters  [0m[38;5;34mrun.sh[0m*  structure.pdb  structure.psf

Run the production#

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

2018-03-17 02:39:56,043 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices
2018-03-17 02:39:56,103 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3
2018-03-17 02:39:56,440 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices
2018-03-17 02:39:56,501 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3
2018-03-17 02:39:56,505 - htmd.queues.localqueue - INFO - Queueing /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/prod
2018-03-17 02:39:56,507 - htmd.queues.localqueue - INFO - Running /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/prod on device 0
2018-03-17 02:44:05,268 - htmd.queues.localqueue - INFO - Completed /data/joao/maindisk/software/repos/Acellera/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()

A Jupyter Widget

Note: Waters are not for clarity