PyMOL Commands for Docking and MD: Beginner Cheat Sheet

Answer: Most docking and molecular dynamics work in PyMOL runs on about fifteen commands. You use load and split_states to bring poses in, select with byres, within and organic to isolate a binding site, show for representations, distance in mode 2 for polar contacts, align, super or cealign to superimpose, and png with ray for figures.

PyMOL has hundreds of commands. A student running docking and MD uses a small, stable subset of them, week after week. This cheat sheet is that subset, taught as typed commands rather than menu clicks, because commands are what you can paste into a .pml script, repeat on the next target, and put in a methods section. Every command below is checked against the official PyMOL reference.

Which PyMOL commands does a docking or MD student actually need?

The working set is small. If you learn the fifteen commands in the table below, you can inspect any docked pose or any simulated structure without touching a menu.

Everything else in PyMOL is a refinement of these. If you are still choosing where PyMOL fits in your training, the computational biology skills roadmap shows which stage it belongs to.

How do you load a receptor and a docked pose into the same session?

Load the receptor first, then the poses, and give both an explicit object name so your later selections are unambiguous.

load receptor.pdb, receptor
remove solvent
h_add polymer

AutoDock Vina writes its results as a multi-model .pdbqt file, which is not a format PyMOL reads cleanly. Convert it first with Open Babel:

obabel vina_out.pdbqt -O poses.sdf

Add -m to that command and Open Babel splits the models into separate numbered files instead of one multi-molecule file. Our guide to molecular file format conversion with Open Babel covers the flags in detail. Back in PyMOL, load the converted file and separate the poses:

load poses.sdf, poses
split_states poses, prefix=pose
delete poses

split_states takes the arguments object [, first [, last [, prefix ]]], so split_states poses, 1, 3, prefix=pose gives you only the top three. That matters because Vina writes its modes in rank order, and you rarely inspect all nine.

For a finished MD run, PyMOL loads the topology and the trajectory in two steps. The official reference gives the GROMACS case directly:

load md_0_1.gro, mytraj
load_traj md_0_1_center.xtc, mytraj

A long trajectory will exhaust your RAM if PyMOL builds every state. Set set defer_builds_mode, 3 before the load_traj line to avoid that. For serious trajectory work most groups still prefer VMD, and our walkthrough on visualising a GROMACS trajectory in VMD explains why. PyMOL is where you go to make the final frame look publishable.

How does PyMOL selection algebra isolate a binding site?

Selection algebra is the single skill that separates a slow PyMOL user from a fast one. Instead of clicking residues, you describe them.

The vocabulary you need is short:

• polymer is the protein or nucleic acid, organic is the small molecule, solvent is water.
• chain A, resi 45-60, resn HIS and name CA match the identifier columns of the file.
• within and around are the distance operators. s1 within 4 of s2 keeps atoms of s1; s1 around 4 returns neighbours and excludes s1 itself.
• byres expands an atom-level hit to the whole residue, which is almost always what you want for a pocket.

Put together, the one line that defines a binding site is:

select site, byres (polymer within 5 of organic)

Watch the priority rule documented on the PyMOL Wiki Selection Algebra page: the by operators bind weakly, so byres s1 or s2 is read as byres (s1 or s2). Parenthesise when you mean otherwise. To see what you actually caught, print it:

iterate site and name CA, print(chain, resi, resn)

That gives you the residue list you will paste into your results table, derived from coordinates rather than from squinting at a figure. Compare it against the residues an interaction-diagram tool reports; our post on 2D protein-ligand interaction diagrams with LigPlot+ and PLIP shows how the two views should agree.

Which representation commands do you need?

Four lines take a raw session to a readable one:

hide everything
show cartoon, polymer
show sticks, organic
show sticks, site and not (name C+N+O)

The last line is the trick worth memorising. It shows pocket side chains while hiding the backbone atoms that would otherwise clutter the image. Then colour by element so the ligand reads correctly:

color grey80, polymer
color yellow, organic and elem C
util.cnc organic

util.cnc colours non-carbon atoms by element and leaves your carbon colour alone. If you want to show the pocket shape, add a translucent surface rather than an opaque one:

show surface, site
set transparency, 0.5

How do you find polar contacts and measure distances?

The distance command does both jobs. Its signature is distance name, selection1, selection2, cutoff, mode, and mode 2 restricts the output to polar contacts:

distance polar_contacts, receptor, pose1, 3.5, mode=2

A 3.5 Angstrom cutoff is the common convention for a hydrogen bond between heavy atoms. Hide the labels with hide labels, polar_contacts if the numbers crowd the figure, and recolour the dashes with color magenta, polar_contacts.

When you need one number rather than a picture, use get_distance, which returns the value instead of creating an object:

get_distance /receptor//A/189/OD1, /pose1////O1

One caveat that costs students time: the PyMOL reference states that distance modes 5 to 8, which cover pi-pi and pi-cation interactions, “are available only in incentive PyMOL”. The free build gives you polar contacts and plain distances, not aromatic interaction detection. Use PLIP or Discovery Studio Visualizer for those, and read your contact list back against the interpretation of docking results before you claim a binding mode.

Also remember that mode 2 needs hydrogens on polar atoms to behave sensibly. A structure straight from the PDB has none, which is why h_add polymer appears in the loading block above.

align, super or cealign: which should you use to superimpose two structures?

All three superimpose a mobile structure onto a target and report an RMSD, but they pair atoms differently, and picking the wrong one is the most common source of a misleading redocking RMSD.

The official PyMOL reference describes align as a command that “performs a sequence alignment followed by a structural superposition”, then runs refinement cycles to reject outliers. Its defaults are an outlier cutoff of 2.0 Angstroms over a maximum of 5 cycles, and the same reference notes that it “does a good job on proteins with decent sequence similarity (identity >30%)”. Below that identity, sequence-based pairing starts to fail.

Read the last row twice. The PyMOL reference flags it explicitly: in cealign the mobile and target arguments are swapped relative to align and super. Getting it backwards still produces an RMSD, just not the one you meant.

cealign implements the combinatorial extension algorithm published by Shindyalov and Bourne in Protein Engineering in 1998 (volume 11, issue 9, pages 739 to 747), which builds an alignment from short aligned fragment pairs rather than from sequence. That is why it survives the twilight zone where sequence-based methods stop working.

For a redocking control, you usually want the all-atom RMSD with no outlier rejection at all, because rejecting outliers is exactly what hides a bad pose:

align pose1, crystal_ligand, cycles=0, transform=0

transform=0 reports the number without moving anything, which keeps your scene intact. Our step-by-step on superimposing two proteins in PyMOL works through the protein case in full.

How do you save a session and a publication-ready image?

Save the session before you save the image, because the session is what lets you or your supervisor reopen the exact scene:

save complex_analysis.pse

Then render. The png command accepts physical units when you also give a dpi value, which is the fastest way to hit a journal figure specification:

bg_color white
set ray_opaque_background, 0
png figure1.png, width=12cm, dpi=300, ray=1

ray=1 runs the ray tracer before writing, and ray_opaque_background, 0 gives a transparent background that drops cleanly onto a poster or a slide. For a presentation look, set ray_trace_mode, 1 adds a black outline to the rendering. Our guide to publication-quality figures in PyMOL covers lighting and depth cue settings.

Save the commands themselves as a .pml file and rerun the whole scene with pymol -cq figure.pml in a terminal. That is how the same figure gets regenerated for every compound in your series without repeating a single click.

Does the free open-source build do all of this?

Almost. The open-source PyMOL repository maintained by Schrodinger describes itself as the “Open-source foundation of the user-sponsored PyMOL molecular visualization system” and is published under a BSD-like license, while the commercial product, called Incentive PyMOL, is sold separately with support. Every command in this cheat sheet runs in the open-source build except the pi interaction distance modes noted earlier.

The practical difference for a student is packaging rather than capability: the paid build ships as a one-click installer with bundled plugins, the free one is installed through conda or built from source. Our walkthrough on installing PyMOL for free covers both routes. Before you recommend the paid build to a lab, check whether anything you rely on is actually Incentive-only.

Troubleshooting: the errors that stop beginners

Frequently asked questions

What is the difference between within and around in PyMOL?

s1 within 4 of s2 returns atoms belonging to s1 that lie within 4 Angstroms of s2. s1 around 4 returns the neighbouring atoms and excludes s1 itself. For a binding pocket you normally want byres (polymer within 5 of organic).

Can PyMOL open AutoDock Vina output directly?

Not reliably. The .pdbqt format carries partial charges and atom types in columns PyMOL does not expect. Convert to SDF or PDB with Open Babel, then use split_states to separate the binding modes.

Which command should I use for a redocking RMSD?

Use align with cycles=0 and transform=0 on the ligand selections. Outlier rejection is on by default and will flatter a pose that is genuinely wrong.

Can PyMOL replace VMD for molecular dynamics analysis?

No. PyMOL reads a GROMACS trajectory through load_traj and renders frames beautifully, but VMD and the GROMACS analysis tools are built for measuring trajectories. Use PyMOL for representative structures and figures.

Is open-source PyMOL good enough for a thesis?

Yes. The open-source build handles loading, selections, contacts, alignment and ray-traced rendering. Only a few features, including the pi-pi and pi-cation distance modes, are restricted to Incentive PyMOL.

How do I run the same PyMOL scene on twenty compounds?

Put the commands in a .pml script and run pymol -cq scene.pml from the shell in a loop. That is the entire reason to learn PyMOL as commands rather than menus.

Where this fits in your workflow

PyMOL sits between docking and reporting. You dock, you inspect the pose in PyMOL, you check the contacts, and you carry the residue list into your results. After MD, you extract a representative structure and come back to PyMOL for the figure. The skills roadmap places it alongside structure preparation, and a residue mutation in PyMOL is often the next thing you will be asked to do.

This guide was written by the StemSkills Lab team, whose members have more than ten years of combined work in sequence and structural bioinformatics, drug discovery and design, and multiscale molecular modeling.