Tutorial: How to plot High-Quality Molecular Electrostatic Potential Surfaces (MEPS)

Intro

In this tutorial, you will learn how to visualize Molecular Electrostatic Potential Surfaces (MEPS) using Chimera 1.17.3. The process involves multiple steps, including generating electronic density data with ORCA, post-processing with a custom Python script, and finally, using Chimera to create high-quality MEPS visualizations.

Table of contents

Required software

The ORCA calculations and python data manipulations have been tested on Ubuntu 20.04 WSL virtual machine (Windows 11). Chimera UCSF has been installed and run on Windows.

An overview on MEPS

What Are Molecular Electrostatic Potential Surfaces (MEPS)?

Molecular Electrostatic Potential Surfaces (MEPS) are graphical representations that map the electrostatic potential (ESP) generated by a molecule’s charge distribution onto a surface. This potential is influenced by contributions from both the molecular nuclei and the surrounding electron density. MEPS provide a visual method to analyze how this distribution influences molecular interactions with other molecules, ions, or electric fields, offering insights into reactivity, intermolecular interactions, and binding affinities.

The electrostatic potential at any point around a molecule represents the energy that a positive test charge would experience at that point. It is determined by the positive influences of the nuclei and the negative contributions from the electron cloud. This makes MEPS useful for assessing chemical behavior, especially in terms of how different regions of a molecule might attract or repel other charged particles.

For an introductory review on MEPS, see “Molecular Electrostatic Potentials: Significance and Applications” and “The electrostatic potential: an overview” by Politzer and Murray, which explore the role of Molecular Electrostatic Potentials in understanding chemical reactivity, noncovalent interactions, and their relationship with electron density.

Physical intuition behind MEPS

Physically, the electrostatic potential at any point around a molecule is the energy experienced by a positive test charge due to the electric field generated by the molecule’s electron cloud and nuclei. A MEPS visually represents this potential, with regions of the surface color-coded to indicate areas of positive and negative potential:

• Positive regions (often colored blue) indicate areas where the molecule has a deficiency of electron density, typically around regions with partially positive atoms or functional groups. These areas are more likely to attract nucleophiles or negatively charged species.
• Negative regions (often colored red) correspond to areas with an excess of electron density, usually around electronegative atoms or groups. These regions are likely to interact with electrophiles or positively charged species.

By analyzing the MEPS, one can predict sites on a molecule that are more likely to participate in certain types of chemical interactions, aiding in the design and understanding of molecular behavior.

Note: MEPS are commonly visualized by mapping the potential onto a surface of constant electron density. However, it is essential to recognize that the MEP is a distinct property from the electron density. It reflects the combined effects of both electrons and nuclei, rather than solely the density of electrons, thereby providing a more comprehensive understanding of molecular interactions.

Applications of MEPS

MEPS are widely used in computational chemistry for several purposes:

• Predicting Reactivity: MEPS can identify nucleophilic and electrophilic sites on a molecule, helping to predict where reactions are likely to occur.
• Understanding Intermolecular Interactions: By comparing the MEPS of different molecules, one can infer how they might interact, which is particularly useful in drug design, host-guest chemistry, and studying molecular recognition.
• Exploring Solvent Effects: MEPS can help explain how molecules interact with solvents or in biological environments, aiding in the understanding of solvation phenomena.
• Designing Catalysts: In catalysis, MEPS can guide the design of catalysts by identifying regions of a molecule that will interact favorably with substrates or other reaction components.

Key Parameters and Considerations When Using MEPS

When generating and interpreting MEPS, several key factors must be considered:

1. Choice of Isosurface: The electron density isosurface used for mapping the electrostatic potential significantly influences the MEPS. A common choice is an isosurface corresponding to a density of 0.002 or 0.004 e/ų, which typically represents the van der Waals surface of the molecule.
2. Level of Theory: The quality of the MEPS depends on the level of quantum chemical theory used in the calculations. Higher levels of theory (e.g., DFT with a good functional and basis set) provide more accurate electron densities and potentials.
3. Grid Resolution: The resolution of the grid used to calculate the electrostatic potential affects the smoothness and accuracy of the MEPS. Finer grids (i.e., higher grid point densities) yield more detailed surfaces but require more computational resources.
4. Color Mapping: The range and choice of colors used to represent different potential values on the MEPS can influence how easily one can interpret the surface. It’s important to use a color scheme that accurately reflects the potential values and aids in visual interpretation.
5. Solvent Effects: Including solvent effects in the calculation (e.g., via implicit solvent models) can significantly alter the electrostatic potential and, consequently, the MEPS. For molecules in solution, it’s crucial to account for these effects to obtain realistic surfaces.
6. Interpretation Caution: While MEPS are powerful tools, they should be interpreted with care. The visual representation of potential can sometimes be misleading if not correlated with chemical intuition or further computational analysis. For example, local minima or maxima in potential do not always correspond directly to chemical reactivity.

Step by step guide

Step 1: Running the ORCA Calculation

The first step involves performing a quantum chemistry calculation using ORCA to obtain the electronic density data. Carefully choose the appropriate geometry based on the specific physical problem you are studying, as this choice can significantly impact the accuracy of your results. Be sure to include the “keepdens” keyword in your ORCA input, as demonstrated in the sample below, to ensure the density data is retained for further analysis:

! M062X ma-def2-TZVP Normalprint Printbasis PrintMOs keepdens
%maxcore 1500
%PAL NPROCS 16 END
%CPCM SMD TRUE
      SMDSOLVENT "WATER"
END
*xyzfile 0 2 your_structure.xyz

Key Points:

• Functional and Basis Set: Choose an appropriate functional and basis set for your molecule (e.g. for anions, use augmented basis sets e.g. ma-def2-TZVP)
• Solvation Model: Unless you plan to study the charge distribution in vacuum, make sure to specify an implicit solvent for your calculation
• Computational Resources: Adjust the %maxcore and %PAL NPROCS directives according to your available computational resources.

Step 2: Generating Gaussian Style Electronic Density Cube File with ORCA_Plot

After completing the ORCA calculation, you will generate a cube file format suitable for Chimera using orca_plot.

Run the following command:

orca_plot your_structure.gbw -i

In the interactive menu, select the following options:

• 1 – Enter type of plot → 2 – (scf) electron density → Confirm
• 4 – Enter the number of grid intervals → Type the value (good quality can be obtained from a value of 80 or higher)
• 5 – Select output file format → 7 – 3D Gaussian cube
• 10 – Generate the plot

This will generate the eldens.cube file.

Although this process can be done manually as shown, it can also be automated with the following one-liner:

echo -e "1\n2\ny\n4\n80\n5\n7\n10\n11" | orca_plot your_structure.gbw -i

Step 3: Post-Processing with MEP Script to Generate ESP Volumetric Data

To create the electrostatic potential (ESP) map, you will use a custom Python script (mep.py) that can be found in the appendix of this tutorial. Ensure the following files from your ORCA calculation are in your working directory:

• your_structure.gbw
• your_structure.scfp
• your_structure_mep.xyz
• your_structure_mep.out

Make the mep.py script executable and run it with:

python mep.py your_structure 80 64

Replace your_structure with the appropriate file name prefix. The parameters 80 and 64 refer to the grid points and number of processes, respectively.

Note: the calculation of the MEP scales fairly well even at high processor counts, so feel free to exploit any available cores.

Step 4: Visualizing MEPS in Chimera

You will now use Chimera to visualize the ESP data. You can follow the following video tutorial for creating the MESP visualization (taken from Plotting molecular electrostatic potential surfaces (MEPS) with ORCA – Henrique’s Blog (henriquecastro.info)):

Suggested settings for three-color colouring:

• Orange: #ffff7ae10000
• Green: #5c29c28f8000
• Light blue: #7ae151ebffff

Keep in mind that you can save the session as a Python file, allowing you to easily restore it later and continue from where you left off.

This process should produce a result similar to the following:

Appendix

"""
(c) 2013, 2019, 2022 Marius Retegan
License: BSD-2-Clause
Description: Create a .cube file of the molecular electrostatic
             potential (MEP) using ORCA.
Run: python mep.py basename npoints nprocs (e.g. python mep.py water 40 2).
Arguments: basename - file name without the extension;
                      this should be the same for the .gbw and .scfp.
           npoints  - number of grid points per side
                      (80 should be fine)
           nprocs - number of processors for parallel calculations
"""
#!/usr/bin/env python
def read_xyz(xyz):
    atoms = []
    x, y, z = [], [], []
    with open(xyz) as fp:
        # Skip the first two lines.
        next(fp)
        next(fp)
        for line in fp:
            data = line.split()
            atoms.append(data[0])
            x.append(float(data[1]))
            y.append(float(data[2]))
            z.append(float(data[3]))
    return atoms, np.array(x), np.array(y), np.array(z)
def read_vpot(vpot):
    v = []
    with open(vpot) as fp:
        next(fp)
        for line in fp:
            data = line.split()
            v.append(float(data[3]))
    return np.array(v)
if __name__ == "__main__":
    import os
    import shutil
    import subprocess
    import sys
    import numpy as np
    ang_to_au = 1.0 / 0.5291772083
    elements = [None,
         "H", "He",
         "Li", "Be",
         "B", "C", "N", "O", "F", "Ne",
         "Na", "Mg",
         "Al", "Si", "P", "S", "Cl", "Ar",
         "K", "Ca",
         "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
         "Ga", "Ge", "As", "Se", "Br", "Kr",
         "Rb", "Sr",
         "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd",
         "In", "Sn", "Sb", "Te", "I", "Xe",
         "Cs", "Ba",
         "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
         "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
         "Tl", "Pb", "Bi", "Po", "At", "Rn",
         "Fr", "Ra",
         "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No",
         "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Uub"]
    basename = sys.argv[1]
    xyz = f"{basename}.xyz"
    if not os.path.isfile(xyz):
        sys.exit("Could not find the .xyz. To quickly generate one for "
                 f"your molecule run: echo 11 | orca_plot {basename}.gbw -i.")
    atoms, x, y, z = read_xyz(xyz)
    try:
        npoints = int(sys.argv[2])
    except ValueError:
        sys.exit(f"Invalid number of points: {sys.argv[2]}")
    try:
        nprocs = int(sys.argv[3])
    except IndexError:
        nprocs = 1
    except ValueError:
        sys.exit(f"Invalid number of cpus: {sys.argv[3]}")
    natoms = len(atoms)
    extent = 7.0
    xmin = x.min() * ang_to_au - extent
    xmax = x.max() * ang_to_au + extent
    ymin = y.min() * ang_to_au - extent
    ymax = y.max() * ang_to_au + extent
    zmin = z.min() * ang_to_au - extent
    zmax = z.max() * ang_to_au + extent
    with open(f"{basename}_mep.inp", "w") as fp:
        fp.write(f"{nprocs:d}\n")
        fp.write(f"{basename}.gbw\n")
        fp.write(f"{basename}.scfp\n")
        fp.write(f"{basename}_mep.xyz\n")
        fp.write(f"{basename}_mep.out\n")
    with open(f"{basename}_mep.xyz", "w") as fp:
        fp.write(f"{npoints**3:d}\n")
        for ix in np.linspace(xmin, xmax, npoints, True):
            for iy in np.linspace(ymin, ymax, npoints, True):
                for iz in np.linspace(zmin, zmax, npoints, True):
                    fp.write(f"{ix:12.6f} {iy:12.6f} {iz:12.6f}\n")
    orca_vpot = shutil.which("orca_vpot")
    if orca_vpot is None:
        sys.exit(f"Could not find the orca_vpot executable.")
    subprocess.check_call([f"{orca_vpot}", f"{basename}_mep.inp"])
    vpot = read_vpot(f"{basename}_mep.out")
    with open(f"{basename}_mep.cube", "w") as fp:
        fp.write("Generated with ORCA\n")
        fp.write(f"Electrostatic potential for {basename}\n")
        fp.write(f"{len(atoms):5d}{xmin:12.6f}{ymin:12.6f}{zmin:12.6f}\n")
        xstep = (xmax - xmin) / float(npoints - 1)
        fp.write(f"{npoints:5d}{xstep:12.6f}{0:12.6f}{0:12.6f}\n")
        ystep = (ymax - ymin) / float(npoints - 1)
        fp.write(f"{npoints:5d}{0:12.6f}{ystep:12.6f}{0:12.6f}\n")
        zstep = (zmax - zmin) / float(npoints - 1)
        fp.write(f"{npoints:5d}{0:12.6f}{0:12.6f}{zstep:12.6f}\n")
        for i, atom in enumerate(atoms):
            index = elements.index(atom)
            xi, yi, zi = x[i] * ang_to_au, y[i] * ang_to_au, z[i] * ang_to_au
            fp.write(f"{index:5d}{0:12.6f}{xi:12.6f}{yi:12.6f}{zi:12.6f}\n")
        m = 0
        n = 0
        vpot = np.reshape(vpot, (npoints, npoints, npoints))
        for ix in range(npoints):
            for iy in range(npoints):
                for iz in range(npoints):
                    fp.write(f"{vpot[ix][iy][iz]:14.5e}")
                    m += 1
                    n += 1
                    if (n > 5):
                        fp.write("\n")
                        n = 0
                if n != 0:
                    fp.write("\n")
                    n = 0

Mattia Felice Palermo

Researcher at CIC energiGUNE

Mattia Felice Palermo, Ph.D. in Computational Chemistry, has experience in industrial R&D with corporate and SME companies. His expertise includes molecular dynamics modeling of liquid crystals and polymers, and predicting electrochemical properties using electronic structure methods. He led the computational chemistry group at Green Energy Storage and now is a researcher at the CIC energiGUNE research center.

• Mattia Felice Palermo

  https://modelinghub.xyz/author/mattiafelicepalermo/

  Tutorial: Parametrizing dihedrals
• Mattia Felice Palermo

  https://modelinghub.xyz/author/mattiafelicepalermo/

  Tutorial: multidimensional dihedral PES in ORCA + visualization
• Mattia Felice Palermo

  https://modelinghub.xyz/author/mattiafelicepalermo/

  Step-by-Step Guide to Installing SLURM with memory limit and core affinity on Ubuntu – Single Node