Protein Ligand Complex MD Setup tutorial using BioExcel Building Blocks (biobb)¶
Based on the official Gromacs tutorial: http://www.mdtutorials.com/gmx/complex/index.html¶
This tutorial aims to illustrate the process of setting up a simulation system containing a protein in complex with a ligand, step by step, using the BioExcel Building Blocks library (biobb). The particular example used is the T4 lysozyme L99A/M102Q protein (PDB code 3HTB), in complex with the 2-propylphenol small molecule (3-letter Code JZ4).
Biobb modules used:
- biobb_io: Tools to fetch biomolecular data from public databases.
- biobb_model: Tools to model macromolecular structures.
- biobb_chemistry: Tools to manipulate chemical data.
- biobb_gromacs: Tools to setup and run Molecular Dynamics simulations.
- biobb_analysis: Tools to analyse Molecular Dynamics trajectories.
- biobb_structure_utils: Tools to modify or extract information from a PDB structure file.
Auxiliar libraries used:
- nb_conda_kernels: Enables a Jupyter Notebook or JupyterLab application in one conda environment to access kernels for Python, R, and other languages found in other environments.
- nglview: Jupyter/IPython widget to interactively view molecular structures and trajectories in notebooks.
- ipywidgets: Interactive HTML widgets for Jupyter notebooks and the IPython kernel.
- os: Python miscellaneous operating system interfaces
- plotly: Python interactive graphing library integrated in Jupyter notebooks.
- simpletraj: Lightweight coordinate-only trajectory reader based on code from GROMACS, MDAnalysis and VMD.
Conda Installation and Launch¶
git clone https://github.com/bioexcel/biobb_wf_protein-complex_md_setup.git
cd biobb_wf_protein-complex_md_setup
conda env create -f conda_env/environment.yml
conda activate biobb_Protein-Complex_MDsetup_tutorial
jupyter nbextension enable python-markdown/main
jupyter-notebook biobb_wf_protein-complex_md_setup/notebooks/biobb_Protein-Complex_MDsetup_tutorial.ipynb
Please execute the following commands before launching the Jupyter Notebook if you experience some
jupyter-nbextension enable --py --user widgetsnbextension
jupyter-nbextension enable --py --user nglview
Pipeline steps:¶
- Input Parameters
- Fetching PDB Structure
- Fix Protein Structure
- Create Protein System Topology
- Create ligand system topology
- Preparing Ligand Restraints
- Create new protein-ligand complex structure file
- Create new protein-ligand complex topology file
- Create Solvent Box
- Fill the Box with Water Molecules
- Adding Ions
- Energetically Minimize the System
- Equilibrate the System (NVT)
- Equilibrate the System (NPT)
- Free Molecular Dynamics Simulation
- Post-processing and Visualizing Resulting 3D Trajectory
- Output Files
- Questions & Comments
Input parameters¶
Input parameters needed:
- pdbCode: PDB code of the protein-ligand complex structure (e.g. 3HTB)
- ligandCode: Small molecule 3-letter code for the ligand structure (e.g. JZ4)
- mol_charge: Charge of the small molecule, needed to add hydrogen atoms.
In [1]:
import nglview
import ipywidgets
import os
import zipfile
pdbCode = "3HTB"
ligandCode = "JZ4"
mol_charge = 0
Fetching PDB structure¶
Downloading PDB structure with the protein-ligand complex from the RCSB PDB database.
Alternatively, a PDB file can be used as starting structure.
Splitting the molecule in three different files:
- proteinFile: Protein structure
- ligandFile: Ligand structure
- complexFile: Protein-ligand complex structure
Building Blocks used:
- Pdb from biobb_io.api.pdb
In [2]:
# Downloading desired PDB file
# Import module
from biobb_io.api.pdb import pdb
# Create properties dict and inputs/outputs
downloaded_pdb = pdbCode+'.orig.pdb'
prop = {
'pdb_code': pdbCode,
'filter': False
}
# Create and launch bb
pdb(output_pdb_path=downloaded_pdb,
properties=prop)
2022-09-15 12:05:53,980 [MainThread ] [INFO ] Executing biobb_io.api.pdb Version: 3.8.0 2022-09-15 12:05:53,980 [MainThread ] [INFO ] Downloading 3htb from: https://www.ebi.ac.uk/pdbe/entry-files/download/pdb3htb.ent 2022-09-15 12:05:54,316 [MainThread ] [INFO ] Writting pdb to: 3HTB.orig.pdb
Out[2]:
0
In [3]:
# Extracting Protein, Ligand and Protein-Ligand Complex to three different files
# Import module
from biobb_structure_utils.utils.extract_heteroatoms import extract_heteroatoms
from biobb_structure_utils.utils.extract_molecule import extract_molecule
from biobb_structure_utils.utils.cat_pdb import cat_pdb
# Create properties dict and inputs/outputs
proteinFile = pdbCode+'.pdb'
ligandFile = ligandCode+'.pdb'
complexFile = pdbCode+'_'+ligandCode+'.pdb'
prop = {
'heteroatoms' : [{"name": "JZ4"}]
}
extract_heteroatoms(input_structure_path=downloaded_pdb,
output_heteroatom_path=ligandFile,
properties=prop)
extract_molecule(input_structure_path=downloaded_pdb,
output_molecule_path=proteinFile)
print(proteinFile, ligandFile, complexFile)
cat_pdb(input_structure1=proteinFile,
input_structure2=ligandFile,
output_structure_path=complexFile)
2022-09-15 12:05:54,584 [MainThread ] [INFO ] Executing biobb_structure_utils.utils.extract_heteroatoms Version: 3.8.0 2022-09-15 12:05:54,620 [MainThread ] [INFO ] Writting pdb to: JZ4.pdb 2022-09-15 12:05:54,621 [MainThread ] [INFO ] Removed: [] 2022-09-15 12:05:54,623 [MainThread ] [INFO ] Executing biobb_structure_utils.utils.extract_molecule Version: 3.8.0 2022-09-15 12:05:54,624 [MainThread ] [INFO ] Creating 0f220499-2c15-487b-81f0-6533cfad03a9 temporary folder 2022-09-15 12:05:54,625 [MainThread ] [INFO ] check_structure -i /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB.orig.pdb -o 3HTB.pdb --force_save --non_interactive command_list --list 0f220499-2c15-487b-81f0-6533cfad03a9/extract_prot.lst 2022-09-15 12:05:55,680 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:05:55,681 [MainThread ] [INFO ] ================================================================================ = BioBB structure checking utility v3.10.1 = = P. Andrio, A. Hospital, G. Bayarri, J.L. Gelpi 2018-22 = ================================================================================ Warning: sequence features only available in mmCIF format or with external fasta input Structure /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB.orig.pdb loaded Title: 2-propylphenol in complex with t4 lysozyme l99a/m102q Experimental method: x-ray diffraction Keywords: hydrolase, glycosidase, bacteriolytic enzyme, antimicrobial Resolution (A): 1.81 Num. models: 1 Num. chains: 1 (A: Protein) Num. residues: 387 Num. residues with ins. codes: 0 Num. residues with H atoms: 0 Num. HETATM residues: 224 Num. ligands or modified residues: 4 Num. water mol.: 220 Num. atoms: 1544 Small mol ligands found PO4 A165 PO4 A166 JZ4 A167 BME A168 Step 1: ligands --remove All Running ligands. Options: --remove All 4 Ligands detected PO4 A165 PO4 A166 JZ4 A167 BME A168 Ligands removed All (4) Step 2: water --remove Yes Running water. Options: --remove Yes 220 Water molecules detected 220 Water molecules removed Command list completed Final Num. models: 1 Final Num. chains: 1 (A: Protein) Final Num. residues: 163 Final Num. residues with ins. codes: 0 Final Num. residues with H atoms: 0 Final Num. HETATM residues: 0 Final Num. ligands or modified residues: 0 Final Num. water mol.: 0 Final Num. atoms: 1300 Structure saved on 3HTB.pdb 2022-09-15 12:05:55,682 [MainThread ] [INFO ] 2022-09-15 12:05:55,695 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:05:55,696 [MainThread ] [INFO ] Removed: [] 3HTB.pdb JZ4.pdb 3HTB_JZ4.pdb 2022-09-15 12:05:55,700 [MainThread ] [INFO ] Executing biobb_structure_utils.utils.cat_pdb Version: 3.8.0 2022-09-15 12:05:55,702 [MainThread ] [INFO ] Removed: []
Out[3]:
0
Visualizing 3D structures¶
Visualizing the generated PDB structures using NGL:
- Protein structure (Left)
- Ligand structure (Center)
- Protein-ligand complex (Right)
In [4]:
# Show structures: protein, ligand and protein-ligand complex
view1 = nglview.show_structure_file(proteinFile)
view1._remote_call('setSize', target='Widget', args=['350px','400px'])
view1.camera='orthographic'
view1
view2 = nglview.show_structure_file(ligandFile)
view2.add_representation(repr_type='ball+stick')
view2._remote_call('setSize', target='Widget', args=['350px','400px'])
view2.camera='orthographic'
view2
view3 = nglview.show_structure_file(complexFile)
view3.add_representation(repr_type='licorice', radius='.5', selection=ligandCode)
view3._remote_call('setSize', target='Widget', args=['350px','400px'])
view3.camera='orthographic'
view3
ipywidgets.HBox([view1, view2, view3])
HBox(children=(NGLWidget(), NGLWidget(), NGLWidget()))
Fix protein structure¶
Checking and fixing (if needed) the protein structure:
- Modeling missing side-chain atoms, modifying incorrect amide assignments, choosing alternative locations.
- Checking for missing backbone atoms, heteroatoms, modified residues and possible atomic clashes.
Building Blocks used:
- FixSideChain from biobb_model.model.fix_side_chain
In [5]:
# Check & Fix Protein Structure
# Import module
from biobb_model.model.fix_side_chain import fix_side_chain
# Create prop dict and inputs/outputs
fixed_pdb = pdbCode+'_fixed.pdb'
# Create and launch bb
fix_side_chain(input_pdb_path=proteinFile,
output_pdb_path=fixed_pdb)
2022-09-15 12:05:55,796 [MainThread ] [INFO ] Executing biobb_model.model.fix_side_chain Version: 3.8.0 2022-09-15 12:05:55,797 [MainThread ] [INFO ] check_structure -i 3HTB.pdb -o 3HTB_fixed.pdb --force_save fixside --fix ALL 2022-09-15 12:05:56,091 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:05:56,092 [MainThread ] [INFO ] ================================================================================ = BioBB structure checking utility v3.10.1 = = P. Andrio, A. Hospital, G. Bayarri, J.L. Gelpi 2018-22 = ================================================================================ Warning: sequence features only available in mmCIF format or with external fasta input Structure 3HTB.pdb loaded Title: Experimental method: unknown Resolution (A): N.A. Num. models: 1 Num. chains: 1 (A: Protein) Num. residues: 163 Num. residues with ins. codes: 0 Num. residues with H atoms: 0 Num. HETATM residues: 0 Num. ligands or modified residues: 0 Num. water mol.: 0 Num. atoms: 1300 Running fixside. Options: --fix ALL No residues with missing or unknown side chain atoms found Structure not modified, saving due to --force_save option Final Num. models: 1 Final Num. chains: 1 (A: Protein) Final Num. residues: 163 Final Num. residues with ins. codes: 0 Final Num. residues with H atoms: 0 Final Num. HETATM residues: 0 Final Num. ligands or modified residues: 0 Final Num. water mol.: 0 Final Num. atoms: 1300 Structure saved on 3HTB_fixed.pdb 2022-09-15 12:05:56,093 [MainThread ] [INFO ] Removed: []
Out[5]:
0
Create protein system topology¶
Building GROMACS topology corresponding to the protein structure.
Force field used in this tutorial is amber99sb-ildn: AMBER parm99 force field with corrections on backbone (sb) and side-chain torsion potentials (ildn). Water molecules type used in this tutorial is spc/e.
Adding hydrogen atoms if missing. Automatically identifying disulfide bridges.
Generating two output files:
- GROMACS structure (gro file)
- GROMACS topology ZIP compressed file containing:
- GROMACS topology top file (top file)
- GROMACS position restraint file/s (itp file/s)
Building Blocks used:
- Pdb2gmx from biobb_gromacs.gromacs.pdb2gmx
In [6]:
# Create Protein system topology
# Import module
from biobb_gromacs.gromacs.pdb2gmx import pdb2gmx
# Create inputs/outputs
output_pdb2gmx_gro = pdbCode+'_pdb2gmx.gro'
output_pdb2gmx_top_zip = pdbCode+'_pdb2gmx_top.zip'
prop = {
'force_field' : 'amber99sb-ildn',
'water_type': 'spce'
}
# Create and launch bb
pdb2gmx(input_pdb_path=fixed_pdb,
output_gro_path=output_pdb2gmx_gro,
output_top_zip_path=output_pdb2gmx_top_zip,
properties=prop)
2022-09-15 12:05:57,545 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.pdb2gmx Version: 3.8.0
2022-09-15 12:05:57,546 [MainThread ] [INFO ] GROMACS Pdb2gmx 20222 version detected
2022-09-15 12:05:57,547 [MainThread ] [INFO ] gmx -nobackup -nocopyright pdb2gmx -f 3HTB_fixed.pdb -o 3HTB_pdb2gmx.gro -p p2g.top -water spce -ff amber99sb-ildn -i posre.itp
2022-09-15 12:05:57,862 [MainThread ] [INFO ] Exit code 0
2022-09-15 12:05:57,864 [MainThread ] [INFO ] Using the Amber99sb-ildn force field in directory amber99sb-ildn.ff
going to rename amber99sb-ildn.ff/aminoacids.r2b
going to rename amber99sb-ildn.ff/dna.r2b
going to rename amber99sb-ildn.ff/rna.r2b
Reading 3HTB_fixed.pdb...
Read '', 1364 atoms
Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.
There are 1 chains and 0 blocks of water and 163 residues with 1364 atoms
chain #res #atoms
1 'A' 163 1364
there were 0 atoms with zero occupancy and 178 atoms with occupancy unequal to one (out of 1364 atoms). Check your pdb file.
Reading residue database... (Amber99sb-ildn)
Processing chain 1 'A' (1364 atoms, 163 residues)
Identified residue MET1 as a starting terminus.
Identified residue ASN163 as a ending terminus.
Checking for duplicate atoms....
Now there are 1300 atoms. Deleted 64 duplicates.
Generating any missing hydrogen atoms and/or adding termini.
Now there are 163 residues with 2614 atoms
Making bonds...
Number of bonds was 2635, now 2634
Generating angles, dihedrals and pairs...
Making cmap torsions...
There are 7252 dihedrals, 522 impropers, 4751 angles
6847 pairs, 2634 bonds and 0 virtual sites
Total mass 18512.384 a.m.u.
Total charge 6.000 e
Writing topology
Writing coordinate file...
--------- PLEASE NOTE ------------
You have successfully generated a topology from: 3HTB_fixed.pdb.
The Amber99sb-ildn force field and the spce water model are used.
--------- ETON ESAELP ------------
2022-09-15 12:05:57,864 [MainThread ] [INFO ] :-) GROMACS - gmx pdb2gmx, 2022.2-conda_forge (-:
Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx
Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial
Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks
Command line:
gmx -nobackup -nocopyright pdb2gmx -f 3HTB_fixed.pdb -o 3HTB_pdb2gmx.gro -p p2g.top -water spce -ff amber99sb-ildn -i posre.itp
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/aminoacids.r2b
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/dna.r2b
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/rna.r2b
there were 0 atoms with zero occupancy and 178 atoms with occupancy unequal to one (out of 1364 atoms). Check your pdb file.
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/atomtypes.atp
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/aminoacids.rtp
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/dna.rtp
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/rna.rtp
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/aminoacids.hdb
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/dna.hdb
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/rna.hdb
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/aminoacids.n.tdb
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/aminoacids.c.tdb
Analysing hydrogen-bonding network for automated assignment of histidine
protonation. 273 donors and 248 acceptors were found.
There are 347 hydrogen bonds
Will use HISD for residue 31
8 out of 8 lines of specbond.dat converted successfully
Special Atom Distance matrix:
MET1 MET6 HIS31 CYS54 CYS97 MET106
SD10 SD56 NE2276 SG446 SG808 SD883
MET6 SD56 0.631
HIS31 NE2276 2.196 1.931
CYS54 SG446 2.781 2.718 1.052
CYS97 SG808 0.820 0.539 2.117 2.897
MET106 SD883 1.900 1.314 1.848 2.841 1.635
MET120 SD982 2.525 2.014 3.286 4.293 2.112 1.520
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/aminoacids.arn
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/dna.arn
Opening force field file /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/share/gromacs/top/amber99sb-ildn.ff/rna.arn
Before cleaning: 6889 pairs
Before cleaning: 7252 dihedrals
GROMACS reminds you: "Wedged as we are between two eternities of idleness, there is no excuse for being idle now" (Anthony Burgess)
2022-09-15 12:05:57,865 [MainThread ] [INFO ] Compressing topology to: 3HTB_pdb2gmx_top.zip
2022-09-15 12:05:57,866 [MainThread ] [INFO ] Ignored file amber99sb-ildn.ff/forcefield.itp
2022-09-15 12:05:57,873 [MainThread ] [INFO ] Ignored file amber99sb-ildn.ff/spce.itp
2022-09-15 12:05:57,873 [MainThread ] [INFO ] Ignored file amber99sb-ildn.ff/ions.itp
2022-09-15 12:05:57,876 [MainThread ] [INFO ] Adding:
2022-09-15 12:05:57,876 [MainThread ] [INFO ] ['p2g.top', 'posre.itp']
2022-09-15 12:05:57,877 [MainThread ] [INFO ] to: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx_top.zip
2022-09-15 12:05:57,878 [MainThread ] [INFO ] Removed: ['p2g.top', 'posre.itp']
Out[6]:
0
Create ligand system topology¶
Building GROMACS topology corresponding to the ligand structure.
Force field used in this tutorial step is amberGAFF: General AMBER Force Field, designed for rational drug design.
- Step 1: Add hydrogen atoms if missing.
- Step 2: Energetically minimize the system with the new hydrogen atoms.
- Step 3: Generate ligand topology (parameters).
Building Blocks used:
- ReduceAddHydrogens from biobb_chemistry.ambertools.reduce_add_hydrogens
- BabelMinimize from biobb_chemistry.babelm.babel_minimize
- AcpypeParamsGMX from biobb_chemistry.acpype.acpype_params_gmx
Step 1: Add hydrogen atoms¶
In [7]:
# Create Ligand system topology, STEP 1
# Reduce_add_hydrogens: add Hydrogen atoms to a small molecule (using Reduce tool from Ambertools package)
# Import module
from biobb_chemistry.ambertools.reduce_add_hydrogens import reduce_add_hydrogens
# Create prop dict and inputs/outputs
output_reduce_h = ligandCode+'.reduce.H.pdb'
prop = {
'nuclear' : 'true'
}
# Create and launch bb
reduce_add_hydrogens(input_path=ligandFile,
output_path=output_reduce_h,
properties=prop)
2022-09-15 12:05:57,888 [MainThread ] [INFO ] reduce -NUClear -OH -ROTNH3 -ALLALT /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4.pdb > JZ4.reduce.H.pdb 2022-09-15 12:05:59,084 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:05:59,085 [MainThread ] [INFO ] reduce: version 3.3 06/02/2016, Copyright 1997-2016, J. Michael Word Processing file: "/Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4.pdb" Database of HETATM connections: "/Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial//dat/reduce_wwPDB_het_dict.txt " VDW dot density = 16/A^2 Orientation penalty scale = 1 (100%) Eliminate contacts within 3 bonds. Ignore atoms with |occupancy| <= 0.01 during adjustments. Waters ignored if B-Factor >= 40 or |occupancy| < 0.66 Aromatic rings in amino acids accept hydrogen bonds. Building or keeping OH & SH Hydrogens. Rotating NH3 Hydrogens. Not processing Met methyls. WARNING: atom H13A from JZ4 will be treated as hydrogen WARNING: atom H14A from JZ4 will be treated as hydrogen WARNING: atom H13A from JZ4 will be treated as hydrogen WARNING: atom H13A from JZ4 will be treated as hydrogen WARNING: atom H14A from JZ4 will be treated as hydrogen WARNING: atom H14A from JZ4 will be treated as hydrogen WARNING: atom H14A from JZ4 will be treated as hydrogen WARNING: atom H14A from JZ4 will be treated as hydrogen WARNING: atom H13A from JZ4 will be treated as hydrogen WARNING: atom H13A from JZ4 will be treated as hydrogen WARNING: atom H14A from JZ4 will be treated as hydrogen WARNING: atom H13A from JZ4 will be treated as hydrogen Singles(size 1): A 167 JZ4 OAB orientation 1: A 167 JZ4 OAB : rot 120: bump=0.000, HB=0.000, total=0.000 Found 0 hydrogens (0 hets) Standardized 0 hydrogens (0 hets) Added 12 hydrogens (12 hets) Removed 0 hydrogens (0 hets) Adjusted 1 group(s) If you publish work which uses reduce, please cite: Word, et. al. (1999) J. Mol. Biol. 285, 1735-1747. For more information see http://kinemage.biochem.duke.edu
Out[7]:
0
Step 2: Energetically minimize the system with the new hydrogen atoms.¶
In [8]:
# Create Ligand system topology, STEP 2
# Babel_minimize: Structure energy minimization of a small molecule after being modified adding hydrogen atoms
# Import module
from biobb_chemistry.babelm.babel_minimize import babel_minimize
# Create prop dict and inputs/outputs
output_babel_min = ligandCode+'.H.min.mol2'
prop = {
'method' : 'sd',
'criteria' : '1e-10',
'force_field' : 'GAFF'
}
# Create and launch bb
babel_minimize(input_path=output_reduce_h,
output_path=output_babel_min,
properties=prop)
2022-09-15 12:05:59,095 [MainThread ] [INFO ] Hydrogens is not correct, assigned default value: False
2022-09-15 12:05:59,096 [MainThread ] [INFO ] Steps is not correct, assigned default value: 2500
2022-09-15 12:05:59,096 [MainThread ] [INFO ] Cut-off is not correct, assigned default value: False
2022-09-15 12:05:59,097 [MainThread ] [INFO ] Rvdw is not correct, assigned default value: 6.0
2022-09-15 12:05:59,097 [MainThread ] [INFO ] Rele is not correct, assigned default value: 10.0
2022-09-15 12:05:59,098 [MainThread ] [INFO ] Frequency is not correct, assigned default value: 10
2022-09-15 12:05:59,098 [MainThread ] [INFO ] obminimize -c 1e-10 -sd -ff GAFF -ipdb /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4.reduce.H.pdb -omol2 > JZ4.H.min.mol2
2022-09-15 12:06:24,611 [MainThread ] [INFO ] Exit code 0
2022-09-15 12:06:24,612 [MainThread ] [INFO ]
A T O M T Y P E S
IDX TYPE RING
1 c3 NO
2 ca AR
3 ca AR
4 ca AR
5 ca AR
6 ca AR
7 ca AR
8 c3 NO
9 c3 NO
10 oh NO
11 ho NO
12 hc NO
13 hc NO
14 ha NO
15 ha NO
16 ha NO
17 hc NO
18 hc NO
19 hc NO
20 hc NO
21 hc NO
22 ha NO
C H A R G E S
IDX CHARGE
1 -0.064979
2 -0.061278
3 -0.058294
4 -0.019907
5 0.120013
6 -0.055115
7 -0.005954
8 -0.024481
9 -0.051799
10 -0.506496
11 0.292144
12 0.026577
13 0.031423
14 0.065403
15 0.061871
16 0.061769
17 0.022987
18 0.022987
19 0.022987
20 0.026577
21 0.031423
22 0.062142
S E T T I N G U P C A L C U L A T I O N S
SETTING UP BOND CALCULATIONS...
SETTING UP ANGLE CALCULATIONS...
SETTING UP TORSION CALCULATIONS...
SETTING UP IMPROPER TORSION CALCULATIONS...
SETTING UP VAN DER WAALS CALCULATIONS...
SETTING UP ELECTROSTATIC CALCULATIONS...
S T E E P E S T D E S C E N T
STEPS = 2500
STEP n E(n) E(n-1)
------------------------------------
0 50.048 ----
10 36.21724 36.73363
20 32.62519 32.90050
30 30.42176 30.60740
40 28.84022 28.98002
50 27.60107 27.71398
60 26.57419 26.66957
70 25.69242 25.77532
80 24.91809 24.99145
90 24.22825 24.29392
100 23.60848 23.66758
110 23.04959 23.10302
120 22.54147 22.59026
130 22.07637 22.12111
140 21.64898 21.69017
150 21.25484 21.29288
160 20.89016 20.92541
170 20.55171 20.58447
180 20.23670 20.26723
190 19.94278 19.97129
200 19.66790 19.69459
210 19.41031 19.43534
220 19.16849 19.19201
230 18.94114 18.96327
240 18.72711 18.74795
250 18.52539 18.54504
260 18.33509 18.35363
270 18.15540 18.17292
280 17.98561 18.00217
290 17.82508 17.84074
300 17.67321 17.68803
310 17.52946 17.54349
320 17.39333 17.40662
330 17.26436 17.27695
340 17.14213 17.15406
350 17.02624 17.03755
360 16.91631 16.92705
370 16.81202 16.82221
380 16.71303 16.72270
390 16.61905 16.62823
400 16.52979 16.53852
410 16.44500 16.45329
420 16.36442 16.37230
430 16.28783 16.29532
440 16.21501 16.22213
450 16.14576 16.15253
460 16.07988 16.08632
470 16.01719 16.02332
480 15.95753 15.96337
490 15.90074 15.90630
500 15.84667 15.85196
510 15.79517 15.80021
520 15.74612 15.75092
530 15.69939 15.70397
540 15.65487 15.65923
550 15.61243 15.61659
560 15.57198 15.57594
570 15.53342 15.53720
580 15.49665 15.50025
590 15.46159 15.46502
600 15.42815 15.43142
610 15.39625 15.39937
620 15.36581 15.36879
630 15.33678 15.33962
640 15.30907 15.31178
650 15.28262 15.28521
660 15.25738 15.25985
670 15.23328 15.23564
680 15.21028 15.21253
690 15.18831 15.19046
700 15.16733 15.16939
710 15.14730 15.14926
720 15.12817 15.13004
730 15.10989 15.11168
740 15.09242 15.09414
750 15.07574 15.07737
760 15.05979 15.06136
770 15.04455 15.04605
780 15.02999 15.03142
790 15.01606 15.01743
800 15.00275 15.00405
810 14.99002 14.99127
820 14.97784 14.97903
830 14.96619 14.96733
840 14.95505 14.95614
850 14.94438 14.94543
860 14.93418 14.93518
870 14.92441 14.92537
880 14.91506 14.91597
890 14.90610 14.90698
900 14.89752 14.89836
910 14.88930 14.89010
920 14.88142 14.88219
930 14.87387 14.87461
940 14.86662 14.86733
950 14.85968 14.86036
960 14.85301 14.85367
970 14.84661 14.84724
980 14.84047 14.84107
990 14.83457 14.83515
1000 14.82890 14.82945
1010 14.82344 14.82398
1020 14.81820 14.81872
1030 14.81315 14.81365
1040 14.80830 14.80877
1050 14.80361 14.80408
1060 14.79910 14.79955
1070 14.79475 14.79518
1080 14.79055 14.79097
1090 14.78650 14.78690
1100 14.78257 14.78296
1110 14.77878 14.77916
1120 14.77511 14.77547
1130 14.77155 14.77190
1140 14.76818 14.76848
1150 14.76618 14.76634
1160 14.76457 14.76473
1170 14.76301 14.76316
1180 14.76147 14.76162
1190 14.75998 14.76012
1200 14.75851 14.75866
1210 14.75708 14.75722
1220 14.75568 14.75582
1230 14.75431 14.75445
1240 14.75297 14.75310
1250 14.75166 14.75179
1260 14.75037 14.75050
1270 14.74911 14.74924
1280 14.74788 14.74800
1290 14.74667 14.74679
1300 14.74548 14.74560
1310 14.74432 14.74443
1320 14.74317 14.74329
1330 14.74205 14.74216
1340 14.74095 14.74106
1350 14.73987 14.73998
1360 14.73881 14.73891
1370 14.73776 14.73787
1380 14.73673 14.73684
1390 14.73572 14.73582
1400 14.73473 14.73483
1410 14.73375 14.73385
1420 14.73278 14.73288
1430 14.73183 14.73193
1440 14.73090 14.73099
1450 14.72997 14.73006
1460 14.72906 14.72915
1470 14.72816 14.72825
1480 14.72727 14.72736
1490 14.72639 14.72648
1500 14.72553 14.72561
1510 14.72493 14.72497
1520 14.72450 14.72454
1530 14.72408 14.72413
1540 14.72367 14.72371
1550 14.72327 14.72331
1560 14.72287 14.72291
1570 14.72247 14.72251
1580 14.72208 14.72212
1590 14.72170 14.72174
1600 14.72132 14.72136
1610 14.72095 14.72099
1620 14.72058 14.72062
1630 14.72022 14.72025
1640 14.71986 14.71989
1650 14.71950 14.71954
1660 14.71915 14.71919
1670 14.71881 14.71884
1680 14.71847 14.71850
1690 14.71813 14.71817
1700 14.71780 14.71783
1710 14.71747 14.71750
1720 14.71715 14.71718
1730 14.71683 14.71686
1740 14.71651 14.71654
1750 14.71619 14.71623
1760 14.71588 14.71592
1770 14.71558 14.71561
1780 14.71527 14.71530
1790 14.71497 14.71500
1800 14.71468 14.71471
1810 14.71438 14.71441
1820 14.71409 14.71412
1830 14.71380 14.71383
1840 14.71351 14.71354
1850 14.71323 14.71326
1860 14.71305 14.71306
1870 14.71290 14.71291
1880 14.71276 14.71277
1890 14.71261 14.71263
1900 14.71247 14.71249
1910 14.71233 14.71235
1920 14.71220 14.71221
1930 14.71206 14.71208
1940 14.71193 14.71194
1950 14.71180 14.71181
1960 14.71167 14.71168
1970 14.71154 14.71155
1980 14.71141 14.71143
1990 14.71129 14.71130
2000 14.71117 14.71118
2010 14.71105 14.71106
2020 14.71093 14.71094
2030 14.71081 14.71082
2040 14.71070 14.71071
2050 14.71058 14.71059
2060 14.71047 14.71048
2070 14.71036 14.71037
2080 14.71025 14.71026
2090 14.71014 14.71015
2100 14.71004 14.71005
2110 14.70993 14.70994
2120 14.70983 14.70984
2130 14.70972 14.70973
2140 14.70962 14.70963
2150 14.70952 14.70953
2160 14.70942 14.70943
2170 14.70932 14.70933
2180 14.70922 14.70923
2190 14.70912 14.70913
2200 14.70903 14.70904
2210 14.70895 14.70895
2220 14.70889 14.70890
2230 14.70884 14.70885
2240 14.70879 14.70880
2250 14.70874 14.70875
2260 14.70870 14.70870
2270 14.70865 14.70865
2280 14.70860 14.70861
2290 14.70856 14.70856
2300 14.70851 14.70852
2310 14.70847 14.70847
2320 14.70842 14.70843
2330 14.70838 14.70839
2340 14.70834 14.70834
2350 14.70830 14.70830
2360 14.70826 14.70826
2370 14.70822 14.70822
2380 14.70818 14.70818
2390 14.70814 14.70814
2400 14.70810 14.70811
2410 14.70806 14.70807
2420 14.70803 14.70803
2430 14.70799 14.70799
2440 14.70796 14.70796
2450 14.70792 14.70792
2460 14.70788 14.70789
2470 14.70785 14.70785
2480 14.70782 14.70782
2490 14.70778 14.70779
2500 14.70775 14.70775
Time: 0.228067seconds. Iterations per second: 10966.1
Out[8]:
0
Visualizing 3D structures¶
Visualizing the small molecule generated PDB structures using NGL:
- Original Ligand Structure (Left)
- Ligand Structure with hydrogen atoms added (with Reduce program) (Center)
- Ligand Structure with hydrogen atoms added (with Reduce program), energy minimized (with Open Babel) (Right)
In [9]:
# Show different structures generated (for comparison)
view1 = nglview.show_structure_file(ligandFile)
view1.add_representation(repr_type='ball+stick')
view1._remote_call('setSize', target='Widget', args=['350px','400px'])
view1.camera='orthographic'
view1
view2 = nglview.show_structure_file(output_reduce_h)
view2.add_representation(repr_type='ball+stick')
view2._remote_call('setSize', target='Widget', args=['350px','400px'])
view2.camera='orthographic'
view2
view3 = nglview.show_structure_file(output_babel_min)
view3.add_representation(repr_type='ball+stick')
view3._remote_call('setSize', target='Widget', args=['350px','400px'])
view3.camera='orthographic'
view3
ipywidgets.HBox([view1, view2, view3])
HBox(children=(NGLWidget(), NGLWidget(), NGLWidget()))
Step 3: Generate ligand topology (parameters).¶
In [10]:
# Create Ligand system topology, STEP 3
# Acpype_params_gmx: Generation of topologies for GROMACS with ACPype
# Import module
from biobb_chemistry.acpype.acpype_params_gmx import acpype_params_gmx
# Create prop dict and inputs/outputs
output_acpype_gro = ligandCode+'params.gro'
output_acpype_itp = ligandCode+'params.itp'
output_acpype_top = ligandCode+'params.top'
output_acpype = ligandCode+'params'
prop = {
'basename' : output_acpype,
'charge' : mol_charge
}
# Create and launch bb
acpype_params_gmx(input_path=output_babel_min,
output_path_gro=output_acpype_gro,
output_path_itp=output_acpype_itp,
output_path_top=output_acpype_top,
properties=prop)
2022-09-15 12:06:24,697 [MainThread ] [INFO ] acpype -i /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4.H.min.mol2 -b JZ4params.VA80jY -n 0 2022-09-15 12:06:34,415 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:34,416 [MainThread ] [INFO ] ========================================================================== | ACPYPE: AnteChamber PYthon Parser interfacE v. 2022.6.6 (c) 2022 AWSdS | ========================================================================== ==> ... charge set to 0 ==> Executing Antechamber... ==> * Antechamber OK * ==> * Parmchk OK * ==> Executing Tleap... ==> * Tleap OK * ==> Removing temporary files... ==> Using OpenBabel v.3.1.0 ==> Writing NEW PDB file ==> Writing CNS/XPLOR files ==> Writing GROMACS files ==> Disambiguating lower and uppercase atomtypes in GMX top file, even if identical. ==> Writing GMX dihedrals for GMX 4.5 and higher. ==> Writing CHARMM files ==> Writing pickle file JZ4params.VA80jY.pkl ==> Removing temporary files... Total time of execution: 9s 2022-09-15 12:06:34,418 [MainThread ] [INFO ] File JZ4params.top succesfully created 2022-09-15 12:06:34,419 [MainThread ] [INFO ] File JZ4params.itp succesfully created 2022-09-15 12:06:34,421 [MainThread ] [INFO ] File JZ4params.gro succesfully created 2022-09-15 12:06:34,422 [MainThread ] [INFO ] Removed temporary folder: JZ4params.VA80jY.acpype
Out[10]:
0
Preparing Ligand Restraints¶
In subsequent steps of the pipeline, such as the equilibration stages of the protein-ligand complex system, it is recommended to apply some restraints to the small molecule, to avoid a possible change in position due to protein repulsion. Position restraints will be applied to the ligand, using a **force constant of 1000 KJ/mol*nm^2** on the three coordinates: x, y and z. In this steps the restriction files will be created and integrated in the ligand topology.
- Step 1: Creating an index file with a new group including just the small molecule heavy atoms.
- Step 2: Generating the position restraints file.
Building Blocks used:
Step 1: Creating an index file for the small molecule heavy atoms¶
In [11]:
# MakeNdx: Creating index file with a new group (small molecule heavy atoms)
from biobb_gromacs.gromacs.make_ndx import make_ndx
# Create prop dict and inputs/outputs
output_ligand_ndx = ligandCode+'_index.ndx'
prop = {
'selection': "0 & ! a H*"
}
# Create and launch bb
make_ndx(input_structure_path=output_acpype_gro,
output_ndx_path=output_ligand_ndx,
properties=prop)
2022-09-15 12:06:34,462 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.make_ndx Version: 3.8.0 2022-09-15 12:06:34,464 [MainThread ] [INFO ] GROMACS MakeNdx 20222 version detected 2022-09-15 12:06:34,464 [MainThread ] [INFO ] echo -e '0 & ! a H*\nq' | gmx -nobackup -nocopyright make_ndx -f JZ4params.gro -o JZ4_index.ndx 2022-09-15 12:06:34,498 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:34,500 [MainThread ] [INFO ] Going to read 0 old index file(s) Analysing residue names: There are: 1 Other residues Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups... 0 System : 22 atoms 1 Other : 22 atoms 2 JZ4 : 22 atoms nr : group '!': not 'name' nr name 'splitch' nr Enter: list groups 'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues 't': atom type '|': or 'keep' nr 'splitat' nr 'h': help 'r': residue 'res' nr 'chain' char "name": group 'case': case sensitive 'q': save and quit 'ri': residue index > Copied index group 0 'System' Found 12 atoms with name H* Complemented group: 10 atoms Merged two groups with AND: 22 10 -> 10 3 System_&_!H* : 10 atoms > 2022-09-15 12:06:34,500 [MainThread ] [INFO ] :-) GROMACS - gmx make_ndx, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx -nobackup -nocopyright make_ndx -f JZ4params.gro -o JZ4_index.ndx Reading structure file GROMACS reminds you: "Your Excellency, I have no need of this hypothesis." (Pierre Laplace, to Napoleon on why his works on celestial mechanics make no mention of God.) 2022-09-15 12:06:34,501 [MainThread ] [INFO ] Removed: []
Out[11]:
0
Step 2: Generating the position restraints file¶
In [12]:
# Genrestr: Generating the position restraints file
from biobb_gromacs.gromacs.genrestr import genrestr
# Create prop dict and inputs/outputs
output_restraints_top = ligandCode+'_posres.itp'
prop = {
'force_constants': "1000 1000 1000",
'restrained_group': "System"
}
# Create and launch bb
genrestr(input_structure_path=output_acpype_gro,
input_ndx_path=output_ligand_ndx,
output_itp_path=output_restraints_top,
properties=prop)
2022-09-15 12:06:34,539 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.genrestr Version: 3.8.0 2022-09-15 12:06:34,540 [MainThread ] [INFO ] GROMACS Genrestr 20222 version detected 2022-09-15 12:06:34,540 [MainThread ] [INFO ] echo "System" | gmx -nobackup -nocopyright genrestr -f JZ4params.gro -o JZ4_posres.itp -fc 1000 1000 1000 -n JZ4_index.ndx 2022-09-15 12:06:34,573 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:34,574 [MainThread ] [INFO ] Select group to position restrain Selected 0: 'System' 2022-09-15 12:06:34,575 [MainThread ] [INFO ] :-) GROMACS - gmx genrestr, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx -nobackup -nocopyright genrestr -f JZ4params.gro -o JZ4_posres.itp -fc 1000 1000 1000 -n JZ4_index.ndx Reading structure file Group 0 ( System) has 22 elements Group 1 ( Other) has 22 elements Group 2 ( JZ4) has 22 elements Group 3 ( System_&_!H*) has 10 elements Select a group: GROMACS reminds you: "Your Excellency, I have no need of this hypothesis." (Pierre Laplace, to Napoleon on why his works on celestial mechanics make no mention of God.) 2022-09-15 12:06:34,575 [MainThread ] [INFO ] Removed: []
Out[12]:
0
Create new protein-ligand complex structure file¶
Building new protein-ligand complex PDB file with:
- The new protein system with fixed problems from Fix Protein Structure step and hydrogens atoms added from Create Protein System Topology step.
- The new ligand system with hydrogens atoms added from Create Ligand System Topology step.
This new structure is needed for GROMACS as it is force field-compliant, it has all the new hydrogen atoms, and the atom names are matching the newly generated protein and ligand topologies.
Building Blocks used:
- GMXTrjConvStr from biobb_analysis.gromacs.gmx_trjconv_str
In [13]:
# biobb analysis module
from biobb_analysis.gromacs.gmx_trjconv_str import gmx_trjconv_str
from biobb_structure_utils.utils.cat_pdb import cat_pdb
# Convert gro (with hydrogens) to pdb (PROTEIN)
proteinFile_H = pdbCode+'_'+ligandCode+'_complex_H.pdb'
prop = {
'selection' : 'System'
}
# Create and launch bb
gmx_trjconv_str(input_structure_path=output_pdb2gmx_gro,
input_top_path=output_pdb2gmx_gro,
output_str_path=proteinFile_H,
properties=prop)
# Convert gro (with hydrogens) to pdb (LIGAND)
ligandFile_H = ligandCode+'_complex_H.pdb'
prop = {
'selection' : 'System'
}
# Create and launch bb
gmx_trjconv_str(input_structure_path=output_acpype_gro,
input_top_path=output_acpype_gro,
output_str_path=ligandFile_H,
properties=prop)
# Concatenating both PDB files: Protein + Ligand
complexFile_H = pdbCode+'_'+ligandCode+'_H.pdb'
# Create and launch bb
cat_pdb(input_structure1=proteinFile_H,
input_structure2=ligandFile_H,
output_structure_path=complexFile_H)
2022-09-15 12:06:34,585 [MainThread ] [INFO ] Executing biobb_analysis.gromacs.gmx_trjconv_str Version: 3.8.0 2022-09-15 12:06:34,585 [MainThread ] [INFO ] echo "System" | gmx trjconv -f /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx.gro -s /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx.gro -o 3HTB_JZ4_complex_H.pdb -nocenter 2022-09-15 12:06:34,629 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:34,630 [MainThread ] [INFO ] Note that major changes are planned in future for trjconv, to improve usability and utility. Select group for output Selected 0: 'System' 2022-09-15 12:06:34,631 [MainThread ] [INFO ] :-) GROMACS - gmx trjconv, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx trjconv -f /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx.gro -s /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx.gro -o 3HTB_JZ4_complex_H.pdb -nocenter Will write pdb: Protein data bank file Group 0 ( System) has 2614 elements Group 1 ( Protein) has 2614 elements Group 2 ( Protein-H) has 1301 elements Group 3 ( C-alpha) has 163 elements Group 4 ( Backbone) has 489 elements Group 5 ( MainChain) has 653 elements Group 6 ( MainChain+Cb) has 805 elements Group 7 ( MainChain+H) has 815 elements Group 8 ( SideChain) has 1799 elements Group 9 ( SideChain-H) has 648 elements Select a group: Reading frames from gro file 'GRoups of Organic Molecules in ACtion for Science', 2614 atoms. Reading frame 0 time 0.000 Precision of /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx.gro is 0.001 (nm) Last frame 0 time 0.000 -> frame 0 time 0.000 Last written: frame 0 time 0.000 GROMACS reminds you: "Your Excellency, I have no need of this hypothesis." (Pierre Laplace, to Napoleon on why his works on celestial mechanics make no mention of God.) 2022-09-15 12:06:34,633 [MainThread ] [INFO ] Executing biobb_analysis.gromacs.gmx_trjconv_str Version: 3.8.0 2022-09-15 12:06:34,634 [MainThread ] [INFO ] echo "System" | gmx trjconv -f /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4params.gro -s /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4params.gro -o JZ4_complex_H.pdb -nocenter 2022-09-15 12:06:34,668 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:34,669 [MainThread ] [INFO ] Note that major changes are planned in future for trjconv, to improve usability and utility. Select group for output Selected 0: 'System' 2022-09-15 12:06:34,669 [MainThread ] [INFO ] :-) GROMACS - gmx trjconv, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx trjconv -f /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4params.gro -s /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4params.gro -o JZ4_complex_H.pdb -nocenter Will write pdb: Protein data bank file Group 0 ( System) has 22 elements Group 1 ( Other) has 22 elements Group 2 ( JZ4) has 22 elements Select a group: Reading frames from gro file 'JZ4params_GMX.gro created by acpype (v: 2022.6.6) on Thu Sep 15 12:06:24 2022', 22 atoms. Reading frame 0 time 0.000 Precision of /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/JZ4params.gro is 0.001 (nm) Last frame 0 time 0.000 -> frame 0 time 0.000 Last written: frame 0 time 0.000 GROMACS reminds you: "Your Excellency, I have no need of this hypothesis." (Pierre Laplace, to Napoleon on why his works on celestial mechanics make no mention of God.) 2022-09-15 12:06:34,672 [MainThread ] [INFO ] Executing biobb_structure_utils.utils.cat_pdb Version: 3.8.0 2022-09-15 12:06:34,674 [MainThread ] [INFO ] Removed: []
Out[13]:
0
Create new protein-ligand complex topology file¶
Building new protein-ligand complex GROMACS topology file with:
- The new protein system topology generated from Create Protein System Topology step.
- The new ligand system topology generated from Create Ligand System Topology step.
NOTE: From this point on, the protein-ligand complex structure and topology generated can be used in a regular MD setup.
Building Blocks used:
- AppendLigand from biobb_gromacs.gromacs_extra.append_ligand (NOTE: link should be updated with the documentation)
In [14]:
# AppendLigand: Append a ligand to a GROMACS topology
# Import module
from biobb_gromacs.gromacs_extra.append_ligand import append_ligand
# Create prop dict and inputs/outputs
output_complex_top = pdbCode+'_'+ligandCode+'_complex.top.zip'
posresifdef = "POSRES_"+ligandCode.upper()
prop = {
'posres_name': posresifdef
}
# Create and launch bb
append_ligand(input_top_zip_path=output_pdb2gmx_top_zip,
input_posres_itp_path=output_restraints_top,
input_itp_path=output_acpype_itp,
output_top_zip_path=output_complex_top,
properties=prop)
2022-09-15 12:06:34,682 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs_extra.append_ligand Version: 3.8.0 2022-09-15 12:06:34,685 [MainThread ] [INFO ] Extracting: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_pdb2gmx_top.zip 2022-09-15 12:06:34,685 [MainThread ] [INFO ] to: 2022-09-15 12:06:34,686 [MainThread ] [INFO ] ['62a771b0-ba3a-4895-9132-b426ffdefa34/p2g.top', '62a771b0-ba3a-4895-9132-b426ffdefa34/posre.itp'] 2022-09-15 12:06:34,686 [MainThread ] [INFO ] Unzipping: 2022-09-15 12:06:34,686 [MainThread ] [INFO ] 3HTB_pdb2gmx_top.zip 2022-09-15 12:06:34,687 [MainThread ] [INFO ] To: 2022-09-15 12:06:34,687 [MainThread ] [INFO ] 62a771b0-ba3a-4895-9132-b426ffdefa34/p2g.top 2022-09-15 12:06:34,687 [MainThread ] [INFO ] 62a771b0-ba3a-4895-9132-b426ffdefa34/posre.itp 2022-09-15 12:06:34,705 [MainThread ] [INFO ] Compressing topology to: 3HTB_JZ4_complex.top.zip 2022-09-15 12:06:34,705 [MainThread ] [INFO ] Ignored file 62a771b0-ba3a-4895-9132-b426ffdefa34/amber99sb-ildn.ff/forcefield.itp 2022-09-15 12:06:34,713 [MainThread ] [INFO ] Ignored file 62a771b0-ba3a-4895-9132-b426ffdefa34/amber99sb-ildn.ff/spce.itp 2022-09-15 12:06:34,713 [MainThread ] [INFO ] Ignored file 62a771b0-ba3a-4895-9132-b426ffdefa34/amber99sb-ildn.ff/ions.itp 2022-09-15 12:06:34,715 [MainThread ] [INFO ] Adding: 2022-09-15 12:06:34,715 [MainThread ] [INFO ] ['62a771b0-ba3a-4895-9132-b426ffdefa34/JZ4_posres.itp', '62a771b0-ba3a-4895-9132-b426ffdefa34/JZ4params.itp', '62a771b0-ba3a-4895-9132-b426ffdefa34/ligand.top', '62a771b0-ba3a-4895-9132-b426ffdefa34/posre.itp'] 2022-09-15 12:06:34,716 [MainThread ] [INFO ] to: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_complex.top.zip 2022-09-15 12:06:34,717 [MainThread ] [INFO ] Removed: ['62a771b0-ba3a-4895-9132-b426ffdefa34']
Out[14]:
0
Create solvent box¶
Define the unit cell for the protein-ligand complex to fill it with water molecules.
Truncated octahedron box is used for the unit cell. This box type is the one which best reflects the geometry of the solute/protein, in this case a globular protein, as it approximates a sphere. It is also convenient for the computational point of view, as it accumulates less water molecules at the corners, reducing the final number of water molecules in the system and making the simulation run faster.
A protein to box distance of 0.8 nm is used, and the protein is centered in the box.
Building Blocks used:
- Editconf from biobb_gromacs.gromacs.editconf
In [15]:
# Editconf: Create solvent box
# Import module
from biobb_gromacs.gromacs.editconf import editconf
# Create prop dict and inputs/outputs
output_editconf_gro = pdbCode+'_'+ligandCode+'_complex_editconf.gro'
prop = {
'box_type': 'octahedron',
'distance_to_molecule': 0.8
}
# Create and launch bb
editconf(input_gro_path=complexFile_H,
output_gro_path=output_editconf_gro,
properties=prop)
2022-09-15 12:06:34,754 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.editconf Version: 3.8.0
2022-09-15 12:06:34,756 [MainThread ] [INFO ] Centering molecule in the box.
2022-09-15 12:06:34,756 [MainThread ] [INFO ] Distance of the box to molecule: 0.80
2022-09-15 12:06:34,756 [MainThread ] [INFO ] Box type: octahedron
2022-09-15 12:06:34,757 [MainThread ] [INFO ] GROMACS Editconf 20222 version detected
2022-09-15 12:06:34,758 [MainThread ] [INFO ] gmx -nobackup -nocopyright editconf -f 3HTB_JZ4_H.pdb -o 3HTB_JZ4_complex_editconf.gro -d 0.8 -bt octahedron -c
2022-09-15 12:06:34,806 [MainThread ] [INFO ] Exit code 0
2022-09-15 12:06:34,807 [MainThread ] [INFO ] Note that major changes are planned in future for editconf, to improve usability and utility.
Read 2636 atoms
Volume: 1322.64 nm^3, corresponds to roughly 595100 electrons
No velocities found
system size : 4.069 4.192 5.073 (nm)
diameter : 5.849 (nm)
center : 2.241 -1.669 -0.935 (nm)
box vectors : 8.398 10.808 14.572 (nm)
box angles : 90.00 90.00 90.00 (degrees)
box volume :1322.64 (nm^3)
shift : 1.483 6.937 3.976 (nm)
new center : 3.725 5.267 3.041 (nm)
new box vectors : 7.449 7.449 7.449 (nm)
new box angles : 70.53 109.47 70.53 (degrees)
new box volume : 318.21 (nm^3)
2022-09-15 12:06:34,808 [MainThread ] [INFO ] :-) GROMACS - gmx editconf, 2022.2-conda_forge (-:
Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx
Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial
Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks
Command line:
gmx -nobackup -nocopyright editconf -f 3HTB_JZ4_H.pdb -o 3HTB_JZ4_complex_editconf.gro -d 0.8 -bt octahedron -c
GROMACS reminds you: "Your Excellency, I have no need of this hypothesis." (Pierre Laplace, to Napoleon on why his works on celestial mechanics make no mention of God.)
2022-09-15 12:06:34,808 [MainThread ] [INFO ] Removed: []
Out[15]:
0
In [16]:
# Solvate: Fill the box with water molecules
from biobb_gromacs.gromacs.solvate import solvate
# Create prop dict and inputs/outputs
output_solvate_gro = pdbCode+'_'+ligandCode+'_solvate.gro'
output_solvate_top_zip = pdbCode+'_'+ligandCode+'_solvate_top.zip'
# Create and launch bb
solvate(input_solute_gro_path=output_editconf_gro,
output_gro_path=output_solvate_gro,
input_top_zip_path=output_complex_top,
output_top_zip_path=output_solvate_top_zip)
2022-09-15 12:06:34,852 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.solvate Version: 3.8.0
2022-09-15 12:06:34,855 [MainThread ] [INFO ] Extracting: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_complex.top.zip
2022-09-15 12:06:34,855 [MainThread ] [INFO ] to:
2022-09-15 12:06:34,855 [MainThread ] [INFO ] ['d24db90e-9670-4851-b315-de8ebeb2269f/JZ4_posres.itp', 'd24db90e-9670-4851-b315-de8ebeb2269f/JZ4params.itp', 'd24db90e-9670-4851-b315-de8ebeb2269f/ligand.top', 'd24db90e-9670-4851-b315-de8ebeb2269f/posre.itp']
2022-09-15 12:06:34,856 [MainThread ] [INFO ] Unzipping:
2022-09-15 12:06:34,856 [MainThread ] [INFO ] 3HTB_JZ4_complex.top.zip
2022-09-15 12:06:34,857 [MainThread ] [INFO ] To:
2022-09-15 12:06:34,857 [MainThread ] [INFO ] d24db90e-9670-4851-b315-de8ebeb2269f/JZ4_posres.itp
2022-09-15 12:06:34,858 [MainThread ] [INFO ] d24db90e-9670-4851-b315-de8ebeb2269f/JZ4params.itp
2022-09-15 12:06:34,858 [MainThread ] [INFO ] d24db90e-9670-4851-b315-de8ebeb2269f/ligand.top
2022-09-15 12:06:34,858 [MainThread ] [INFO ] d24db90e-9670-4851-b315-de8ebeb2269f/posre.itp
2022-09-15 12:06:34,859 [MainThread ] [INFO ] GROMACS Solvate 20222 version detected
2022-09-15 12:06:34,859 [MainThread ] [INFO ] gmx -nobackup -nocopyright solvate -cp 3HTB_JZ4_complex_editconf.gro -cs spc216.gro -o 3HTB_JZ4_solvate.gro -p d24db90e-9670-4851-b315-de8ebeb2269f/ligand.top
2022-09-15 12:06:35,064 [MainThread ] [INFO ] Exit code 0
2022-09-15 12:06:35,065 [MainThread ] [INFO ]
WARNING: Masses and atomic (Van der Waals) radii will be guessed
based on residue and atom names, since they could not be
definitively assigned from the information in your input
files. These guessed numbers might deviate from the mass
and radius of the atom type. Please check the output
files if necessary. Note, that this functionality may
be removed in a future GROMACS version. Please, consider
using another file format for your input.
NOTE: From version 5.0 gmx solvate uses the Van der Waals radii
from the source below. This means the results may be different
compared to previous GROMACS versions.
++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
A. Bondi
van der Waals Volumes and Radii
J. Phys. Chem. 68 (1964) pp. 441-451
-------- -------- --- Thank You --- -------- --------
Adding line for 9531 solvent molecules with resname (SOL) to topology file (d24db90e-9670-4851-b315-de8ebeb2269f/ligand.top)
2022-09-15 12:06:35,065 [MainThread ] [INFO ] :-) GROMACS - gmx solvate, 2022.2-conda_forge (-:
Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx
Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial
Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks
Command line:
gmx -nobackup -nocopyright solvate -cp 3HTB_JZ4_complex_editconf.gro -cs spc216.gro -o 3HTB_JZ4_solvate.gro -p d24db90e-9670-4851-b315-de8ebeb2269f/ligand.top
Reading solute configuration
Reading solvent configuration
Initialising inter-atomic distances...
Generating solvent configuration
Will generate new solvent configuration of 5x4x4 boxes
Solvent box contains 36060 atoms in 12020 residues
Removed 5094 solvent atoms due to solvent-solvent overlap
Removed 2373 solvent atoms due to solute-solvent overlap
Sorting configuration
Found 1 molecule type:
SOL ( 3 atoms): 9531 residues
Generated solvent containing 28593 atoms in 9531 residues
Writing generated configuration to 3HTB_JZ4_solvate.gro
Output configuration contains 31229 atoms in 9695 residues
Volume : 318.205 (nm^3)
Density : 1000.12 (g/l)
Number of solvent molecules: 9531
Processing topology
GROMACS reminds you: "When the universe has expanded, time will contract" (Franz Ferdinand)
2022-09-15 12:06:35,066 [MainThread ] [INFO ] Compressing topology to: 3HTB_JZ4_solvate_top.zip
2022-09-15 12:06:35,067 [MainThread ] [INFO ] Ignored file d24db90e-9670-4851-b315-de8ebeb2269f/amber99sb-ildn.ff/forcefield.itp
2022-09-15 12:06:35,074 [MainThread ] [INFO ] Ignored file d24db90e-9670-4851-b315-de8ebeb2269f/amber99sb-ildn.ff/spce.itp
2022-09-15 12:06:35,074 [MainThread ] [INFO ] Ignored file d24db90e-9670-4851-b315-de8ebeb2269f/amber99sb-ildn.ff/ions.itp
2022-09-15 12:06:35,077 [MainThread ] [INFO ] Adding:
2022-09-15 12:06:35,077 [MainThread ] [INFO ] ['d24db90e-9670-4851-b315-de8ebeb2269f/JZ4_posres.itp', 'd24db90e-9670-4851-b315-de8ebeb2269f/JZ4params.itp', 'd24db90e-9670-4851-b315-de8ebeb2269f/ligand.top', 'd24db90e-9670-4851-b315-de8ebeb2269f/posre.itp']
2022-09-15 12:06:35,078 [MainThread ] [INFO ] to: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_solvate_top.zip
2022-09-15 12:06:35,079 [MainThread ] [INFO ] Removed: ['d24db90e-9670-4851-b315-de8ebeb2269f']
Out[16]:
0
Visualizing 3D structure¶
Visualizing the protein-ligand complex with the newly added solvent box using NGL
Note the octahedral box filled with water molecules surrounding the protein structure, which is centered right in the middle of the box.
In [17]:
#Show protein
view = nglview.show_structure_file(output_solvate_gro)
view.clear_representations()
view.add_representation(repr_type='cartoon', selection='protein', color='sstruc')
view.add_representation(repr_type='licorice', radius='.5', selection=ligandCode)
view.add_representation(repr_type='line', linewidth='1', selection='SOL', opacity='.3')
view._remote_call('setSize', target='Widget', args=['','600px'])
view.camera='orthographic'
view
NGLWidget()
Step 1: Creating portable binary run file for ion generation¶
In [18]:
# Grompp: Creating portable binary run file for ion generation
from biobb_gromacs.gromacs.grompp import grompp
# Create prop dict and inputs/outputs
prop = {
'mdp':{
'nsteps':'5000'
},
'simulation_type':'minimization',
'maxwarn': 1
}
output_gppion_tpr = pdbCode+'_'+ligandCode+'_complex_gppion.tpr'
# Create and launch bb
grompp(input_gro_path=output_solvate_gro,
input_top_zip_path=output_solvate_top_zip,
output_tpr_path=output_gppion_tpr,
properties=prop)
2022-09-15 12:06:35,209 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.grompp Version: 3.8.0 2022-09-15 12:06:35,212 [MainThread ] [INFO ] Extracting: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_solvate_top.zip 2022-09-15 12:06:35,213 [MainThread ] [INFO ] to: 2022-09-15 12:06:35,213 [MainThread ] [INFO ] ['d14e4218-8c9a-4bff-8cde-518183e188f1/JZ4_posres.itp', 'd14e4218-8c9a-4bff-8cde-518183e188f1/JZ4params.itp', 'd14e4218-8c9a-4bff-8cde-518183e188f1/ligand.top', 'd14e4218-8c9a-4bff-8cde-518183e188f1/posre.itp'] 2022-09-15 12:06:35,213 [MainThread ] [INFO ] Unzipping: 2022-09-15 12:06:35,214 [MainThread ] [INFO ] 3HTB_JZ4_solvate_top.zip 2022-09-15 12:06:35,214 [MainThread ] [INFO ] To: 2022-09-15 12:06:35,214 [MainThread ] [INFO ] d14e4218-8c9a-4bff-8cde-518183e188f1/JZ4_posres.itp 2022-09-15 12:06:35,215 [MainThread ] [INFO ] d14e4218-8c9a-4bff-8cde-518183e188f1/JZ4params.itp 2022-09-15 12:06:35,215 [MainThread ] [INFO ] d14e4218-8c9a-4bff-8cde-518183e188f1/ligand.top 2022-09-15 12:06:35,215 [MainThread ] [INFO ] d14e4218-8c9a-4bff-8cde-518183e188f1/posre.itp 2022-09-15 12:06:35,216 [MainThread ] [INFO ] GROMACS Grompp 20222 version detected 2022-09-15 12:06:35,216 [MainThread ] [INFO ] gmx -nobackup -nocopyright grompp -f 07b1d1dd-9c47-4592-8ecc-a02aeb0d884b/grompp.mdp -c 3HTB_JZ4_solvate.gro -r 3HTB_JZ4_solvate.gro -p d14e4218-8c9a-4bff-8cde-518183e188f1/ligand.top -o 3HTB_JZ4_complex_gppion.tpr -po mdout.mdp -maxwarn 1 2022-09-15 12:06:35,476 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:35,477 [MainThread ] [INFO ] Setting the LD random seed to -172007507 Generated 2556 of the 2556 non-bonded parameter combinations Generated 2556 of the 2556 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein_chain_A' Excluding 3 bonded neighbours molecule type 'JZ4params' Excluding 2 bonded neighbours molecule type 'SOL' ++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++ J. S. Hub, B. L. de Groot, H. Grubmueller, G. Groenhof Quantifying Artifacts in Ewald Simulations of Inhomogeneous Systems with a Net Charge J. Chem. Theory Comput. 10 (2014) pp. 381-393 -------- -------- --- Thank You --- -------- -------- Analysing residue names: There are: 163 Protein residues There are: 1 Other residues There are: 9531 Water residues Analysing Protein... Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups... The largest distance between excluded atoms is 0.419 nm Calculating fourier grid dimensions for X Y Z Using a fourier grid of 64x64x64, spacing 0.116 0.116 0.116 Estimate for the relative computational load of the PME mesh part: 0.18 This run will generate roughly 2 Mb of data 2022-09-15 12:06:35,477 [MainThread ] [INFO ] :-) GROMACS - gmx grompp, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx -nobackup -nocopyright grompp -f 07b1d1dd-9c47-4592-8ecc-a02aeb0d884b/grompp.mdp -c 3HTB_JZ4_solvate.gro -r 3HTB_JZ4_solvate.gro -p d14e4218-8c9a-4bff-8cde-518183e188f1/ligand.top -o 3HTB_JZ4_complex_gppion.tpr -po mdout.mdp -maxwarn 1 Ignoring obsolete mdp entry 'ns-type' Generating 1-4 interactions: fudge = 0.5 NOTE 1 [file ligand.top, line 24863]: System has non-zero total charge: 6.000000 Total charge should normally be an integer. See http://www.gromacs.org/Documentation/Floating_Point_Arithmetic for discussion on how close it should be to an integer. WARNING 1 [file ligand.top, line 24863]: You are using Ewald electrostatics in a system with net charge. This can lead to severe artifacts, such as ions moving into regions with low dielectric, due to the uniform background charge. We suggest to neutralize your system with counter ions, possibly in combination with a physiological salt concentration. Number of degrees of freedom in T-Coupling group rest is 65091.00 NOTE 2 [file 07b1d1dd-9c47-4592-8ecc-a02aeb0d884b/grompp.mdp]: Removing center of mass motion in the presence of position restraints might cause artifacts. When you are using position restraints to equilibrate a macro-molecule, the artifacts are usually negligible. There were 2 notes There was 1 warning GROMACS reminds you: "When the universe has expanded, time will contract" (Franz Ferdinand) 2022-09-15 12:06:35,480 [MainThread ] [INFO ] Removed: ['d14e4218-8c9a-4bff-8cde-518183e188f1', '07b1d1dd-9c47-4592-8ecc-a02aeb0d884b', 'mdout.mdp']
Out[18]:
0
Step 2: Adding ions to neutralize the system and reach a 0.05 molar concentration¶
Replace solvent molecules with ions to neutralize the system and reaching a 0.05 molar ionic concentration
In [19]:
# Genion: Adding ions to reach a 0.05 molar concentration
from biobb_gromacs.gromacs.genion import genion
# Create prop dict and inputs/outputs
prop={
'neutral':True,
'concentration':0.05
}
output_genion_gro = pdbCode+'_'+ligandCode+'_genion.gro'
output_genion_top_zip = pdbCode+'_'+ligandCode+'_genion_top.zip'
# Create and launch bb
genion(input_tpr_path=output_gppion_tpr,
output_gro_path=output_genion_gro,
input_top_zip_path=output_solvate_top_zip,
output_top_zip_path=output_genion_top_zip,
properties=prop)
2022-09-15 12:06:35,524 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.genion Version: 3.8.0 2022-09-15 12:06:35,527 [MainThread ] [INFO ] Extracting: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_solvate_top.zip 2022-09-15 12:06:35,527 [MainThread ] [INFO ] to: 2022-09-15 12:06:35,528 [MainThread ] [INFO ] ['f769f185-1a01-4d3c-89ec-75ca69ba59b5/JZ4_posres.itp', 'f769f185-1a01-4d3c-89ec-75ca69ba59b5/JZ4params.itp', 'f769f185-1a01-4d3c-89ec-75ca69ba59b5/ligand.top', 'f769f185-1a01-4d3c-89ec-75ca69ba59b5/posre.itp'] 2022-09-15 12:06:35,528 [MainThread ] [INFO ] Unzipping: 2022-09-15 12:06:35,529 [MainThread ] [INFO ] 3HTB_JZ4_solvate_top.zip 2022-09-15 12:06:35,529 [MainThread ] [INFO ] To: 2022-09-15 12:06:35,530 [MainThread ] [INFO ] f769f185-1a01-4d3c-89ec-75ca69ba59b5/JZ4_posres.itp 2022-09-15 12:06:35,530 [MainThread ] [INFO ] f769f185-1a01-4d3c-89ec-75ca69ba59b5/JZ4params.itp 2022-09-15 12:06:35,530 [MainThread ] [INFO ] f769f185-1a01-4d3c-89ec-75ca69ba59b5/ligand.top 2022-09-15 12:06:35,531 [MainThread ] [INFO ] f769f185-1a01-4d3c-89ec-75ca69ba59b5/posre.itp 2022-09-15 12:06:35,531 [MainThread ] [INFO ] To reach up 0.05 mol/litre concentration 2022-09-15 12:06:35,532 [MainThread ] [INFO ] GROMACS Genion 20222 version detected 2022-09-15 12:06:35,532 [MainThread ] [INFO ] echo "SOL" | gmx -nobackup -nocopyright genion -s 3HTB_JZ4_complex_gppion.tpr -o 3HTB_JZ4_genion.gro -p f769f185-1a01-4d3c-89ec-75ca69ba59b5/ligand.top -neutral -conc 0.05 -seed 1993 2022-09-15 12:06:35,623 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:35,625 [MainThread ] [INFO ] Will try to add 10 NA ions and 16 CL ions. Select a continuous group of solvent molecules Selected 15: 'SOL' Processing topology Replacing 26 solute molecules in topology file (f769f185-1a01-4d3c-89ec-75ca69ba59b5/ligand.top) by 10 NA and 16 CL ions. 2022-09-15 12:06:35,625 [MainThread ] [INFO ] :-) GROMACS - gmx genion, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx -nobackup -nocopyright genion -s 3HTB_JZ4_complex_gppion.tpr -o 3HTB_JZ4_genion.gro -p f769f185-1a01-4d3c-89ec-75ca69ba59b5/ligand.top -neutral -conc 0.05 -seed 1993 Reading file 3HTB_JZ4_complex_gppion.tpr, VERSION 2022.2-conda_forge (single precision) Reading file 3HTB_JZ4_complex_gppion.tpr, VERSION 2022.2-conda_forge (single precision) Group 0 ( System) has 31229 elements Group 1 ( Protein) has 2614 elements Group 2 ( Protein-H) has 1301 elements Group 3 ( C-alpha) has 163 elements Group 4 ( Backbone) has 489 elements Group 5 ( MainChain) has 653 elements Group 6 ( MainChain+Cb) has 805 elements Group 7 ( MainChain+H) has 815 elements Group 8 ( SideChain) has 1799 elements Group 9 ( SideChain-H) has 648 elements Group 10 ( Prot-Masses) has 2614 elements Group 11 ( non-Protein) has 28615 elements Group 12 ( Other) has 22 elements Group 13 ( JZ4) has 22 elements Group 14 ( Water) has 28593 elements Group 15 ( SOL) has 28593 elements Group 16 ( non-Water) has 2636 elements Select a group: Number of (3-atomic) solvent molecules: 9531 Using random seed 1993. Replacing solvent molecule 1541 (atom 7259) with NA Replacing solvent molecule 8441 (atom 27959) with NA Replacing solvent molecule 8991 (atom 29609) with NA Replacing solvent molecule 2646 (atom 10574) with NA Replacing solvent molecule 7436 (atom 24944) with NA Replacing solvent molecule 7473 (atom 25055) with NA Replacing solvent molecule 4459 (atom 16013) with NA Replacing solvent molecule 2713 (atom 10775) with NA Replacing solvent molecule 1704 (atom 7748) with NA Replacing solvent molecule 263 (atom 3425) with NA Replacing solvent molecule 9220 (atom 30296) with CL Replacing solvent molecule 1491 (atom 7109) with CL Replacing solvent molecule 4138 (atom 15050) with CL Replacing solvent molecule 1264 (atom 6428) with CL Replacing solvent molecule 1799 (atom 8033) with CL Replacing solvent molecule 5081 (atom 17879) with CL Replacing solvent molecule 8928 (atom 29420) with CL Replacing solvent molecule 1377 (atom 6767) with CL Replacing solvent molecule 1462 (atom 7022) with CL Replacing solvent molecule 6146 (atom 21074) with CL Replacing solvent molecule 2812 (atom 11072) with CL Replacing solvent molecule 6808 (atom 23060) with CL Replacing solvent molecule 8111 (atom 26969) with CL Replacing solvent molecule 2599 (atom 10433) with CL Replacing solvent molecule 1194 (atom 6218) with CL Replacing solvent molecule 3086 (atom 11894) with CL GROMACS reminds you: "When the universe has expanded, time will contract" (Franz Ferdinand) 2022-09-15 12:06:35,626 [MainThread ] [INFO ] Compressing topology to: 3HTB_JZ4_genion_top.zip 2022-09-15 12:06:35,627 [MainThread ] [INFO ] Ignored file f769f185-1a01-4d3c-89ec-75ca69ba59b5/amber99sb-ildn.ff/forcefield.itp 2022-09-15 12:06:35,634 [MainThread ] [INFO ] Ignored file f769f185-1a01-4d3c-89ec-75ca69ba59b5/amber99sb-ildn.ff/spce.itp 2022-09-15 12:06:35,635 [MainThread ] [INFO ] Ignored file f769f185-1a01-4d3c-89ec-75ca69ba59b5/amber99sb-ildn.ff/ions.itp 2022-09-15 12:06:35,638 [MainThread ] [INFO ] Adding: 2022-09-15 12:06:35,638 [MainThread ] [INFO ] ['f769f185-1a01-4d3c-89ec-75ca69ba59b5/JZ4_posres.itp', 'f769f185-1a01-4d3c-89ec-75ca69ba59b5/JZ4params.itp', 'f769f185-1a01-4d3c-89ec-75ca69ba59b5/ligand.top', 'f769f185-1a01-4d3c-89ec-75ca69ba59b5/posre.itp'] 2022-09-15 12:06:35,639 [MainThread ] [INFO ] to: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_genion_top.zip 2022-09-15 12:06:35,640 [MainThread ] [INFO ] Removed: ['f769f185-1a01-4d3c-89ec-75ca69ba59b5']
Out[19]:
0
Visualizing 3D structure¶
Visualizing the protein-ligand complex with the newly added ionic concentration using NGL
In [20]:
#Show protein
view = nglview.show_structure_file(output_genion_gro)
view.clear_representations()
view.add_representation(repr_type='cartoon', selection='protein', color='sstruc')
view.add_representation(repr_type='licorice', radius='.5', selection=ligandCode)
view.add_representation(repr_type='ball+stick', selection='NA')
view.add_representation(repr_type='ball+stick', selection='CL')
view._remote_call('setSize', target='Widget', args=['','600px'])
view.camera='orthographic'
view
NGLWidget()
Energetically minimize the system¶
Energetically minimize the protein-ligand complex till reaching a desired potential energy.
- Step 1: Creating portable binary run file for energy minimization
- Step 2: Energetically minimize the protein-ligand complex till reaching a force of 500 kJ mol-1 nm-1.
- Step 3: Checking energy minimization results. Plotting energy by time during the minimization process.
Building Blocks used:
- Grompp from biobb_gromacs.gromacs.grompp
- Mdrun from biobb_gromacs.gromacs.mdrun
- GMXEnergy from biobb_analysis.gromacs.gmx_energy
Step 1: Creating portable binary run file for energy minimization¶
Method used to run the energy minimization is a steepest descent, with a **maximum force of 500 KJ/mol*nm^2, and a minimization step size of 1fs. The maximum number of steps** to perform if the maximum force is not reached is 5,000 steps.
In [21]:
# Grompp: Creating portable binary run file for mdrun
from biobb_gromacs.gromacs.grompp import grompp
# Create prop dict and inputs/outputs
prop = {
'mdp':{
'nsteps':'5000',
'emstep': 0.01,
'emtol':'500'
},
'simulation_type':'minimization'
}
output_gppmin_tpr = pdbCode+'_'+ligandCode+'_gppmin.tpr'
# Create and launch bb
grompp(input_gro_path=output_genion_gro,
input_top_zip_path=output_genion_top_zip,
output_tpr_path=output_gppmin_tpr,
properties=prop)
2022-09-15 12:06:35,788 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.grompp Version: 3.8.0 2022-09-15 12:06:35,792 [MainThread ] [INFO ] Extracting: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks/3HTB_JZ4_genion_top.zip 2022-09-15 12:06:35,792 [MainThread ] [INFO ] to: 2022-09-15 12:06:35,793 [MainThread ] [INFO ] ['e10ba0fb-fdfc-4598-a332-8912ad996dd7/JZ4_posres.itp', 'e10ba0fb-fdfc-4598-a332-8912ad996dd7/JZ4params.itp', 'e10ba0fb-fdfc-4598-a332-8912ad996dd7/ligand.top', 'e10ba0fb-fdfc-4598-a332-8912ad996dd7/posre.itp'] 2022-09-15 12:06:35,793 [MainThread ] [INFO ] Unzipping: 2022-09-15 12:06:35,793 [MainThread ] [INFO ] 3HTB_JZ4_genion_top.zip 2022-09-15 12:06:35,794 [MainThread ] [INFO ] To: 2022-09-15 12:06:35,794 [MainThread ] [INFO ] e10ba0fb-fdfc-4598-a332-8912ad996dd7/JZ4_posres.itp 2022-09-15 12:06:35,794 [MainThread ] [INFO ] e10ba0fb-fdfc-4598-a332-8912ad996dd7/JZ4params.itp 2022-09-15 12:06:35,794 [MainThread ] [INFO ] e10ba0fb-fdfc-4598-a332-8912ad996dd7/ligand.top 2022-09-15 12:06:35,795 [MainThread ] [INFO ] e10ba0fb-fdfc-4598-a332-8912ad996dd7/posre.itp 2022-09-15 12:06:35,796 [MainThread ] [INFO ] GROMACS Grompp 20222 version detected 2022-09-15 12:06:35,796 [MainThread ] [INFO ] gmx -nobackup -nocopyright grompp -f ad6b4ef7-fa59-4922-a0b4-a61f23f39e44/grompp.mdp -c 3HTB_JZ4_genion.gro -r 3HTB_JZ4_genion.gro -p e10ba0fb-fdfc-4598-a332-8912ad996dd7/ligand.top -o 3HTB_JZ4_gppmin.tpr -po mdout.mdp -maxwarn 10 2022-09-15 12:06:36,056 [MainThread ] [INFO ] Exit code 0 2022-09-15 12:06:36,057 [MainThread ] [INFO ] Setting the LD random seed to 2004877054 Generated 2556 of the 2556 non-bonded parameter combinations Generated 2556 of the 2556 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Protein_chain_A' Excluding 3 bonded neighbours molecule type 'JZ4params' Excluding 2 bonded neighbours molecule type 'SOL' Excluding 1 bonded neighbours molecule type 'NA' Excluding 1 bonded neighbours molecule type 'CL' Analysing residue names: There are: 163 Protein residues There are: 1 Other residues There are: 9505 Water residues There are: 26 Ion residues Analysing Protein... Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups... Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups... The largest distance between excluded atoms is 0.419 nm Calculating fourier grid dimensions for X Y Z Using a fourier grid of 64x64x64, spacing 0.116 0.116 0.116 Estimate for the relative computational load of the PME mesh part: 0.18 This run will generate roughly 2 Mb of data 2022-09-15 12:06:36,058 [MainThread ] [INFO ] :-) GROMACS - gmx grompp, 2022.2-conda_forge (-: Executable: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial/bin.SSE2/gmx Data prefix: /Users/gbayarri/opt/anaconda3/envs/biobb_Protein-Complex_MDsetup_tutorial Working dir: /Users/gbayarri/biobb_wf_protein-complex_md_setup/biobb_wf_protein-complex_md_setup/notebooks Command line: gmx -nobackup -nocopyright grompp -f ad6b4ef7-fa59-4922-a0b4-a61f23f39e44/grompp.mdp -c 3HTB_JZ4_genion.gro -r 3HTB_JZ4_genion.gro -p e10ba0fb-fdfc-4598-a332-8912ad996dd7/ligand.top -o 3HTB_JZ4_gppmin.tpr -po mdout.mdp -maxwarn 10 Ignoring obsolete mdp entry 'ns-type' Generating 1-4 interactions: fudge = 0.5 Number of degrees of freedom in T-Coupling group rest is 65013.00 NOTE 1 [file ad6b4ef7-fa59-4922-a0b4-a61f23f39e44/grompp.mdp]: Removing center of mass motion in the presence of position restraints might cause artifacts. When you are using position restraints to equilibrate a macro-molecule, the artifacts are usually negligible. There was 1 note GROMACS reminds you: "I am rarely happier than when spending an entire day programming my computer to perform automatically a task that it would otherwise take me a good ten seconds to do by hand." (Douglas Adams) 2022-09-15 12:06:36,060 [MainThread ] [INFO ] Removed: ['e10ba0fb-fdfc-4598-a332-8912ad996dd7', 'ad6b4ef7-fa59-4922-a0b4-a61f23f39e44', 'mdout.mdp']
Out[21]:
0
Step 2: Running Energy Minimization¶
Running energy minimization using the tpr file generated in the previous step.
In [ ]:
# Mdrun: Running minimization
from biobb_gromacs.gromacs.mdrun import mdrun
# Create prop dict and inputs/outputs
output_min_trr = pdbCode+'_'+ligandCode+'_min.trr'
output_min_gro = pdbCode+'_'+ligandCode+'_min.gro'
output_min_edr = pdbCode+'_'+ligandCode+'_min.edr'
output_min_log = pdbCode+'_'+ligandCode+'_min.log'
# Create and launch bb
mdrun(input_tpr_path=output_gppmin_tpr,
output_trr_path=output_min_trr,
output_gro_path=output_min_gro,
output_edr_path=output_min_edr,
output_log_path=output_min_log)
2022-09-15 12:06:36,099 [MainThread ] [INFO ] Executing biobb_gromacs.gromacs.mdrun Version: 3.8.0 2022-09-15 12:06:36,100 [MainThread ] [INFO ] GROMACS Mdrun 20222 version detected 2022-09-15 12:06:36,100 [MainThread ] [INFO ] gmx -nobackup -nocopyright mdrun -s 3HTB_JZ4_gppmin.tpr -o 3HTB_JZ4_min.trr -c 3HTB_JZ4_min.gro -e 3HTB_JZ4_min.edr -g 3HTB_JZ4_min.log
Step 3: Checking Energy Minimization results¶
Checking energy minimization results. Plotting potential energy by time during the minimization process.
In [ ]:
# GMXEnergy: Getting system energy by time
from biobb_analysis.gromacs.gmx_energy import gmx_energy
# Create prop dict and inputs/outputs
output_min_ene_xvg = pdbCode+'_'+ligandCode+'_min_ene.xvg'
prop = {
'terms': ["Potential"]
}
# Create and launch bb
gmx_energy(input_energy_path=output_min_edr,
output_xvg_path=output_min_ene_xvg,
properties=prop)
In [ ]:
import plotly
import plotly.graph_objs as go
#Read data from file and filter energy values higher than 1000 Kj/mol^-1
with open(output_min_ene_xvg,'r') as energy_file:
x,y = map(
list,
zip(*[
(float(line.split()[0]),float(line.split()[1]))
for line in energy_file
if not line.startswith(("#","@"))
if float(line.split()[1]) < 1000
])
)
plotly.offline.init_notebook_mode(connected=True)
fig = ({
"data": [go.Scatter(x=x, y=y)],
"layout": go.Layout(title="Energy Minimization",
xaxis=dict(title = "Energy Minimization Step"),
yaxis=dict(title = "Potential Energy KJ/mol-1")
)
})
plotly.offline.iplot(fig)
Equilibrate the system (NVT)¶
Equilibrate the protein-ligand complex system in NVT ensemble (constant Number of particles, Volume and Temperature). To avoid temperature coupling problems, a new "system" group will be created including the protein + the ligand to be assigned to a single thermostatting group.
- Step 1: Creating an index file with a new group including the protein-ligand complex.
- Step 2: Creating portable binary run file for system equilibration
- Step 3: Equilibrate the protein-ligand complex with NVT ensemble.
- Step 4: Checking NVT Equilibration results. Plotting system temperature by time during the NVT equilibration process.
Building Blocks used:
- MakeNdx from biobb_gromacs.gromacs.make_ndx
- Grompp from biobb_gromacs.gromacs.grompp
- Mdrun from biobb_gromacs.gromacs.mdrun
- GMXEnergy from biobb_analysis.gromacs.gmx_energy
Step 1: Creating an index file with a new group including the protein-ligand complex¶
In [ ]:
# MakeNdx: Creating index file with a new group (protein-ligand complex)
from biobb_gromacs.gromacs.make_ndx import make_ndx
# Create prop dict and inputs/outputs
output_complex_ndx = pdbCode+'_'+ligandCode+'_index.ndx'
prop = {
'selection': "\"Protein\"|\"Other\""
}
# Create and launch bb
make_ndx(input_structure_path=output_min_gro,
output_ndx_path=output_complex_ndx,
properties=prop)
Step 2: Creating portable binary run file for system equilibration (NVT)¶
Note that for the purposes of temperature coupling, the protein-ligand complex (Protein_Other) is considered as a single entity.
In [ ]:
# Grompp: Creating portable binary run file for NVT System Equilibration
from biobb_gromacs.gromacs.grompp import grompp
# Create prop dict and inputs/outputs
output_gppnvt_tpr = pdbCode+'_'+ligandCode+'gppnvt.tpr'
prop = {
'mdp':{
'nsteps':'5000',
'tc-grps': 'Protein_Other Water_and_ions',
'Define': '-DPOSRES -D' + posresifdef
},
'simulation_type':'nvt'
}
# Create and launch bb
grompp(input_gro_path=output_min_gro,
input_top_zip_path=output_genion_top_zip,
input_ndx_path=output_complex_ndx,
output_tpr_path=output_gppnvt_tpr,
properties=prop)
Step 3: Running NVT equilibration¶
In [ ]:
# Mdrun: Running NVT System Equilibration
from biobb_gromacs.gromacs.mdrun import mdrun
# Create prop dict and inputs/outputs
output_nvt_trr = pdbCode+'_'+ligandCode+'_nvt.trr'
output_nvt_gro = pdbCode+'_'+ligandCode+'_nvt.gro'
output_nvt_edr = pdbCode+'_'+ligandCode+'_nvt.edr'
output_nvt_log = pdbCode+'_'+ligandCode+'_nvt.log'
output_nvt_cpt = pdbCode+'_'+ligandCode+'_nvt.cpt'
# Create and launch bb
mdrun(input_tpr_path=output_gppnvt_tpr,
output_trr_path=output_nvt_trr,
output_gro_path=output_nvt_gro,
output_edr_path=output_nvt_edr,
output_log_path=output_nvt_log,
output_cpt_path=output_nvt_cpt)
Step 4: Checking NVT Equilibration results¶
Checking NVT Equilibration results. Plotting system temperature by time during the NVT equilibration process.
In [ ]:
# GMXEnergy: Getting system temperature by time during NVT Equilibration
from biobb_analysis.gromacs.gmx_energy import gmx_energy
# Create prop dict and inputs/outputs
output_nvt_temp_xvg = pdbCode+'_'+ligandCode+'_nvt_temp.xvg'
prop = {
'terms': ["Temperature"]
}
# Create and launch bb
gmx_energy(input_energy_path=output_nvt_edr,
output_xvg_path=output_nvt_temp_xvg,
properties=prop)
In [ ]:
import plotly
import plotly.graph_objs as go
# Read temperature data from file
with open(output_nvt_temp_xvg,'r') as temperature_file:
x,y = map(
list,
zip(*[
(float(line.split()[0]),float(line.split()[1]))
for line in temperature_file
if not line.startswith(("#","@"))
])
)
plotly.offline.init_notebook_mode(connected=True)
fig = ({
"data": [go.Scatter(x=x, y=y)],
"layout": go.Layout(title="Temperature during NVT Equilibration",
xaxis=dict(title = "Time (ps)"),
yaxis=dict(title = "Temperature (K)")
)
})
plotly.offline.iplot(fig)
Equilibrate the system (NPT)¶
Equilibrate the protein-ligand complex system in NPT ensemble (constant Number of particles, Pressure and Temperature) .
- Step 1: Creating portable binary run file for system equilibration
- Step 2: Equilibrate the protein-ligand complex with NPT ensemble.
- Step 3: Checking NPT Equilibration results. Plotting system pressure and density by time during the NPT equilibration process.
Building Blocks used:
- Grompp from biobb_gromacs.gromacs.grompp
- Mdrun from biobb_gromacs.gromacs.mdrun
- GMXEnergy from biobb_analysis.gromacs.gmx_energy
Step 1: Creating portable binary run file for system equilibration (NPT)¶
In [ ]:
# Grompp: Creating portable binary run file for (NPT) System Equilibration
from biobb_gromacs.gromacs.grompp import grompp
# Create prop dict and inputs/outputs
output_gppnpt_tpr = pdbCode+'_'+ligandCode+'_gppnpt.tpr'
prop = {
'mdp':{
'type': 'npt',
'nsteps':'5000',
'tc-grps': 'Protein_Other Water_and_ions',
'Define': '-DPOSRES -D' + posresifdef
},
'simulation_type':'npt'
}
# Create and launch bb
grompp(input_gro_path=output_nvt_gro,
input_top_zip_path=output_genion_top_zip,
input_ndx_path=output_complex_ndx,
output_tpr_path=output_gppnpt_tpr,
input_cpt_path=output_nvt_cpt,
properties=prop)
Step 2: Running NPT equilibration¶
In [ ]:
# Mdrun: Running NPT System Equilibration
from biobb_gromacs.gromacs.mdrun import mdrun
# Create prop dict and inputs/outputs
output_npt_trr = pdbCode+'_'+ligandCode+'_npt.trr'
output_npt_gro = pdbCode+'_'+ligandCode+'_npt.gro'
output_npt_edr = pdbCode+'_'+ligandCode+'_npt.edr'
output_npt_log = pdbCode+'_'+ligandCode+'_npt.log'
output_npt_cpt = pdbCode+'_'+ligandCode+'_npt.cpt'
# Create and launch bb
mdrun(input_tpr_path=output_gppnpt_tpr,
output_trr_path=output_npt_trr,
output_gro_path=output_npt_gro,
output_edr_path=output_npt_edr,
output_log_path=output_npt_log,
output_cpt_path=output_npt_cpt)
Step 3: Checking NPT Equilibration results¶
Checking NPT Equilibration results. Plotting system pressure and density by time during the NPT equilibration process.
In [ ]:
# GMXEnergy: Getting system pressure and density by time during NPT Equilibration
from biobb_analysis.gromacs.gmx_energy import gmx_energy
# Create prop dict and inputs/outputs
output_npt_pd_xvg = pdbCode+'_'+ligandCode+'_npt_PD.xvg'
prop = {
'terms': ["Pressure","Density"]
}
# Create and launch bb
gmx_energy(input_energy_path=output_npt_edr,
output_xvg_path=output_npt_pd_xvg,
properties=prop)
In [ ]:
import plotly
from plotly import subplots
import plotly.graph_objs as go
# Read pressure and density data from file
with open(output_npt_pd_xvg,'r') as pd_file:
x,y,z = map(
list,
zip(*[
(float(line.split()[0]),float(line.split()[1]),float(line.split()[2]))
for line in pd_file
if not line.startswith(("#","@"))
])
)
plotly.offline.init_notebook_mode(connected=True)
trace1 = go.Scatter(
x=x,y=y
)
trace2 = go.Scatter(
x=x,y=z
)
fig = subplots.make_subplots(rows=1, cols=2, print_grid=False)
fig.append_trace(trace1, 1, 1)
fig.append_trace(trace2, 1, 2)
fig['layout']['xaxis1'].update(title='Time (ps)')
fig['layout']['xaxis2'].update(title='Time (ps)')
fig['layout']['yaxis1'].update(title='Pressure (bar)')
fig['layout']['yaxis2'].update(title='Density (Kg*m^-3)')
fig['layout'].update(title='Pressure and Density during NPT Equilibration')
fig['layout'].update(showlegend=False)
plotly.offline.iplot(fig)
Free Molecular Dynamics Simulation¶
Upon completion of the two equilibration phases (NVT and NPT), the system is now well-equilibrated at the desired temperature and pressure. The position restraints can now be released. The last step of the protein-ligand complex MD setup is a short, free MD simulation, to ensure the robustness of the system.
- Step 1: Creating portable binary run file to run a free MD simulation.
- Step 2: Run short MD simulation of the protein-ligand complex.
- Step 3: Checking results for the final step of the setup process, the free MD run. Plotting Root Mean Square deviation (RMSd) and Radius of Gyration (Rgyr) by time during the free MD run step.
Building Blocks used:
- Grompp from biobb_gromacs.gromacs.grompp
- Mdrun from biobb_gromacs.gromacs.mdrun
- GMXRms from biobb_analysis.gromacs.gmx_rms
- GMXRgyr from biobb_analysis.gromacs.gmx_rgyr
Step 1: Creating portable binary run file to run a free MD simulation¶
In [ ]:
# Grompp: Creating portable binary run file for mdrun
from biobb_gromacs.gromacs.grompp import grompp
# Create prop dict and inputs/outputs
prop = {
'mdp':{
#'nsteps':'500000' # 1 ns (500,000 steps x 2fs per step)
#'nsteps':'5000' # 10 ps (5,000 steps x 2fs per step)
'nsteps':'25000' # 50 ps (25,000 steps x 2fs per step)
},
'simulation_type':'free'
}
output_gppmd_tpr = pdbCode+'_'+ligandCode + '_gppmd.tpr'
# Create and launch bb
grompp(input_gro_path=output_npt_gro,
input_top_zip_path=output_genion_top_zip,
output_tpr_path=output_gppmd_tpr,
input_cpt_path=output_npt_cpt,
properties=prop)
Step 2: Running short free MD simulation¶
In [ ]:
# Mdrun: Running free dynamics
from biobb_gromacs.gromacs.mdrun import mdrun
# Create prop dict and inputs/outputs
output_md_trr = pdbCode+'_'+ligandCode+'_md.trr'
output_md_gro = pdbCode+'_'+ligandCode+'_md.gro'
output_md_edr = pdbCode+'_'+ligandCode+'_md.edr'
output_md_log = pdbCode+'_'+ligandCode+'_md.log'
output_md_cpt = pdbCode+'_'+ligandCode+'_md.cpt'
# Create and launch bb
mdrun(input_tpr_path=output_gppmd_tpr,
output_trr_path=output_md_trr,
output_gro_path=output_md_gro,
output_edr_path=output_md_edr,
output_log_path=output_md_log,
output_cpt_path=output_md_cpt)
Step 3: Checking free MD simulation results¶
Checking results for the final step of the setup process, the free MD run. Plotting Root Mean Square deviation (RMSd) and Radius of Gyration (Rgyr) by time during the free MD run step. RMSd against the experimental structure (input structure of the pipeline) and against the minimized and equilibrated structure (output structure of the NPT equilibration step).
In [ ]:
# GMXRms: Computing Root Mean Square deviation to analyse structural stability
# RMSd against minimized and equilibrated snapshot (backbone atoms)
from biobb_analysis.gromacs.gmx_rms import gmx_rms
# Create prop dict and inputs/outputs
output_rms_first = pdbCode+'_'+ligandCode+'_rms_first.xvg'
prop = {
'selection': 'Backbone'
}
# Create and launch bb
gmx_rms(input_structure_path=output_gppmd_tpr,
input_traj_path=output_md_trr,
output_xvg_path=output_rms_first,
properties=prop)
In [ ]:
# GMXRms: Computing Root Mean Square deviation to analyse structural stability
# RMSd against experimental structure (backbone atoms)
from biobb_analysis.gromacs.gmx_rms import gmx_rms
# Create prop dict and inputs/outputs
output_rms_exp = pdbCode+'_'+ligandCode+'_rms_exp.xvg'
prop = {
'selection': 'Backbone'
}
# Create and launch bb
gmx_rms(input_structure_path=output_gppmin_tpr,
input_traj_path=output_md_trr,
output_xvg_path=output_rms_exp,
properties=prop)
In [ ]:
import plotly
import plotly.graph_objs as go
# Read RMS vs first snapshot data from file
with open(output_rms_first,'r') as rms_first_file:
x,y = map(
list,
zip(*[
(float(line.split()[0]),float(line.split()[1]))
for line in rms_first_file
if not line.startswith(("#","@"))
])
)
# Read RMS vs experimental structure data from file
with open(output_rms_exp,'r') as rms_exp_file:
x2,y2 = map(
list,
zip(*[
(float(line.split()[0]),float(line.split()[1]))
for line in rms_exp_file
if not line.startswith(("#","@"))
])
)
trace1 = go.Scatter(
x = x,
y = y,
name = 'RMSd vs first'
)
trace2 = go.Scatter(
x = x,
y = y2,
name = 'RMSd vs exp'
)
data = [trace1, trace2]
plotly.offline.init_notebook_mode(connected=True)
fig = ({
"data": data,
"layout": go.Layout(title="RMSd during free MD Simulation",
xaxis=dict(title = "Time (ps)"),
yaxis=dict(title = "RMSd (nm)")
)
})
plotly.offline.iplot(fig)
In [ ]:
# GMXRgyr: Computing Radius of Gyration to measure the protein compactness during the free MD simulation
from biobb_analysis.gromacs.gmx_rgyr import gmx_rgyr
# Create prop dict and inputs/outputs
output_rgyr = pdbCode+'_'+ligandCode+'_rgyr.xvg'
prop = {
'selection': 'Backbone'
}
# Create and launch bb
gmx_rgyr(input_structure_path=output_gppmin_tpr,
input_traj_path=output_md_trr,
output_xvg_path=output_rgyr,
properties=prop)
In [ ]:
import plotly
import plotly.graph_objs as go
# Read Rgyr data from file
with open(output_rgyr,'r') as rgyr_file:
x,y = map(
list,
zip(*[
(float(line.split()[0]),float(line.split()[1]))
for line in rgyr_file
if not line.startswith(("#","@"))
])
)
plotly.offline.init_notebook_mode(connected=True)
fig = ({
"data": [go.Scatter(x=x, y=y)],
"layout": go.Layout(title="Radius of Gyration",
xaxis=dict(title = "Time (ps)"),
yaxis=dict(title = "Rgyr (nm)")
)
})
plotly.offline.iplot(fig)
Post-processing and Visualizing resulting 3D trajectory¶
Post-processing and Visualizing the protein-ligand complex system MD setup resulting trajectory using NGL
- Step 1: Imaging the resulting trajectory, stripping out water molecules and ions and correcting periodicity issues.
- Step 2: Generating a dry structure, removing water molecules and ions from the final snapshot of the MD setup pipeline.
- Step 3: Visualizing the imaged trajectory using the dry structure as a topology.
Building Blocks used:
- GMXImage from biobb_analysis.gromacs.gmx_image
- GMXTrjConvStr from biobb_analysis.gromacs.gmx_trjconv_str
Step 1: Imaging the resulting trajectory.¶
Stripping out water molecules and ions and correcting periodicity issues
In [ ]:
# GMXImage: "Imaging" the resulting trajectory
# Removing water molecules and ions from the resulting structure
from biobb_analysis.gromacs.gmx_image import gmx_image
# Create prop dict and inputs/outputs
output_imaged_traj = pdbCode+'_imaged_traj.trr'
prop = {
'center_selection': 'Protein_Other',
'output_selection': 'Protein_Other',
'pbc' : 'mol',
'center' : True
}
# Create and launch bb
gmx_image(input_traj_path=output_md_trr,
input_top_path=output_gppmd_tpr,
input_index_path=output_complex_ndx,
output_traj_path=output_imaged_traj,
properties=prop)
Step 2: Generating the output dry structure.¶
Removing water molecules and ions from the resulting structure
In [ ]:
# GMXTrjConvStr: Converting and/or manipulating a structure
# Removing water molecules and ions from the resulting structure
# The "dry" structure will be used as a topology to visualize
# the "imaged dry" trajectory generated in the previous step.
from biobb_analysis.gromacs.gmx_trjconv_str import gmx_trjconv_str
# Create prop dict and inputs/outputs
output_dry_gro = pdbCode+'_md_dry.gro'
prop = {
'selection': 'Protein_Other'
}
# Create and launch bb
gmx_trjconv_str(input_structure_path=output_md_gro,
input_top_path=output_gppmd_tpr,
input_index_path=output_complex_ndx,
output_str_path=output_dry_gro,
properties=prop)
Step 3: Visualizing the generated dehydrated trajectory.¶
Using the imaged trajectory (output of the Post-processing step 1) with the dry structure (output of the Post-processing step 2) as a topology.
In [ ]:
# Show trajectory
view = nglview.show_simpletraj(nglview.SimpletrajTrajectory(output_imaged_traj, output_dry_gro), gui=True)
view
Output files¶
Important Output files generated:
- {{output_md_gro}}: Final structure (snapshot) of the MD setup protocol.
- {{output_md_trr}}: Final trajectory of the MD setup protocol.
- {{output_md_cpt}}: Final checkpoint file, with information about the state of the simulation. It can be used to restart or continue a MD simulation.
- {{output_gppmd_tpr}}: Final tpr file, GROMACS portable binary run input file. This file contains the starting structure of the MD setup free MD simulation step, together with the molecular topology and all the simulation parameters. It can be used to extend the simulation.
- {{output_genion_top_zip}}: Final topology of the MD system. It is a compressed zip file including a topology file (.top) and a set of auxiliar include topology files (.itp).
Analysis (MD setup check) output files generated:
- {{output_rms_first}}: Root Mean Square deviation (RMSd) against minimized and equilibrated structure of the final free MD run step.
- {{output_rms_exp}}: Root Mean Square deviation (RMSd) against experimental structure of the final free MD run step.
- {{output_rgyr}}: Radius of Gyration of the final free MD run step of the setup pipeline.
Questions & Comments¶
Questions, issues, suggestions and comments are really welcome!
GitHub issues:
BioExcel forum: