To detect steric clashes in a molecular structure, you can use a combination of MDAnalysis for atom selection and distance calculations. A common approach involves calculating the distance between all atom pairs and comparing the sum of their Van der Waals (VdW) radii against the actual distance, adjusted by an allowance and a cutoff threshold.
For detecting clashes, cutoff values of 0.4-1.0 Å and allowance values of 0.2-0.6 Å are generally considered reasonable.
Note: The implementation below requires MDAnalysis and MDAnalysis.lib.distances to be available in your environment.
import MDAnalysis as mda
import MDAnalysis.topology.tables as tables
import MDAnalysis.lib.distances as distances
# Load structure
structure = "1UBQ.pdb"
ubiquitin = mda.Universe(structure)
# Define a function to check for clashes
def check_for_clashes(atoms, cutoff=0.6, allowance=0.4, vdwradii=None):
clashes = set()
positions = atoms.positions
if not vdwradii:
vdwradii = tables.vdwradii
for a in atoms:
if not a.type in vdwradii.keys():
print(f"{a.type} has no defined vdw radius")
return clashes
# Calculate all-to-all distance matrix
dist = distances.self_distance_array(positions)
N = len(positions)
# Generate pairs for the upper triangle
pairs = []
for i in range(N):
for j in range(i + 1, N):
pairs.append((i, j))
for x, (d, (i, j)) in enumerate(zip(dist, pairs)):
a1, a2 = atoms[i], atoms[j]
r1, r2 = vdwradii[a1.type], vdwradii[a2.type]
overlap = (r1 + r2) - dist[x] - allowance
if overlap >= cutoff:
clashes.add(frozenset([i+1, j+1]))
return clashes
# Example usage with specific residues
gln49 = ubiquitin.select_atoms('resid 49')
arg42 = ubiquitin.select_atoms('resid 42')
clash_group = arg42 + gln49
clashes = check_for_clashes(clash_group, cutoff=0.0, allowance=0.0)
