Skip to content

Commit 83963d4

Browse files
bobbypatonbobbypaton
committed
add Charry-Tkatchenko VDW radii, --radii & --measure CLI options, fix grid sterimol bugs
- Add Charry-Tkatchenko free-atom VDW radii (JCTC 2024) to constants.py - Add --radii CLI option to select between bondi and charry-tkatchenko radii sets - Add --measure CLI option to select classic (Verloop) or grid-based sterimol - Force Bondi radii when --sambvca is used - Add Atom2 column to output table when sterimol is active - Format sterimol-only output as a table matching the combined vbur+sterimol style - Fix grid sterimol L: recompute grid bounds after molecule rotation so the grid covers the full rotated molecule extent - Fix density grid sterimol: rotate occupied grid points to match molecular orientation so L, Bmin, Bmax are measured along the correct axes - Extract get_rotation_angles() in calculator.py for reuse - Add analysis/ directory with radii comparison script and 3-panel figure
1 parent e248cb3 commit 83963d4

9 files changed

Lines changed: 1141 additions & 58 deletions

File tree

‎.gitignore‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,3 +31,5 @@ htmlcov/
3131

3232
# Claude Code
3333
.claude/
34+
ct4c00784_si_001.pdf
35+
TODO

‎analysis/README.md‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
# Analysis
2+
3+
Comparison scripts and figures for DBSTEP steric parameter calculations.
4+
5+
## Radii Comparison: Bondi vs Charry-Tkatchenko
6+
7+
**Script:** `compare_radii.py`
8+
9+
Generates a 3-panel figure comparing the two VDW radii sets available in DBSTEP:
10+
11+
- **(a)** Scatter plot of Bondi vs Charry-Tkatchenko radii for all 54 shared elements. The Charry-Tkatchenko radii (free-atom, derived from dipole polarizability) are systematically larger than Bondi radii for most elements, particularly for transition metals (e.g. Cu, Zn, Pd, Pt).
12+
- **(b)** Molecular volumes computed with both radii sets for all 22 molecules in `dbstep/data/`.
13+
- **(c)** Percent buried volumes (%V_Bur) at R = 3.5 A for the same molecules.
14+
15+
**References:**
16+
- Bondi: *J. Phys. Chem.* **1964**, *68*, 441; Mantina et al. *J. Phys. Chem. A* **2009**, *113*, 5806
17+
- Charry-Tkatchenko: *J. Chem. Theory Comput.* **2024**, *20*, 7844-7855 (R_vdW^free[alpha] from Table S1, SI)
18+
19+
### Running
20+
21+
```bash
22+
uv run python analysis/compare_radii.py
23+
```
24+
25+
Outputs `analysis/radii_comparison.png`.
26+
27+
![Radii Comparison](radii_comparison.png)

‎analysis/compare_radii.py‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
"""
2+
Compare Bondi and Charry-Tkatchenko VDW radii sets.
3+
4+
Generates a 3-panel figure:
5+
(a) Bondi vs Charry-Tkatchenko radii for all shared elements
6+
(b) Molecular volumes with both radii sets for molecules in data/
7+
(c) Buried volumes with both radii sets for molecules in data/
8+
9+
Usage:
10+
python analysis/compare_radii.py
11+
"""
12+
13+
import os, sys
14+
import numpy as np
15+
import matplotlib.pyplot as plt
16+
17+
# ensure the package is importable when running from the repo root
18+
sys.path.insert(0, os.path.join(os.path.dirname(__file__), ".."))
19+
20+
from dbstep.Dbstep import dbstep
21+
from dbstep.constants import bondi, charry_tkatchenko, periodic_table
22+
23+
24+
def get_shared_elements():
25+
"""Return elements present in both radii dicts, ordered by atomic number."""
26+
return [el for el in periodic_table
27+
if el in bondi and el in charry_tkatchenko and el not in ("Bq", "")]
28+
29+
30+
def compute_volumes(data_dir):
31+
"""Compute Mol_Vol and %V_Bur with both radii sets for all xyz files in data_dir."""
32+
files = sorted(f for f in os.listdir(data_dir) if f.endswith(".xyz"))
33+
names, bondi_mvol, ct_mvol, bondi_vbur, ct_vbur = [], [], [], [], []
34+
35+
for f in files:
36+
path = os.path.join(data_dir, f)
37+
label = os.path.splitext(f)[0]
38+
39+
r_bondi = dbstep(path, volume=True, radii="bondi", quiet=True)
40+
r_ct = dbstep(path, volume=True, radii="charry-tkatchenko", quiet=True)
41+
42+
names.append(label)
43+
bondi_mvol.append(r_bondi.occ_vol)
44+
ct_mvol.append(r_ct.occ_vol)
45+
bondi_vbur.append(r_bondi.bur_vol)
46+
ct_vbur.append(r_ct.bur_vol)
47+
48+
return names, bondi_mvol, ct_mvol, bondi_vbur, ct_vbur
49+
50+
51+
def main():
52+
data_dir = os.path.join(os.path.dirname(__file__), "..", "dbstep", "data")
53+
54+
# --- Panel (a): radii comparison ---
55+
elements = get_shared_elements()
56+
b_radii = [bondi[el] for el in elements]
57+
ct_radii = [charry_tkatchenko[el] for el in elements]
58+
59+
# --- Panels (b) & (c): volume comparison ---
60+
print("Computing volumes for molecules in data/ ...")
61+
names, bondi_mvol, ct_mvol, bondi_vbur, ct_vbur = compute_volumes(data_dir)
62+
63+
# --- Plot ---
64+
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
65+
66+
# Panel (a): radii scatter
67+
ax = axes[0]
68+
lim = [0.8, max(max(b_radii), max(ct_radii)) + 0.2]
69+
ax.plot(lim, lim, "k--", lw=0.8, alpha=0.5)
70+
ax.scatter(b_radii, ct_radii, s=30, c="steelblue", edgecolors="k", linewidths=0.5, zorder=3)
71+
# label notable outliers
72+
for i, el in enumerate(elements):
73+
if abs(b_radii[i] - ct_radii[i]) > 0.4:
74+
ax.annotate(el, (b_radii[i], ct_radii[i]), fontsize=7, ha="left",
75+
xytext=(4, 4), textcoords="offset points")
76+
ax.set_xlabel("Bondi radius (Å)")
77+
ax.set_ylabel("Charry-Tkatchenko radius (Å)")
78+
ax.set_title("(a) VDW Radii Comparison")
79+
ax.set_xlim(lim)
80+
ax.set_ylim(lim)
81+
ax.set_aspect("equal")
82+
83+
# Panel (b): molecular volume parity
84+
ax = axes[1]
85+
mvol_lim = [0, max(max(bondi_mvol), max(ct_mvol)) * 1.1]
86+
ax.plot(mvol_lim, mvol_lim, "k--", lw=0.8, alpha=0.5)
87+
ax.scatter(bondi_mvol, ct_mvol, s=30, c="steelblue", edgecolors="k", linewidths=0.5, zorder=3)
88+
for i, name in enumerate(names):
89+
ax.annotate(name, (bondi_mvol[i], ct_mvol[i]), fontsize=6, ha="left",
90+
xytext=(4, 4), textcoords="offset points")
91+
ax.set_xlabel("Mol. Vol. Bondi (ų)")
92+
ax.set_ylabel("Mol. Vol. Charry-Tkatchenko (ų)")
93+
r2_mvol = np.corrcoef(bondi_mvol, ct_mvol)[0, 1] ** 2
94+
ax.set_title(f"(b) Molecular Volume ($R^2$ = {r2_mvol:.4f})")
95+
ax.set_xlim(mvol_lim)
96+
ax.set_ylim(mvol_lim)
97+
ax.set_aspect("equal")
98+
99+
# Panel (c): buried volume parity
100+
ax = axes[2]
101+
vbur_lim = [0, max(max(bondi_vbur), max(ct_vbur)) * 1.1]
102+
ax.plot(vbur_lim, vbur_lim, "k--", lw=0.8, alpha=0.5)
103+
ax.scatter(bondi_vbur, ct_vbur, s=30, c="steelblue", edgecolors="k", linewidths=0.5, zorder=3)
104+
for i, name in enumerate(names):
105+
ax.annotate(name, (bondi_vbur[i], ct_vbur[i]), fontsize=6, ha="left",
106+
xytext=(4, 4), textcoords="offset points")
107+
ax.set_xlabel("%V_Bur Bondi")
108+
ax.set_ylabel("%V_Bur Charry-Tkatchenko")
109+
r2_vbur = np.corrcoef(bondi_vbur, ct_vbur)[0, 1] ** 2
110+
ax.set_title(f"(c) Buried Volume ($R^2$ = {r2_vbur:.4f})")
111+
ax.set_xlim(vbur_lim)
112+
ax.set_ylim(vbur_lim)
113+
ax.set_aspect("equal")
114+
115+
plt.tight_layout()
116+
out_path = os.path.join(os.path.dirname(__file__), "radii_comparison.png")
117+
fig.savefig(out_path, dpi=200, bbox_inches="tight")
118+
print(f"Figure saved to {out_path}")
119+
plt.close()
120+
121+
122+
if __name__ == "__main__":
123+
main()

‎analysis/radii_comparison.png‎

195 KB
Loading

‎dbstep/Dbstep.py‎

Lines changed: 44 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
from optparse import OptionParser
88

99
from dbstep import sterics, parse_data, calculator, writer
10-
from dbstep.constants import periodic_table, bondi, metals
10+
from dbstep.constants import periodic_table, bondi, charry_tkatchenko, metals
1111

1212
class dbstep:
1313
"""
@@ -37,8 +37,9 @@ def __init__(self, *args, **kwargs):
3737
self.options = kwargs["options"]
3838
else:
3939
self.options = set_options(kwargs)
40-
# SambVca mode: scale VDW radii by 1.17 and exclude H atoms
40+
# SambVca mode: Bondi radii scaled by 1.17, exclude H atoms
4141
if hasattr(self.options, 'sambvca') and self.options.sambvca:
42+
self.options.radii = "bondi"
4243
self.options.SCALE_VDW = 1.17
4344
self.options.noH = True
4445

@@ -58,9 +59,6 @@ def __init__(self, *args, **kwargs):
5859
# flag volume if buried shell requested
5960
if options.vshell:
6061
options.volume = True
61-
# if measuring volume, need to measure from grid
62-
if options.volume:
63-
options.measure = "grid"
6462

6563
origin = np.array([0, 0, 0])
6664
self._get_spec_atoms(options)
@@ -73,6 +71,10 @@ def __init__(self, *args, **kwargs):
7371
# Rotate molecule to align atom1-atom2 bond along Z-axis
7472
self._orient_molecule(mol, options)
7573

74+
# Recompute grid bounds after rotation so the grid covers the full rotated molecule
75+
if options.sterimol and options.surface == "vdw":
76+
[x_min, x_max, y_min, y_max, z_min, z_max, _] = sterics.max_dim(mol.CARTESIANS, mol.RADII, options)
77+
7678
# Parse scan range
7779
r_min, r_max, r_intervals, strip_width = self._parse_scan(options)
7880

@@ -108,12 +110,13 @@ def _assign_surface(self, mol, file, options, origin):
108110
x_min = x_max = y_min = y_max = z_min = z_max = 0.0
109111

110112
if options.surface == "vdw":
111-
# generate Bondi radii from atom types (default 2.0 for elements not in the Bondi table)
113+
# Select radii set based on options
114+
radii_dict = charry_tkatchenko if options.radii == "charry-tkatchenko" else bondi
112115
for atom in mol.ATOMTYPES:
113-
if atom not in periodic_table and atom not in bondi:
116+
if atom not in periodic_table and atom not in radii_dict:
114117
print("\n UNABLE TO GENERATE VDW RADII FOR ATOM: ", atom)
115118
exit()
116-
mol.RADII = [bondi.get(atom, 2.0) for atom in mol.ATOMTYPES]
119+
mol.RADII = [radii_dict.get(atom, 2.0) for atom in mol.ATOMTYPES]
117120
mol.RADII = np.array(mol.RADII) * options.SCALE_VDW
118121

119122
# Translate molecule to place atom1 at the origin
@@ -156,11 +159,15 @@ def _assign_surface(self, mol, file, options, origin):
156159
return x_min, x_max, y_min, y_max, z_min, z_max
157160

158161
def _orient_molecule(self, mol, options):
159-
"""Rotate molecule to align atom1-atom2 bond along the Z-axis."""
162+
"""Rotate molecule to align atom1-atom2 bond along the Z-axis.
163+
Stores rotation angles in options.rotation for later use."""
164+
options.rotation = None
160165
if not options.sterimol:
161166
return
162167
point = calculator.point_vec(mol.CARTESIANS, options.spec_atom_2)
163168
if len(mol.CARTESIANS) > 1 and not options.norot:
169+
# Compute rotation angles
170+
options.rotation = calculator.get_rotation_angles(mol.CARTESIANS, options.spec_atom_1, point, options.atom3)
164171
if options.surface == "vdw":
165172
mol.CARTESIANS = calculator.rotate_mol(mol.CARTESIANS, options.spec_atom_1, point, options.verbose, options.atom3)
166173
elif options.surface == "density":
@@ -196,14 +203,14 @@ def _build_grid(self, mol, name, options, origin, x_min, x_max, y_min, y_max, z_
196203
y_vals = np.linspace(y_min, y_max, int(1 + round((y_max - y_min) / options.grid)))
197204
z_vals = np.linspace(z_min, z_max, int(1 + round((z_max - z_min) / options.grid)))
198205

199-
if options.measure == "grid":
200-
if options.volume and not options.sterimol:
201-
# Fast path: skip full grid construction for volume-only calculations
206+
if options.volume or options.measure == "grid":
207+
if options.volume and not (options.sterimol and options.measure == "grid"):
208+
# Fast path: skip full grid construction when grid sterimol not needed
202209
occ_grid, occ_vol = sterics.occupied_direct(mol.CARTESIANS, mol.RADII, origin, x_vals, y_vals, z_vals, options)
203210
point_tree = None
204211
grid_axes = (x_vals, y_vals, z_vals)
205212
else:
206-
# Standard path: full grid needed for sterimol
213+
# Standard path: full grid needed for grid-based sterimol
207214
grid = np.array(np.meshgrid(x_vals, y_vals, z_vals)).T.reshape(-1, 3)
208215
occ_grid, point_tree, occ_vol = sterics.occupied(grid, mol.CARTESIANS, mol.RADII, origin, options)
209216

@@ -216,6 +223,9 @@ def _build_grid(self, mol, name, options, origin, x_min, x_max, y_min, y_max, z_
216223
# Use 'ij' indexing so grid order matches cube file density order (x-slow, y-mid, z-fast)
217224
grid = np.array(np.meshgrid(x_vals, y_vals, z_vals, indexing='ij')).reshape(3, -1).T
218225
occ_grid, occ_vol, point_tree = sterics.occupied_dens(grid, mol.DENSITY, options)
226+
# Rotate occupied grid to match the molecular orientation (bond along Z)
227+
if options.rotation is not None:
228+
occ_grid = calculator.apply_rotation(occ_grid, options.rotation)
219229
if options.volume:
220230
grid, point_tree = sterics.resize_grid(x_max, y_max, z_max, x_min, y_min, z_min, options, mol)
221231

@@ -227,7 +237,9 @@ def _print_column_header(self, options):
227237
return
228238
fw = dbstep._file_col_width
229239
if options.volume and options.sterimol:
230-
header = " {:>{fw}} {:>6} {:>6} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format("File", "Atom", "R/Å", "Mol_Vol", "%V_Bur", "%S_Bur", "Bmin", "Bmax", "L", fw=fw)
240+
header = " {:>{fw}} {:>6} {:>6} {:>6} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10}".format("File", "Atom1", "Atom2", "R/Å", "Mol_Vol", "%V_Bur", "%S_Bur", "Bmin", "Bmax", "L", fw=fw)
241+
elif options.sterimol:
242+
header = " {:>{fw}} {:>6} {:>6} {:>10} {:>10} {:>10}".format("File", "Atom1", "Atom2", "Bmin", "Bmax", "L", fw=fw)
231243
elif options.volume:
232244
header = " {:>{fw}} {:>6} {:>6} {:>10} {:>10} {:>10}".format("File", "Atom", "R/Å", "Mol_Vol", "%V_Bur", "%S_Bur", fw=fw)
233245
else:
@@ -267,12 +279,11 @@ def _compute(self, mol, file, options, origin, occ_grid, occ_vol, point_tree, gr
267279
if options.sterimol:
268280
if options.measure == "grid":
269281
L, Bmax, Bmin, cyl = sterics.get_cube_sterimol(occ_grid, rad, options.grid, strip_width, options.pos)
270-
elif options.measure == "classic":
271-
if options.surface == "vdw":
272-
L, Bmax, Bmin, cyl = sterics.get_classic_sterimol(mol.CARTESIANS, mol.RADII, mol.ATOMTYPES)
273-
elif options.surface == "density":
274-
print(" Can't use classic Sterimol with the isodensity surface. Either use VDW radii (--surface vdw) or use grid Sterimol (--sterimol grid)")
275-
exit()
282+
elif options.surface == "vdw":
283+
L, Bmax, Bmin, cyl = sterics.get_classic_sterimol(mol.CARTESIANS, mol.RADII, mol.ATOMTYPES)
284+
else:
285+
print(" Can't use classic Sterimol with the isodensity surface. Use --measure grid or --surface vdw")
286+
exit()
276287
Bmin_list.append(Bmin)
277288
Bmax_list.append(Bmax)
278289
if options.pymol:
@@ -283,19 +294,17 @@ def _compute(self, mol, file, options, origin, occ_grid, occ_vol, point_tree, gr
283294
if options.pymol:
284295
spheres.append(" SPHERE, 0.000, 0.000, 0.000, {:5.3f},".format(rad))
285296
if not options.quiet:
286-
print(" {:>{fw}} {:>6} {:6.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f}".format(fname, options.spec_atom_1, rad, occ_vol, bur_vol, bur_shell, Bmin, Bmax, L, fw=fw))
297+
atom2_str = ",".join(str(a) for a in options.spec_atom_2)
298+
print(" {:>{fw}} {:>6} {:>6} {:6.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f} {:10.2f}".format(fname, options.spec_atom_1, atom2_str, rad, occ_vol, bur_vol, bur_shell, Bmin, Bmax, L, fw=fw))
287299
elif options.volume:
288300
if options.pymol:
289301
spheres.append(" SPHERE, 0.000, 0.000, 0.000, {:5.3f},".format(rad))
290302
if not options.quiet:
291303
print(" {:>{fw}} {:>6} {:6.2f} {:10.2f} {:10.2f} {:10.2f}".format(fname, options.spec_atom_1, rad, occ_vol, bur_vol, bur_shell, fw=fw))
292304
elif options.sterimol:
293-
if not options.scan:
294-
if not options.quiet:
295-
print(" {} / Bmin: {:5.2f} / Bmax: {:5.2f} / L: {:5.2f}".format(file, Bmin, Bmax, L))
296-
else:
297-
if not options.quiet:
298-
print(" {} / R: {:5.2f} / Bmin: {:5.2f} / Bmax: {:5.2f} ".format(file, rad, Bmin, Bmax))
305+
if not options.quiet:
306+
atom2_str = ",".join(str(a) for a in options.spec_atom_2)
307+
print(" {:>{fw}} {:>6} {:>6} {:10.2f} {:10.2f} {:10.2f}".format(fname, options.spec_atom_1, atom2_str, Bmin, Bmax, L, fw=fw))
299308

300309
# Store results on self
301310
if options.measure == "grid":
@@ -399,6 +408,7 @@ def set_options(kwargs):
399408
"vshell": ["vshell", False],
400409
"pymol": ["pymol", False],
401410
"quiet": ["quiet", False],
411+
"radii": ["radii", "bondi"],
402412
"sambvca": ["sambvca", False],
403413
"gridsize": ["gridsize", False],
404414
"measure": ["measure", "classic"],
@@ -446,11 +456,13 @@ def main():
446456
parser.add_option("--norot", dest="norot", action="store_true", help="Do not rotate the molecule (use if structures have been pre-aligned)", default=False)
447457
parser.add_option("--pos", dest="pos", action="store_true", help="Measure Sterimol parameters in positive direction (from atom1 toward atom2)", default=False)
448458
parser.add_option("--quiet", dest="quiet", action="store_true", help="Suppress all print output", default=False)
459+
parser.add_option("--radii", dest="radii", action="store", choices=["bondi", "charry-tkatchenko"], help="VDW radii set: bondi or charry-tkatchenko (default: bondi)", default="bondi")
449460
parser.add_option("-r", dest="radius", action="store", help="Radius of sphere in Angstrom (default: 3.5)", default=3.5, type=float, metavar="radius")
450461
parser.add_option("--sambvca", dest="sambvca", action="store_true", help="Use SambVca 2.1 defaults: scale VDW radii by 1.17 and exclude H atoms", default=False)
451462
parser.add_option("--scalevdw", dest="SCALE_VDW", action="store", help="Scaling factor for VDW radii (default: 1.0)", type=float, default=1.0, metavar="SCALE_VDW")
452463
parser.add_option("--scan", dest="scan", action="store", help="Scan over a range of radii, format: rmin:rmax:interval", default=False, metavar="scan")
453464
parser.add_option("-s", "--sterimol", dest="sterimol", action="store_true", help="Compute Sterimol parameters (L, Bmin, Bmax)", default=False)
465+
parser.add_option("--measure", dest="measure", action="store", choices=["classic", "grid"], help="Sterimol method: classic (Verloop, default) or grid-based", default="classic", metavar="measure")
454466
parser.add_option("--surface", dest="surface", action="store", choices=["vdw", "density"], help="Surface type: Bondi VDW radii or density cube file (default: vdw)", default="vdw", metavar="surface")
455467
parser.add_option("-b", "--vbur", dest="volume", action="store_true", help="Calculate buried volume of input molecule", default=False)
456468
parser.add_option("-v", "--verbose", dest="verbose", action="store_true", help="Print verbose output", default=False)
@@ -459,11 +471,9 @@ def main():
459471
parser.add_option("--vshell", dest="vshell", action="store", help="Calculate buried volume of hollow sphere with given shell width; use -r to set radius", default=False, type=float, metavar="width")
460472
(options, args) = parser.parse_args()
461473

462-
# Sterimol defaults to classic; volume forces grid mode internally
463-
options.measure = "classic"
464-
465-
# SambVca mode: scale VDW radii by 1.17 and exclude H atoms
474+
# SambVca mode: Bondi radii scaled by 1.17, exclude H atoms
466475
if options.sambvca:
476+
options.radii = "bondi"
467477
options.SCALE_VDW = 1.17
468478
options.noH = True
469479

@@ -502,8 +512,9 @@ def main():
502512
if options.sterimol:
503513
print(" Sterimol parameters will be generated using {} mode".format("grid-based" if options.measure == "grid" else "classic"))
504514
if options.surface == "vdw":
505-
print(" Using a Cartesian grid-spacing of {:5.4f} Angstrom".format(options.grid))
506-
print(" Bondi atomic radii will be scaled by {}".format(options.SCALE_VDW))
515+
print(" Using a Cartesian grid-spacing of {:5.4f} Angstrom".format(options.grid))
516+
radii_label = "Charry-Tkatchenko" if options.radii == "charry-tkatchenko" else "Bondi"
517+
print(" {} atomic radii will be scaled by {}".format(radii_label, options.SCALE_VDW))
507518
print(" Hydrogen atoms are {}\n".format("excluded" if options.noH else "included"))
508519
else:
509520
print(" Using {} isodensity surface with cutoff value of {:5.4f} au".format(options.surface, options.isoval))

0 commit comments

Comments
 (0)