DFT-based Steric Parameters
Allows a user to compute steric parameters from chemical structures.
Calculate Sterimol parameters1 (L, Bmin, Bmax), %Buried Volume2, Sterimol2Vec and Vol2Vec parameters
- Compute requested steric parameters from molecular structure files with input options:
-sor--sterimol- Sterimol Parameters (L, Bmin, Bmax)-bor--vbur- Percent Buried Volume-sor--sterimolAND--scan [rmin:rmax:interval]- Sterimol2Vec Parameters-bor--vburAND--scan [rmin:rmax:interval]- Vol2Vec Parameters
-r- Adjust radius of percent buried volume measurements (default 3.5 Angstrom)--dp [n]- Number of decimal places in the printed results (default 2)--cutoff [Å or auto]- Ignore atoms farther than this from atom1.autokeeps exactly the atoms that can occupy the buried-volume sphere, so %V_Bur is unchanged while the grid stays small for large systems (clusters, proteins). With a cutoff, the molecular volume column reports only the kept atoms (MolVol_cut) and Sterimol L is bounded by the cutoff.--quiet- Suppress all printed output (results are still available from the Python object)- Exclude atoms from steric measurement with
--exclude [atom indices]option (no spaces, separated by commas) - Sterimol parameters can be computed using the classic (Verloop) definition from van der Waals radii, or using a three-dimensional grid (default is classic).
- Change measurement type with
--measure ['classic' or 'grid']. Grid-based measurement is automatically used for scans and density surfaces. - Grid point spacing can be adjusted (default spacing is 0.05 Angstrom), adjust with
--grid [# in Angstrom] --pos- Only measure Sterimol parameters in the positive direction (from atom1 toward atom2) in grid mode--atom3 [idx]- Align a third atom to the positive x direction to fully define the molecular orientation--norot- Skip the alignment step for structures that are already aligned along the z-axis
- Change measurement type with
- Three sets of VDW radii are available:
--radii bondi(default) - Bondi radii--radii charry-tkatchenko- Charry-Tkatchenko free-atom radii derived from dipole polarizability--radii cpk- the CPK radii of Verloop's original Sterimol program (and of the earlier patonlab/sterimol code). These depend on the bonding environment: atoms are typed from DFT-D3 style coordination numbers (sp3 vs. aromatic carbon, single- vs. double-bonded oxygen, tetrahedral vs. planar nitrogen, ...), so no bond perception is needed and all input formats work. Elements without a CPK type (metals, B, Si, ...) fall back to Bondi. With this set, classic Sterimol values reproduce the original Fortran program to 0.01 Å; the assigned types are available asmetadata["cpk_type"]on the Python object--scalevdw [factor]- Scale the chosen radii (default 1.0)
- Steric parameters can be measured from electron density .cube files generated by Gaussian (see Gaussian cubegen for information on how to generate these)
- The
--surface densitycommand (default vdw) with a .cube input file will measure sterics from density values read in from the file. - Density values read from the cube file greater than a default cutoff of 0.0016 determine if a molecule is occupying that point in space, this can be changed with
--isoval [number]
- The
- Multi-structure
.xyzand.sdffiles are processed structure by structure, with each row labelled by the structure's name (comment line / SDF title) -tor--tensor- Return a 3D binary occupancy tensor of the aligned molecule (requires--atom1,--atom2and--atom3); add--saveto write it to a.npyfile and--pymolto visualize the voxels--noH- exclude hydrogen atoms from steric measurements--nometals- exclude metal atoms from steric measurements--sambvca- SambVca 2.1 mode (Bondi radii scaled by 1.17, H atoms excluded)--cone- Tolman cone angle and metal-centroid distance of a ligand in a metal complex, with the ligand's Sterimol parameters from the metal (see Metal complexes)--boltzmann- Boltzmann-weighted parameters over conformer ensembles: multi-record SDF/xyz files, or one QM output file per conformer (energies read by cclib);--energy-windowdrops high-energy conformers (see Conformer ensembles)
- Compute graph-based steric contributions in layers spanning outward from a reference functional group with the following input options:
--2d- Toggle 2D measurements on--fg- Specify an atom or functional group to use as a reference as a SMILES string--maxpath- The number of layers to measure. A connectivity matrix is used to compute the shortest path to each atom from the reference functional group.--2d-type- The type of steric contributions to use. Options include Crippen molar refractivities or McGowan volume
See CHANGELOG.md for what changed in each release.
- Python 3.10 or greater
- Non-standard dependencies will be installed along with DBSTEP, but include numpy, scipy, and cclib.
uv add dbstep
- Install using conda
conda install -c conda-forge dbstep - Or using pip
pip install dbstep
After installation, the dbstep command is available directly on the command line, or run as a module with python -m dbstep.
git clone https://github.com/patonlab/DBSTEP.git
cd DBSTEP
uv sync --extra dev
Please reference the DOI of our Zenodo repository with:
Luchini, G.; Patterson, T.; Paton, R. S. DBSTEP: DFT Based Steric Parameters. 2022, DOI: 10.5281/zenodo.4702097
DBSTEP reads .xyz (single or multi-structure), .sdf/.mol (V2000, single or multi-structure), .pdb/.ent (Protein Data Bank, single or multi-MODEL) and Gaussian .com/.gjf input files natively, and Gaussian 16 cube files containing volumetric density information. Quantum chemistry output files are parsed with the cclib module; for the list of supported programs see their documentation here. When used from a Python script, DBSTEP can also read coordinates from RDKit mol objects that carry a 3D conformer.
With a PDB file you can pick a residue and atoms by name instead of file indices. The measurement is centred on the chosen atom and only the atoms within reach of the buried-volume sphere are kept (--cutoff auto is switched on automatically), so a whole protein runs in seconds.
>>>dbstep 1a8o.pdb --residue A:186 --vbur --sterimol --nowater
File Atom1 Atom2 R/Å MolVol_cut %V_Bur %S_Bur Bmin Bmax L
-----------------------------------------------------------------------------------------------------------
1a8o.pdb A:186 THR CA CB 3.50 533.54 65.95 0.00 5.54 7.12 7.19
-----------------------------------------------------------------------------------------------------------
-
--residue A:45selects chain A residue 45 (append an insertion code:A:45A;45alone works when only one chain has it; several residues separated by commas form one selection) -
--atom CAnames atom1 within the residue (default CA);--atom2and--atom3accept atom names too (default atom2: CB, falling back to HA or N for glycine) -
--nowaterdrops water molecules;--nohetdrops other hetero groups such as ligands and ions (modified residues with a peptide backbone, e.g. MSE, stay);--chain Akeeps a single chain -
--exclude-selflets the residue occupy no volume, so only its environment is measured (a pocket size);--self-onlykeeps only the residue, which equals extracting it to its own file -
Sterimol L measured from CA through the surroundings is a distance to the nearest steric wall along the CA→CB direction, not a substituent length, so do not compare it with Verloop values
-
Crystal structures usually lack hydrogens; results differ from a protonated model.
--sambvca(heavy atoms only, Bondi × 1.17) is a consistent choice for raw PDB files -
--decomposesplits %V_Bur between the residues whose atoms fill the sphere (a grid point covered by several residues is shared equally, so the contributions add up to the total); printed after the table and, with--csv out.csv, written toout_contributions.csv. In Python the object carriescontributions, a dict from residue label to percent -
--residue allruns every polymer residue in turn (waters and ligands are skipped, modified residues such as MSE are included; combine with--chain) -
--csv results.csvwrites one row per file, frame, residue and radius with the columnsfile, frame, structure, residue, atom1, atom2, radius, mol_vol, percent_vbur, percent_sbur, bmin, bmax, L, cone_angle, metal_centroid, population, path; this works for any input, not only PDB files
A multi-record .sdf from a conformer search (AQME, CREST, RDKit) is run record by record; --boltzmann adds a population to every conformer and a final boltzmann row with the population-weighted %V_Bur, %S_Bur, L, Bmin and Bmax. Energies are read from an SDF data field such as <Energy> (also E, G, dG and similar; name another one with --boltzmann TAG), or from a floating-point number in the comment line of a multi-frame .xyz (CREST style). --energy-units (kcal, kJ, hartree, eV; default kcal/mol) and --temperature (default 298.15 K) control the weights, and a structure without an energy is an error rather than silently dropped.
>>>dbstep ether_conformers.sdf --atom1 3 --atom2 2 --sterimol --vbur --boltzmann --csv ether.csv
File Atom1 Atom2 R/Å Mol_Vol %V_Bur %S_Bur Bmin Bmax L
-----------------------------------------------------------------------------------------------------------------------
ether 44 3 2 3.50 79.69 39.98 0.00 1.98 3.32 4.09
ether 12 3 2 3.50 79.74 39.74 0.00 1.99 4.28 4.08
ether 6 3 2 3.50 79.81 39.63 0.00 2.00 4.23 4.14
ether_conformers.sdf boltzmann 3 2 3.50 79.69 39.96 0.00 1.98 3.40 4.09
-----------------------------------------------------------------------------------------------------------------------
Boltzmann populations at 298.15 K (kcal/mol): ether 44 0.923, ether 12 0.071, ether 6 0.006
From Python: runs = db.all_frames("ether_conformers.sdf", atom1=3, atom2=2, sterimol=True, volume=True) then dbstep.ensemble.boltzmann_average(runs, temperature=298.15, units="kcal") returns the summary rows and sets population, energy (kcal/mol), energy_rel and in_window on each run.
One QM output per conformer. When every input file holds a single structure (Gaussian, ORCA, ... outputs read by cclib, plain or gzipped, or single-structure xyz files with an energy in the comment line), --boltzmann pools all the files into one ensemble and adds an ensemble boltzmann row. Energies are taken from the output itself: the Gibbs free energy after a frequency calculation, otherwise the last SCF energy (--boltzmann E or --boltzmann G to choose), always in hartree whatever --energy-units says. This is the workflow of the earlier wSterimol code (its conformer generation is covered by AQME or CREST, and its example reproduces: wL 6.33, wB1 1.79, wB5 3.75 with --radii cpk at 298 K). --energy-window 3.0 leaves out conformers more than 3 kcal/mol above the lowest; they stay in the table and the CSV with population 0.
>>>dbstep pentane_*.out --sterimol --atom1 1 --atom2 3 --radii cpk --boltzmann --temperature 298 --energy-window 1.0
File Atom1 Atom2 Bmin Bmax L
-------------------------------------------------------------------------
pentane_1.out 1 3 1.93 4.50 4.63
pentane_10.out 1 3 1.90 4.10 5.89
pentane_18.out 1 3 1.68 2.74 6.99
...
ensemble boltzmann 1 3 1.80 3.74 5.94
-------------------------------------------------------------------------
Boltzmann populations at 298.00 K (hartree): pentane_1.out 0.115, pentane_10.out 0.123, pentane_18.out 0.277, ...
Energies of ensemble: SCF energy read from the output files
Energy window 1.00 kcal/mol: 7 of 9 structures weighted
In Python, run the files one by one and pass the list to boltzmann_average(runs, window=1.0, label="pentane").
Multi-frame .xyz, multi-record .sdf and multi-MODEL .pdb files are trajectories: every frame is measured in turn (the neighbourhood crop is recomputed per frame, so frames cost the same as single structures). --frames start:stop:stride selects frames with Python slice rules on the 0-based index, e.g. --frames 0:1000:10 or --frames ::5, and --csv collects the time series:
dbstep md_frames.pdb --residue A:45 --vbur --nowater --frames ::10 --csv vbur_A45.csv
From Python, db.all_frames("md_frames.pdb", frames="::10", residue="A:45", volume=True, nowater=True) returns one object per frame; the frame and structure columns of results identify it. Binary trajectory formats (DCD, XTC, ...) are not read yet; convert to multi-MODEL PDB or multi-frame XYZ first.
From Python: db.dbstep("1a8o.pdb", residue="A:186", atom="CA", atom2="CB", volume=True, nowater=True, exclude_self=True); the object also exposes atoms, coords and metadata for the atoms that were actually measured, and results, the list of row dictionaries that --csv writes. db.all_residues("1a8o.pdb", volume=True, nowater=True) returns one such object per residue.
To execute the program:
-
Run from the command line with:
dbstep file --atom1 a1idx --atom2 a2idx(orpython -m dbstep ...) -
Run in a Python program by importing:
import dbstep.Dbstep as db(example below)
import dbstep.Dbstep as db
# Create DBSTEP object
mol = db.dbstep(file, atom1=atom1, atom2=atom2, sterimol=True)
# Grab Sterimol Parameters
L = mol.L
Bmin = mol.Bmin
Bmax = mol.BmaxDBSTEP currently takes a coordinate file (see information on appropriate file types above) along with reference atoms and other input options for steric measurement. Sterimol parameters are measured and output to the user using the --sterimol argument, volume parameters can be requested with the --vbur option.
Atoms are specified by referring to the index of an atom in a coordinate file, (ex: "2", referencing the second atom in the file, with indexing starting at 1).
For Sterimol parameters, two atoms need to be specified using the arguments --atom1 [atom1idx] and --atom2 [atom2idx]. The L parameter is measured starting from the specified atom1 coordinates, extending through the atom1-atom2 axis until the end of the molecule is reached. The Bmin and Bmax molecular width parameters are measured on the axis perpendicular to L.
For buried volume parameters, only the --atom1 [atom] argument is necessary to specify.
If no atoms are specified, the first two atoms in the file will be used as reference.
--cone measures one ligand of a metal complex the way the earlier patonlab/sterimol code did for half-sandwich complexes: the Tolman cone angle and the metal-to-centroid distance, together with the ligand's Sterimol parameters measured from the metal along the metal-centroid axis (every other ligand is excluded, the metal itself occupies no volume). --atom1 is the metal and --atom2 the ring atoms (comma separated) or the donor atom of the ligand; with neither given DBSTEP takes the single metal, the largest ring bound to it or, failing that, the nearest donor atom. The ligand is the bonded fragment containing the axis atoms once the metal is removed. Each ligand atom subtends the half angle alpha + asin(r/d) at the metal; the cone angle is twice the mean, over the ligand's sectors (the wedges around each ring atom, or the substituent branches of a donor atom), of the largest half angle in the sector, Tolman's construction for unsymmetrical ligands. Tolman's tabulated values came from CPK models, so --radii cpk is the closest match; the cone angle depends on the chosen radii like every other DBSTEP quantity.
>>>dbstep RhCpMe5Cl2PMe3.log --cone --radii cpk
File Atom1 Atom2 Bmin Bmax L Cone/° M-Cent/Å
-----------------------------------------------------------------------------------------
RhCpMe5Cl2PMe3.log 1 3,4,5,24,25 3.91 4.30 4.05 173.97 1.83
-----------------------------------------------------------------------------------------
Cone angle of RhCpMe5Cl2PMe3.log: apex atom 1, ligand of 25 atoms (axis atoms 3,4,5,24,25), sector half angles 89.3, 87.6, 81.0, 89.3, 87.6
>>>dbstep RhCpMe5Cl2PMe3.log --cone --atom2 17 --radii cpk # the PMe3 ligand through its P atom
-b adds the buried volume of the ligand around the metal, --boltzmann averages cone angles over conformers, and --csv writes cone_angle and metal_centroid columns. In Python the object carries cone_angle, metal_centroid, cone_sectors (half angle per sector) and ligand_atoms. For dimers or several metals, choose the centre with --atom1.
The same measurements are available inside PyMOL, on the structures you have open, with atoms picked by PyMOL selections and the results drawn on the molecule (this replaces the visual side of wSterimol). Install DBSTEP into the Python that PyMOL uses (pip install dbstep; the open-source pymol-open-source wheel on PyPI works too), then in PyMOL or in your .pymolrc:
import dbstep.pymol_plugin
| Command | What it does |
|---|---|
dbstep_sterimol atom1, atom2 [, atom3, radii=bondi, measure=classic, selection=, name=sterimol] |
L, Bmin and Bmax from atom1 along the atom1-atom2 axis, drawn as a blue L axis and green/red Bmin/Bmax circles in the object's own frame |
dbstep_vbur atom1 [, radius=3.5, radii=bondi, selection=, decompose=1] |
%V_bur in a sphere around atom1; with decompose every residue's B-factor is set to its contribution and the contributing residues are coloured white to red |
dbstep_cone metal [, ligand, radii=cpk] |
Tolman cone angle, metal-centroid distance and the ligand's Sterimol parameters, drawn as a translucent cone plus the distance |
dbstep_vdw object [, radii=bondi, scale=1.0] |
Translucent van der Waals copy object_vdw with DBSTEP's radii (Bondi, Charry-Tkatchenko or CPK atom types) |
dbstep_ensemble files, atom1, atom2 [, radii, temperature, window, vbur=1] |
Boltzmann-weighted parameters over conformers (QM outputs or a multi-structure file), loaded into PyMOL with their populations as titles and a transparency that fades out minor conformers |
dbstep_conformers folder [, pattern=*.pdb] |
Load every structure of a folder into one group |
dbstep_style |
White background and light settings for figures |
selection restricts the atoms that are measured; atoms named as atom1/atom2 outside it are kept as zero-radius ghosts, so dbstep_vbur /1a8o//A/186/CA, 3.5, polymer and not resi 186 measures the pocket around Thr186 without waters and without the residue itself. Every command returns the DBSTEP run object, so from the PyMOL Python prompt the full results (results, contributions, cone_sectors, ...) are available. From Python outside PyMOL the same in-memory route is dbstep.Dbstep.from_coords(atoms, coords, ...).
A notebook covering the protein, trajectory and conformer-ensemble workflows is at examples/proteins_and_conformers.ipynb.
Examples for obtaining Sterimol, Sterimol2Vec, Percent Buried Volume and Vol2Vec parameter sets are shown below (all example files found in dbstep/data/ directory).
-
Sterimol Parameters for Ethane
Obtain the Sterimol parameters for an ethane molecule along the C2-C5 bond on the command line:
>>>python -m dbstep dbstep/data/Et.xyz --sterimol --atom1 2 --atom2 5
File Atom1 Atom2 Bmin Bmax L
-------------------------------------------------------
Et.xyz 2 5 1.99 2.13 3.24
-------------------------------------------------------
A visualization of these parameters can be shown in PyMOL using the two output files created by DBSTEP, showing the L parameter in blue, Bmin parameter in green and Bmax parameter in red.
-
Sterimol2Vec Parameters for Ph
The
--scanargument is formatted asrmin:rmax:intervalwhere rmin is the distance from the center along the L axis to start measurements, rmax dictates when to stop measurements, and interval is the frequency of measurements. In this case the length of the molecule (~6A) is measured in 1.0A intervals
>>>python -m dbstep dbstep/data/Ph.xyz --sterimol --atom1 1 --atom2 2 --scan 0.0:6.0:1.0
File Atom1 Atom2 Bmin Bmax L
-------------------------------------------------------
Ph.xyz 1 2 1.65 3.16 1.00
Ph.xyz 1 2 1.65 3.16 2.00
Ph.xyz 1 2 1.65 3.16 3.00
Ph.xyz 1 2 1.65 3.16 4.00
Ph.xyz 1 2 1.65 3.16 5.00
Ph.xyz 1 2 1.65 3.11 5.95
Ph.xyz 1 2 1.15 1.17 5.95
L parameter is 5.95 Ang
-------------------------------------------------------
Displayed in PyMOL, each new Bmin and Bmax axis is added along the L axis.
-
Percent Buried Volume
%Vb is measured by constructing a sphere (typically with a 3.5A radius) around the center atom and measuring how much of the sphere is occupied by the molecule. Output will include the sphere radius, percent buried volume (%V_Bur) and percent buried shell volume (%S_Bur) (zero in all cases unless a scan is being done simultaneously).
>>>python -m dbstep dbstep/data/1Nap.xyz --atom1 2 --vbur
File Atom R/Å Mol_Vol %V_Bur %S_Bur
---------------------------------------------------------
1Nap.xyz 2 3.50 118.65 41.77 0.00
---------------------------------------------------------
For percent buried volume, the PyMOL script will overlay an appropriate sized sphere where measurement took place.
-
Vol2Vec Parameters
When invoking the --vbur and --scan parameters simultaneously, vol2vec parameters can be obtained. In this case, a scan is performed using spheres with radii from 2.0A to 4.0A in 0.5A increments.
>>>python -m dbstep dbstep/data/CHiPr2.xyz --atom1 1 -b --scan 2.0:4.0:0.5
File Atom R/Å Mol_Vol %V_Bur %S_Bur
-----------------------------------------------------------
CHiPr2.xyz 1 2.00 116.50 58.27 49.52
CHiPr2.xyz 1 2.50 116.50 53.53 46.22
CHiPr2.xyz 1 3.00 116.50 48.78 38.14
CHiPr2.xyz 1 3.50 116.50 43.37 29.17
CHiPr2.xyz 1 4.00 116.50 36.73 16.82
-----------------------------------------------------------
-
2D Additive sterics
To calculate 2d graph-based additive sterics, the arguments --2d --fg --maxpath and --2d-type can be used. An input file listing SMILES strings of desired molecule measurements is necessary for calculation. The --fg argument specifies a SMILES string that is common in all provided SMILES inputs to use as a reference point for layer 0. A connectivity matrix will then be used to find atoms 1, 2, 3... N bonds away where N is the max path length specified with the --maxpath argument. One of two types of measurements will be summed at each layer, either Crippen molar refractivities or McGowan volumes, computed for each atom. This can be changed with the --2d-type argument.
>>>python -m dbstep dbstep/data/smiles.txt --2d --fg "C(O)=O" --maxpath 5 --2d-type mcgowan
where smiles.txt looks like:
CC(O)=O
CCC(O)=O
CCCC(O)=O
CCCCC(O)=O
CC(C)C(O)=O
CCC(C)C(O)=O
The output will then be written to the file "smiles_2d_output.csv" in the format:
| 0_mcgowan | 1_mcgowan | 2_mcgowan | 3_mcgowan | 4_mcgowan | Structure |
|---|---|---|---|---|---|
| 6.51 | 19.52 | 0 | 0 | 0 | CC(=O)O |
| 6.51 | 14.09 | 19.52 | 0 | 0 | CCC(=O)O |
| 6.51 | 14.09 | 14.09 | 19.52 | 0 | CCCC(=O)O |
| 6.51 | 14.09 | 14.09 | 14.09 | 19.52 | CCCCC(=O)O |
| 6.51 | 8.66 | 39.04 | 0 | 0 | CC(C)C(=O)O |
| 6.51 | 8.66 | 33.61 | 19.52 | 0 | CCC(C)C(=O)O |
This work is developed by Guilian Luchini, Toby Patterson and Robert Paton and is supported by the NSF Center for Computer-Assisted Synthesis, grant number CHE-1925607