from htmd.ui import *
from moleculekit.config import config
from moleculekit.util import maxDistance
from htmd.home import home
from os.path import join
config(viewer='ngl')
2024-12-05 13:52:37,309 - numexpr.utils - INFO - Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
2024-12-05 13:52:37,310 - 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+23.g4c98f4164.dirty).
Using docking to generate starting poses for simulations#
Download the files for this tutorial from this link
Dock the protein with the ligand#
datadir = join(home(dataDir="building-protein-ligand"))
prot = Molecule(join(datadir, "trypsin.pdb"))
prot.center()
lig = Molecule(join(datadir, "BEN.cif"))
poses, scores = dock(prot, lig)
==============================
*** Open Babel Warning in PerceiveBondOrders
Failed to kekulize aromatic bonds in OBMol::PerceiveBondOrders (title is /tmp/tmpmhkzxhom/protein.pdb)
1 molecule converted
2024-12-05 13:52:42,850 - htmd.dock - INFO - Charges detected in ligand and will be used for docking.
1 molecule converted
################################################################# # If you used AutoDock Vina in your work, please cite: # # # # O. Trott, A. J. Olson, # # AutoDock Vina: improving the speed and accuracy of docking # # with a new scoring function, efficient optimization and # # multithreading, Journal of Computational Chemistry 31 (2010) # # 455-461 # # # # DOI 10.1002/jcc.21334 # # # # Please see http://vina.scripps.edu for more information. # ################################################################# WARNING: The search space volume > 27000 Angstrom^3 (See FAQ) Detected 20 CPUs WARNING: at low exhaustiveness, it may be impossible to utilize all CPUs Reading input ... done. Setting up the scoring function ... done. Analyzing the binding site ... done. Using random seed: 1238370264 Performing search ... 0% 10 20 30 40 50 60 70 80 90 100% |----|----|----|----|----|----|----|----|----|----| *********************************************** done. Refining results ... done. mode | affinity | dist from best mode | (kcal/mol) | rmsd l.b.| rmsd u.b. -----+------------+----------+---------- 1 -4.9 0.000 0.000 2 -4.9 0.042 1.620 3 -4.6 24.177 24.823 4 -4.5 30.437 31.604 5 -4.4 21.893 23.170 6 -4.3 28.748 29.753 7 -4.3 29.399 30.379 8 -4.3 2.472 3.609 9 -4.3 21.328 22.539 10 -4.3 21.329 22.516 11 -4.2 21.484 22.800 12 -4.1 29.504 30.504 13 -4.1 18.311 18.874 14 -4.1 24.407 25.335 15 -4.1 3.897 4.965 16 -4.0 24.067 25.048 17 -4.0 22.026 23.022 18 -3.8 22.494 23.397 19 -3.7 17.615 18.842 20 -3.7 30.998 32.036 Writing output ... done.
20 molecules converted
20 files output. The first is /tmp/tmpmhkzxhom/output_1.pdb
Visualize the docked poses#
mol = Molecule()
mol.append(prot)
for i, p in enumerate(poses):
mol.append(p)
mol.view(sel='protein', style='NewCartoon', hold=True)
mol.view(sel='resname BEN', style='Licorice', color=1)
NGLWidget()
Build systems from docked poses#
molbuilt = []
for i, p in enumerate(poses):
prot = Molecule(join(datadir, "trypsin.pdb"))
prot.filter('chain A and (protein or water or resname CA)')
prot.set('segid', 'P', sel='protein and noh')
prot.set('segid', 'W', sel='water')
prot.set('segid', 'CA', sel='resname CA')
prot.center()
D = maxDistance(prot, 'all')
ligand = p
ligand.set('segid','L')
mol = Molecule(name='combo')
mol.append(prot)
mol.append(ligand)
D = D + 15
smol = solvate(mol, minmax=[[-D, -D, -D], [D, D, D]])
topos = [join(datadir, "BEN.cif")]
params = [join(datadir, "BEN.frcmod")]
molbuilt.append(amber.build(smol, topo=topos, param=params, outdir=f'./docked/build/{i+1}/', saltconc=0.15))
if i==1: # For demonstration purposes lets only build the two first
break
2024-12-05 13:53:44,922 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.
2024-12-05 13:53:44,952 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-12-05 13:53:45,252 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
2024-12-05 13:53:46,435 - htmd.builder.solvate - INFO - 19297 water molecules were added to the system.
2024-12-05 13:53:48,964 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-12-05 13:53:48,981 - htmd.builder.builder - INFO - 6 disulfide bonds were added
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 26, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 42, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 117, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 184, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 174, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 198, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 149, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 163, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 8, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 138, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 110, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 211, insertion: '', segid: 'P'>
2024-12-05 13:53:50,234 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 13:53:53,664 - htmd.builder.amber - INFO - Finished building.
2024-12-05 13:53:54,470 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-12-05 13:53:59,698 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 108 ions for the given salt concentration 0.15 M.
2024-12-05 13:54:11,829 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 13:54:14,825 - htmd.builder.amber - INFO - Finished building.
2024-12-05 13:54:15,564 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
/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:54:22,933 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224
2024-12-05 13:54:23,071 - moleculekit.molecule - INFO - Removed 9 atoms. 1692 atoms remaining in the molecule.
2024-12-05 13:54:23,100 - htmd.builder.solvate - INFO - Using water pdb file at: /home/sdoerr/Work/htmd/htmd/share/solvate/wat.pdb
2024-12-05 13:54:23,346 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2
2024-12-05 13:54:24,556 - htmd.builder.solvate - INFO - 19297 water molecules were added to the system.
2024-12-05 13:54:27,067 - htmd.builder.amber - INFO - Detecting disulfide bonds.
2024-12-05 13:54:27,075 - htmd.builder.builder - INFO - 6 disulfide bonds were added
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 26, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 42, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 117, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 184, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 174, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 198, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 149, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 163, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 8, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 138, insertion: '', segid: 'P'>
Disulfide Bond between: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 110, insertion: '', segid: 'P'>
and: UniqueResidueID<resname: 'CYS', chain: 'A', resid: 211, insertion: '', segid: 'P'>
2024-12-05 13:54:28,367 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 13:54:31,500 - htmd.builder.amber - INFO - Finished building.
2024-12-05 13:54:32,248 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-12-05 13:54:37,561 - htmd.builder.ionize - INFO - Adding 9 anions + 0 cations for neutralizing and 108 ions for the given salt concentration 0.15 M.
2024-12-05 13:54:49,812 - htmd.builder.amber - INFO - Starting the build.
2024-12-05 13:54:52,936 - htmd.builder.amber - INFO - Finished building.
2024-12-05 13:54:53,745 - moleculekit.writers - WARNING - Field "resid" of PDB overflows. Your data will be truncated to 4 characters.
2024-12-05 13:55:01,315 - moleculekit.tools.sequencestructuralalignment - INFO - Alignment #0 was done on 223 residues: 2-224
Equilibrate the build systems#
from acemd.protocols import setup_equilibration
builds = sorted(glob('docked/build/*/'))
for i, bd in enumerate(builds):
setup_equilibration(bd, f'docked/equil/{i+1}/', run="10ns", temperature=298)
2024-12-05 13:55:26,560 - acemd - INFO - Copying docked/build/1/structure.pdb to docked/equil/1/structure.pdb
2024-12-05 13:55:26,564 - acemd - INFO - Copying docked/build/1/structure.prmtop to docked/equil/1/structure.prmtop
2024-12-05 13:55:29,000 - acemd - INFO - Copying docked/build/2/structure.pdb to docked/equil/2/structure.pdb
2024-12-05 13:55:29,003 - acemd - INFO - Copying docked/build/2/structure.prmtop to docked/equil/2/structure.prmtop
mdx = LocalGPUQueue()
mdx.submit(glob('./docked/equil/*/'))
mdx.wait()
Create the production folder#
from acemd.protocols import setup_production
equils = sorted(glob('docked/equil/*/'))
for i, eq in enumerate(equils):
setup_production(eq, f'docked/generators/{i+1}/', run="50ns", temperature=298)
mdx = LocalGPUQueue()
mdx.submit(glob('./docked/generators/*/'))
mdx.wait()