Structural Bioinformatics

How to Calculate Redocking RMSD for AutoDock Vina: Atom Mapping, Symmetry, and Common Errors

Calculate redocking RMSD in the receptor frame, check atom mapping and symmetry, and avoid misleading ligand alignment with a tested RDKit example.

This guide is for researchers calculating reference-pose recovery in their own docking workflow. It provides a runnable, deliberately strict RDKit example and a check against BioChemIntelli's retained 1IEP case. It is not a comparison of different compounds, a binding-affinity calculation, or a substitute for validating a screening protocol.

Decide which RMSD question you are asking

Three calculations can produce a number called RMSD while answering different questions.

Calculation What stays fixed? Useful question
Distance between reported docking modes Their shared receptor frame How different are the returned poses?
Ligand-to-reference redocking RMSD The receptor frame after any justified receptor alignment Did the search recover the reference placement and conformation?
Independently superposed ligand RMSD Only the ligand's internal shape is compared after fitting How similar are two conformations after removing relative position and orientation?

The third calculation can be useful, but it is not the second. A ligand translated out of its pocket could still superpose perfectly on its reference.

RDKit distinguishes CalcRMS, which evaluates coordinates in place, from alignment routines such as GetBestRMS, which optimize superposition. Choose the operation from the scientific question, not from whichever function returns the smallest value. RDKit molecular alignment API.

For interpreting the original result columns, use the separate guide to Vina affinity, RMSD, and pose selection.

Preserve the reference before preparing docking inputs

Retain the experimental ligand coordinates before generating another conformer, optimizing the ligand, or preparing its docking input. A newly generated low-energy conformer is not the experimental reference.

Keep four artifacts distinguishable:

  1. The original complex, including the selected model, chain, ligand instance and alternate-location decision.
  2. The reference ligand with experimentally observed heavy-atom coordinates and reviewed chemical identity.
  3. The prepared docking receptor and ligand.
  4. The output poses and the complete execution record.

If preparation preserved the receptor coordinates, a separate receptor fit may be unnecessary. Check that assumption. If coordinates changed, align appropriate corresponding receptor atoms, then apply the same transformation to the associated ligand coordinates. Do not fit the ligand independently afterward.

For different receptor structures, document the atom selection used for the fit and the local structural differences. Cross-docking into another conformation is a different experiment from cognate redocking.

Establish chemical identity before atom correspondence

An atom index is a position in a file, not a chemical identity. Exporting a ligand can reorder its atoms without changing its structure. Symmetry can also permit more than one chemically equivalent correspondence.

The spyrmsd study explains why index-only matching can inflate RMSD, while unrestricted distance-based assignment can pair atoms in ways that violate molecular connectivity. Graph-based correspondence addresses that problem by respecting the molecular structure. Symmetry correction is not permission to exchange arbitrary atoms simply because they are nearby. Meli and Biggin, 2020.

Before comparing coordinates, confirm the same covalent structure, protonation state, tautomer and specified stereochemistry. Do not neutralize or change bond orders merely to make a failed comparison run.

The example below deliberately requires the same RDKit molecular representation after removing ordinary explicit hydrogens. It considers graph-preserving alternatives with chirality checking. It does not implement a separate resonance-equivalence policy for differently represented terminal groups. If your study requires that policy, declare and test it explicitly rather than mixing the resulting values with this strict metric.

Export the poses without reconstructing chemistry by guesswork

For docking outputs that retain Meeko's required SMILES and atom-mapping information, Meeko can export the ligand poses to SDF:

mk_export.py vina_out.pdbqt -s poses.sdf

Keep the original PDBQT and SDF files. An arbitrary PDBQT without the required identity information is not interchangeable with this input. The documented export route is intended to preserve molecular information rather than infer it solely from coordinates. AutoDock Vina export example, Meeko official usage.

Check that the number and order of exported records correspond to the retained output modes. Do not assume that the requested maximum number of modes equals the number actually written.

Run a strict fixed-frame RMSD calculation

Save the following as fixed_frame_rmsd.py. The example was checked with Python 3.12, RDKit 2025.09.6 and Meeko 0.7.1 for the preceding export. These are the tested versions, not a claim that they are the newest releases.

The script reads one reference SDF record and a multi-record pose SDF. Coordinates must already be in the same frame and expressed in angstroms. It never aligns, minimizes or repairs them. A canonical-identity check precedes full graph matching, and a mapping-limit error stops the calculation instead of returning an incomplete minimum. The matching options are documented in the RDKit molecule API.

"""Strict same-state, fixed-frame heavy-atom RMSD for the editorial example.

No coordinate alignment, state standardization or file repair is performed.
RDKit 2025.09.6 was used for the retained verification. This is an external
analysis example, not an added MolNexus or BioChemIntelli platform feature.
"""

import argparse
import json
import math

from rdkit import Chem, rdBase


def heavy_molecule(mol):
    if mol is None or mol.GetNumConformers() != 1:
        raise ValueError("Expected one readable molecule with one conformer")
    result = Chem.RemoveHs(Chem.Mol(mol))
    if not result.GetNumAtoms() or any(
        atom.GetAtomicNum() <= 1 for atom in result.GetAtoms()
    ):
        raise ValueError("Empty molecule, dummy atoms or retained hydrogen")
    if len(Chem.GetMolFrags(result)) != 1:
        raise ValueError("Use a reviewed single ligand, not a salt or mixture")
    for atom in result.GetAtoms():
        atom.SetAtomMapNum(0)
    conf = result.GetConformer()
    if not conf.Is3D() or not all(
        math.isfinite(value) for xyz in conf.GetPositions() for value in xyz
    ):
        raise ValueError("Expected finite 3D coordinates in angstroms")
    return result


def fixed_frame_rmsd(reference, pose, max_maps=100000):
    if max_maps < 1:
        raise ValueError("max_maps must be positive")
    ref, probe = heavy_molecule(reference), heavy_molecule(pose)
    if Chem.MolToSmiles(ref, isomericSmiles=True) != Chem.MolToSmiles(
        probe, isomericSmiles=True
    ):
        raise ValueError("Different molecular identity or chemical state")
    maps = ref.GetSubstructMatches(
        probe, uniquify=False, useChirality=True, maxMatches=max_maps + 1
    )
    if not maps or len(maps) > max_maps:
        raise ValueError("No complete mapping, or mapping limit exceeded")
    xyz_ref = ref.GetConformer().GetPositions()
    xyz_pose = probe.GetConformer().GetPositions()

    def squared_error(mapping):
        return sum(
            sum((xyz_pose[i][k] - xyz_ref[j][k]) ** 2 for k in range(3))
            for i, j in enumerate(mapping)
        )

    best_map = min(maps, key=squared_error)
    return {
        "rmsd_angstrom": math.sqrt(squared_error(best_map) / ref.GetNumAtoms()),
        "heavy_atom_count": ref.GetNumAtoms(),
        "mapping_count": len(maps),
        "probe_to_reference_atom_map": list(enumerate(best_map)),
    }


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("reference_sdf")
    parser.add_argument("poses_sdf")
    args = parser.parse_args()
    references = list(Chem.SDMolSupplier(args.reference_sdf, removeHs=False))
    if len(references) != 1:
        raise ValueError("Reference SDF must contain exactly one record")
    poses = list(Chem.SDMolSupplier(args.poses_sdf, removeHs=False))
    if not poses:
        raise ValueError("Pose SDF contains no records")
    results = [
        {"sdf_record": i, **fixed_frame_rmsd(references[0], pose)}
        for i, pose in enumerate(poses, start=1)
    ]
    print(json.dumps({
        "rdkit_version": rdBase.rdkitVersion,
        "metric": "strict same-state symmetry-aware heavy-atom fixed-frame RMSD",
        "alignment": "none; input coordinate frames must already agree",
        "atom_indices": "zero-based indices after RemoveHs, before any atom reordering",
        "results": results,
    }, indent=2))


if __name__ == "__main__":
    main()

Run it in an environment containing RDKit:

python fixed_frame_rmsd.py reference.sdf poses.sdf

The output includes the RMSD, heavy-atom count, number of candidate mappings and the selected atom pairs. Those indices refer to the molecules after hydrogen removal. Preserve the input files and software version with the output so that the mapping remains interpretable.

The calculation takes the square root of the mean squared distance across matched heavy atoms, selecting the smallest value among the admitted mappings. It does not optimize any coordinates. It also does not establish that a shared coordinate frame, chemical representation or experimental reference was scientifically appropriate; those remain input-review responsibilities.

Check the calculation with a retained 1IEP example

BioChemIntelli retains an AutoDock Vina 1.2.7 run based on the c-Abl kinase and imatinib tutorial case. The original experimental complex is RCSB PDB 1IEP.

For this retained run, the receptor box center was 15.19, 53.903, 16.917 Å, with dimensions 20, 20, 20 Å, Vina scoring, exhaustiveness 32, seed 42, one CPU, a maximum of nine requested modes and an energy range of 3 kcal/mol. The saved output contains five poses.

Recalculating those saved poses against the retained reference SDF gives:

Retained mode Fixed-frame heavy-atom RMSD to reference, Å
1 0.375
2 1.501
3 12.289
4 12.119
5 12.499

All five comparisons use 37 heavy atoms and four admitted mappings. These are measured results from the retained files, not newly generated docking predictions or expected values for every installation. The top pose's nonzero reference RMSD also illustrates why its zero-valued Vina mode-to-mode columns cannot answer the reference-recovery question.

This is one pose-recovery control. It does not establish screening enrichment, selectivity, performance on another target, or a general accuracy advantage for a graphical interface.

Test the failure modes before trusting a batch

A small set of deliberately constructed controls makes the metric easier to audit.

Control Expected behavior
Reference compared with an identical copy RMSD zero within numerical precision
Same atoms reordered, coordinates attached to the same atoms Same result after graph mapping
Identical reference translated exactly 5 Å along one axis Fixed-frame RMSD of 5 Å, not zero
Explicit hydrogens added without changing heavy-atom coordinates Same heavy-atom result
Different protonation state, stereoisomer or covalent identity Rejected by this strict example
Mapping count exceeds the declared cap Error requiring review, not a reported approximate minimum

For a reference-recovery benchmark, report the RMSD of the top-ranked pose separately from the best RMSD among all generated poses. Selecting the pose with the lowest reference RMSD uses knowledge unavailable for a genuinely unknown binding mode.

Keep the acceptance rule visible. A chosen cutoff is a decision convention for a defined test, not evidence that every pose below it is chemically plausible or that the compound binds. Follow reference recovery with the broader docking-pose audit and target-specific controls.

Where MolNexus fits

MolNexus provides a Windows desktop workflow for preparation review, Vina execution, pose inspection, exports and local history. The Python calculation in this article is a separate analysis step, not a claim that MolNexus automatically validates poses or includes this script.

If you are evaluating a desktop workflow, decide in advance which reference control, files, settings and exports you need to inspect. A favorable RMSD alone cannot decide whether the software fits your working process.

Plan a controlled two-run workflow evaluation.

Frequently asked questions

References

  1. Center for Computational Structural Biology. Basic docking AutoDock Vina documentation Ligand preparation, standard result table, and export section
  2. RDKit contributors. RDKit::MolAlign namespace reference Official RDKit API documentation CalcRMS versus GetBestRMS/alignMol
  3. Rocco Meli, Philip C. Biggin. spyrmsd: symmetry-corrected RMSD calculations in Python Journal of Cheminformatics 12, 49 (2020) DOI: 10.1186/s13321-020-00455-2 Standard RMSD; Hungarian algorithm; Graph isomorphisms
  4. Forli Lab. Meeko: interface for AutoDock Official Meeko repository README Usage: ligand parameterization and output export
  5. RDKit contributors. rdkit.Chem.rdchem module Official RDKit Python documentation Mol.GetSubstructMatches, conformers and atom operations
  6. RCSB Protein Data Bank. 1IEP: Crystal structure of the c-Abl kinase domain in complex with STI-571 RCSB PDB, 1IEP Structure title, ligand and primary citation
Continue exploring
View topic →