Electrostatic Surface in PyMOL with APBS: Step by Step

To make an electrostatic surface in PyMOL, convert your PDB file to PQR with PDB2PQR (which adds hydrogens, charges and radii at a chosen pH), solve the Poisson-Boltzmann equation on a grid with APBS, then colour the molecular surface by the resulting potential map. Use a fixed symmetric range such as -5 to +5 kT/e: red is negative, blue is positive.

That is the whole method in one paragraph. The rest of this guide, written by the StemSkills Lab team (10+ years in structural bioinformatics, drug design and molecular modelling), gives you the exact commands, the settings that change the picture, a no-install alternative in ChimeraX, and fixes for the errors students hit most often. The worked example uses ubiquitin, PDB 1UBQ, because it is small, complete and runs in seconds. The same steps apply to your own docking receptor.

What does an electrostatic surface actually show?

It shows the electrostatic potential that a probe charge would feel just outside the protein surface. Regions where the potential is negative (usually rich in Asp and Glu) are coloured red. Regions where it is positive (rich in Lys and Arg) are coloured blue. White is close to neutral.

For a docking or binding-site project this matters because charge complementarity often explains a pose. A carboxylate on your ligand sitting in a blue pocket is a story you can defend in a results section. A positively charged amine pushed against a red wall is the same story. If you are still locating the pocket itself, start with our guide on how to find the protein binding site for docking, then come back here to characterise it.

The potential from APBS is reported in units of kT/e. The APBS documentation gives the conversion: at about 298 K, 1 kT/e is about 25.6 mV. So a surface patch at +5 kT/e sits at roughly +128 mV.

Which tools do I need, and which route should I pick?

You need three things: a charge and radius assignment step (PDB2PQR), a Poisson-Boltzmann solver (APBS) and a viewer (PyMOL or ChimeraX). There are five practical routes, and they are not equally accurate.

Our recommendation for students: use the command-line route for anything that goes into a report, because every parameter is written in a file you can attach to your methods section. Use ChimeraX coulombic when you want a 30-second check during a docking session. The web route is the APBS-PDB2PQR web server, which runs both steps and returns the map for download.

How do I prepare the PDB file before calculating electrostatics?

Electrostatics are sensitive to what is in the file, so decide deliberately what stays.

• Waters: remove them. Crystal waters are not part of the continuum solvent model. PDB2PQR’s --drop-water flag does this for you.
• Ligands and ions: PDB2PQR’s force fields cover proteins and nucleic acids, not drug-like ligands. If you want the receptor’s own potential (the usual case for explaining a docked pose), remove the ligand. If the ligand must be included, see the troubleshooting section below.
• Missing atoms and loops: the PyMOL wiki page for the plugin notes that PDB2PQR will model up to 10% of missing side-chain atoms. It will not build missing loops, and a missing loop near your pocket removes real charges from the picture. Fix those first with our guide to fixing missing residues and loops in a PDB structure.
• Protonation: histidines and buried Asp/Glu change the map. Pick the pH that matches your experiment. Our post on protonation states for docking and MD explains why this choice matters.

How do I run PDB2PQR from the command line?

PDB2PQR is a Python package. Install it with pip, then check that the executable works:

pip install pdb2pqr
pdb2pqr --help

The package installs the command under two names, pdb2pqr and pdb2pqr30, which run the same program. The basic syntax from the PDB2PQR documentation is pdb2pqr [options] --ff={forcefield} {input} {output.pqr}, where the input can be a local file or a four-character PDB ID that is fetched for you.

For the ubiquitin example:

pdb2pqr --ff=AMBER --titration-state-method=propka --with-ph=7.0 \
        --drop-water --keep-chain --apbs-input=1ubq.in \
        1ubq 1ubq.pqr

What each flag does (checked against the current PDB2PQR source):

• --ff=AMBER chooses the charges and radii. Accepted values are AMBER, CHARMM, PARSE, TYL06, PEOEPB and SWANSON. The default is PARSE, which was parameterised for continuum solvation and is also a sound choice. Whichever you use, name it in your methods.
• --titration-state-method=propka predicts pKa values with PROPKA and sets titratable residues for the pH given by --with-ph (default 7.0). The other option is pkaani. Without this flag, residues get their standard states.
• --drop-water removes waters before processing.
• --keep-chain keeps chain IDs in the PQR, which makes later selections in PyMOL easier.
• --apbs-input=1ubq.in writes a ready-to-run APBS input file sized to your molecule.

Read the log. PDB2PQR reports residues it repaired, atoms it could not place and heavy atoms that were missing. Those warnings are your quality control.

How do I run APBS and what settings does it use?

Download APBS from the official APBS releases page and make sure the apbs binary is on your PATH. Then run the input file PDB2PQR wrote:

apbs 1ubq.in

The command form apbs [options] input-file is the one in the APBS documentation. The generated template uses the mg-auto multigrid method (a coarse grid, then a focused fine grid) and writes the potential as an OpenDX map. Open 1ubq.in in a text editor before you run it. In the current template you will see:

• pdie 2.0: protein interior dielectric.
• sdie 78.54: solvent dielectric (water).
• temp 298.15: temperature in kelvin.
• srad 1.4: solvent probe radius in angstroms.
• No ion lines, because the template’s ionic strength is zero.
• write pot dx ...: the name of the output map. APBS adds the .dx extension, so with this template the file is usually 1ubq.pqr.dx. Check the line in your own file.

If you want physiological salt, add a pair of ion lines (one positive, one negative) to the elec block, with the concentration in molar, for example 0.150 for 150 mM. Check the exact syntax against the ELEC section of the APBS input documentation for your version. Salt screens the potential, so the colours become less intense. Report the ionic strength you used.

The science behind this step is the Poisson-Boltzmann equation, solved numerically. The canonical references are Baker et al. (2001), “Electrostatics of nanosystems: application to microtubules and the ribosome”, PNAS 98(18):10037-10041, and the update paper by Jurrus et al. (2018), “Improvements to the APBS biomolecular solvation software suite”, Protein Science 27(1):112-128. Cite both if you use APBS and PDB2PQR in a thesis.

How do I colour the surface in PyMOL?

There are two ways. The plugin is quicker. The manual commands give you more control and a script you can rerun.

Route A: the APBS Electrostatics plugin

1. Load or fetch the structure: fetch 1ubq, async=0.
2. Open Plugin > APBS Electrostatics.
3. Pick the selection (for example polymer) and click Run.

According to the PyMOL wiki, the plugin runs three steps in order: prepare the molecule (PDB2PQR by default), calculate the map with APBS, and show the coloured molecular surface. For standard proteins it is described as a one-click job. Two options are worth knowing. Focus Selection limits the calculation to, say, a binding site to save time. Grid spacing trades resolution for speed and memory. After the run you can change the colour range from the object panel (the ramp object’s levels menu).

Route B: load your own map and set the ramp by hand

load 1ubq.pqr
load 1ubq.pqr.dx, apbs_map
ramp_new apbs_ramp, apbs_map, [-5, 0, 5]
set surface_ramp_above_mode, 1
set surface_color, apbs_ramp, 1ubq
show surface, 1ubq

Three details make this figure correct rather than just colourful:

• A fixed, symmetric range. [-5, 0, 5] forces white at zero and the same scale at both ends. The default range for ramp_new is [-1, 0, 1], which saturates almost every protein. When you compare two proteins or two mutants, use the same range for both or the comparison means nothing.
• surface_ramp_above_mode, 1 colours the surface by the potential one solvent radius above it, on the solvent accessible surface. The PyMOL wiki explains that a ligand or solvent atom never reaches the molecular surface itself, so this is the more meaningful value to show.
• The default colours are red, white and blue, so negative is red and positive is blue without extra arguments. Add a colour bar (the ramp object shows one) and state the unit, kT/e, in the legend.

For lighting, ray tracing and export settings, follow our PyMOL publication-quality figures guide.

Can I do this in ChimeraX without installing APBS?

Yes. ChimeraX has a built-in command:

open 1ubq
coulombic protein

The ChimeraX coulombic documentation lists these defaults: a distance-dependent dielectric (distDep true, dielectric 4.0), an offset of 1.4 Å from each surface vertex, the red-white-blue palette and range -10,10. If hydrogens are missing, it protonates a temporary copy with addh. Nonstandard residues get Antechamber charges (AM1-BCC by default). Using the protein specifier, as above, keeps ligands and waters out of the calculation.

The important difference is physics. coulombic applies Coulomb’s law to the partial charges, with no continuum solvent and no ions. It is fast and often agrees with APBS on the overall pattern. Do not present it as a Poisson-Boltzmann result. The unit also differs: ChimeraX reports kcal/(mol·e), and at 298 K, ±10 kcal/(mol·e) is roughly ±17 kT/e. So a ChimeraX figure and an APBS figure with “the same numbers” on the colour bar are not on the same scale.

ChimeraX can also show an APBS map. Open the .dx file next to the structure and run color electrostatic on the surface with the map option. The ChimeraX color documentation describes it as the same as color sample, but with an offset of 1.4 Å and a range of -10,10 by default. The documentation explains that the 1.4 Å offset shows the values that any atoms binding to the outside of the surface would “see”. To pick between the two programs in general, see PyMOL vs VMD vs ChimeraX.

How do I read the electrostatic surface next to a docked ligand?

Load your docked complex, colour the receptor surface by potential, and show the ligand as sticks. Then make the surface partly transparent (set transparency, 0.4 in PyMOL) so the pose stays visible. Now ask three questions:

1. Do charged ligand groups face complementary patches? An acid in a blue region and an amine in a red region support the pose. The reverse is a warning sign worth a second look at your docking settings.
2. Is the pocket entrance charged? A strongly charged rim can guide a charged ligand in, even if the pocket floor is neutral.
3. Does the pattern hold at other pH values? Rerun PDB2PQR at pH 6.0 and 8.0. If a histidine in the pocket flips the colour, say so in your discussion.

Treat the figure as supporting evidence. It explains a pose; it does not score it. Pair it with the interaction analysis from our guide on how to interpret molecular docking results.

Why did my electrostatics calculation fail? (Troubleshooting)

1. The plugin says it cannot find apbs or pdb2pqr. The plugin calls external programs. Check in a terminal that apbs --help and pdb2pqr --help both run. Then point the plugin’s program location settings at those paths, or start PyMOL from a shell where they are on the PATH. On many systems a conda environment that is active in your terminal is not visible to a PyMOL started from the desktop.

2. The ligand disappears from the PQR. PDB2PQR assigns parameters only to residues its force field knows, so drug-like HETATM records get no charges. If you need the ligand in the calculation, PDB2PQR’s --ligand option accepts a single MOL2 file whose name matches one ligand in the PDB. Otherwise, assign charges with another tool and load a PQR or MOL2 into the plugin with Prepare Molecule unchecked, as the PyMOL wiki describes.

3. “Too few atoms present to reconstruct or cap residue”. This happens when a truncated file contains a stray backbone atom of the next residue. The PyMOL wiki’s example is 2XWU, chain B, residue 954. Delete that atom (remove /2xwu//B/954) and rerun. When larger pieces are missing, rebuild them before you calculate anything.

4. The whole surface is one saturated colour. The range is too narrow. The ramp_new default of -1 to +1 kT/e does this. Set [-5, 0, 5] or edit the ramp levels in the plugin. If the surface is entirely grey or white instead, the surface is probably not coloured by the ramp at all: rerun set surface_color, apbs_ramp, <object> on the right object name.

5. APBS runs out of memory on a large complex. The template sizes the grid for the whole molecule, and the grid grows with the molecule’s volume. In the plugin, increase Grid spacing (the wiki’s virus-capsid example uses 2.0 Å) or use Focus Selection on your pocket. On the command line, analyse the chains around the binding site instead of the whole assembly.

6. Nucleic acids or modified residues are rejected. PDB2PQR expects RNA residue names RA, RC, RG and RU. The wiki gives a one-line fix in PyMOL: alter polymer & resn A+C+G+U, resn = "R" + resn. For modified amino acids such as selenomethionine (MSE), mutate them to the standard residue or supply parameters. Note that ChimeraX coulombic already treats MSE with MET charges.

Frequently asked questions

What colour range should I use for an electrostatic surface?

A symmetric range centred on zero, commonly -5 to +5 kT/e for APBS maps. Keep the same range for every structure you compare, and print the unit on the colour bar.

Is red positive or negative?

In the default red-white-blue scheme used by PyMOL ramps and ChimeraX, red is negative potential and blue is positive. Say so in the figure legend anyway, because reviewers ask.

Is ChimeraX coulombic as good as APBS?

It is a good quick check of the overall charge pattern. It has no solvent screening or ions and uses a different unit, kcal/(mol·e). For a thesis figure or any quantitative comparison, use PDB2PQR and APBS.

Should I include the docked ligand in the calculation?

Usually no. To explain why a ligand binds, you want the receptor’s own potential, so compute it on the protein alone and display the ligand on top. Include the ligand only when you are asking how it changes the field, and then you must give it charges.

Can I calculate electrostatics on an AlphaFold model?

Yes, with care. Low-confidence regions can place charged side chains in arbitrary positions, so check per-residue confidence first using our guide to AlphaFold pLDDT and PAE scores, and trim or flag disordered tails before running PDB2PQR.

Electrostatic surfaces are one of the binding-site skills covered in our learn molecular docking pillar, and they are one step on the computational biology skills roadmap.