Create a .cube file of the molecular electrostatic potential using ORCA.

mep.py

"""
(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(["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

water.inp

# ORCA input file used to generate the files needed for the MEP calculation.
! bp86 def2-svp keepdens xyzfile

* xyz 0 1
O          0.00000        0.00000        0.11779
H          0.00000        0.75545       -0.47116
H          0.00000       -0.75545       -0.47116
*

Thanks for making a fantastic script. This might be my config specific, but changing line 128 as follows permits multiprocessing MPI calls on a wider range of architectures (i.e. without the "you must use the full path of the executable" error):

    subprocess.check_call([orca_vpot, f"{basename}_mep.inp"])

I am unsure if this is a typo or not (!)

As mentioned orca 5.0 version doesnot create .scfp file rather creates .density file. So, how do I use it to calculate electrostatic potential of a given system

Hi @ssarix I hope you were able to calculate MEP can you guide me through. As initially I used an older orca version to calculate MEP but that it taking quite a long time. Can you guide me how to calculate MEP using orca version 5.0 which generates densities file instead of .scfp?

Hi, i also use version orca 5.0.4 . Today i have use the out from a calculation done using "! B3LYP Freq tightopt 6-31G* ". and _python d:\orca5\mep.py lactal5a_k2 80. lactal5a_k2.* are the files for and from the orca calculation. Resulting file are lactal5a_k2_mep. inp < ...out, ...xyz and ....cube. The .cube is opend in chimerax and give the following image:

@AJ2397 Did you add the keepdensity keyword in your input? Please see the example on the water molecule. I just ran it with version 5.0.4, and it works.

@AJ2397 Add keepdens to your input first. After the calculation is done you must create two cube files to visualize MEP maps. First one is electron density cube. You can create this cube via terminal using this command:
orca_plot filename.gbw -i
Choose 1. Electron density and gaussian cube as cube type. This creates density cube.
After that you can use mep.py script to create electrostatic potential cube. The command is:
python mep.py filename 80

Filename does not need an extension there. Then you can open this two cubes via VMD or Chemcraft to visualize MEP maps.

@ssarix thankyou for your reply but i want to calculated MEP between two points like I actually want numbers. Orca version 4.0 generated .scfp file and using this and 3 other files i could generate a vpot.out file. But with version 5.0 there is error as this .scfp file is not created? Can you suggest something in this regard.

@AJ2397 I did not quite understand what are you trying to achieve. What is the error code that you encounter ?

@AJ2397 If you want to quantitatively analyze the MESP, then I would recommend you use Multiwfn (free). You can convert ORCA gbw files to molden files which can then be read into Multiwfn, which will compute the MESP and provide statistical analysis of the surface automatically. It also allows you to compute the MESP at specific points, among other things. Conveniently, all of this is covered in the extensive tutorials provided by the developer. Highly recommended, just make sure to use the “developer” version as the “official” version is very old and not as efficient. I cannot comment on the relative efficiency of it versus this Python script, though, but it is definitely more flexible. I also think you can hook it up to ORCA to compute the MESP more efficiently, but you’ll have to read the forums to get that working.

…

On Thu, Feb 15, 2024 at 9:51 AM ssarix ***@***.***> wrote: ***@***.**** commented on this gist. ------------------------------ @AJ2397 <https://github.com/AJ2397> I did not quite understand what are you trying to achieve. What is the error code that you encounter ? — Reply to this email directly, view it on GitHub <https://gist.github.com/mretegan/5501553#gistcomment-4903032> or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABA73NKBQCZCVYLAUHBRGJ3YTYOGZBFKMF2HI4TJMJ2XIZLTSKBKK5TBNR2WLJDUOJ2WLJDOMFWWLO3UNBZGKYLEL5YGC4TUNFRWS4DBNZ2F6YLDORUXM2LUPGBKK5TBNR2WLJDHNFZXJJDOMFWWLK3UNBZGKYLEL52HS4DFVRZXKYTKMVRXIX3UPFYGLK2HNFZXIQ3PNVWWK3TUUZ2G64DJMNZZDAVEOR4XAZNEM5UXG5FFOZQWY5LFU42TKMBRGU2THJ3UOJUWOZ3FOKTGG4TFMF2GK> . You are receiving this email because you commented on the thread. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub> .