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. Make sure it is added to your PATH.")

    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

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
*

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> .

Hi sir,
i was generating files for MEP.I could generate the cube file. But when i run python mep. py basename 80 command the following error comes.
can anyone help me to sort it out.I am Using ORCA 4.2.1.
mep.py:21: Warning: 'with' will become a reserved keyword in Python 2.6
File "mep.py", line 21
with open(xyz) as fp:
^
SyntaxError: invalid syntax

@mretegan I noticed that when running mep.py with the multiprocess flag, the orca_vpot command fails since it's not called with the full path:

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!     ERROR (ORCA_MAIN): For parallel runs                     !!!
!!!     ORCA has to be called with full pathname                 !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[file orca_tools/qcsys.cpp, line 40]:

[file orca_tools/qcsys.cpp, line 40]:

Applying the following modification fixes it:

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

Hope it helps!

Mattia

Hi mretegan thank you so much for the amazing code it was working for an older versions of orca however orca 6 now is released and i tested but unfortunately it doesn't work anymore it gives an error can you update/revise it based on the new release features of orca6 thank you so much.

@mattiafelice-palermo Thanks for the fix.

@ochemguru Can you try again using the updated script? There was an error in the code, but it was unrelated to the new version.

@mretegan thank you for the reply i tried the new version and still giving the same problem here is the error :
ORCA ELECTROSTATIC POTENTIAL GENERATION

GBW file ... 1.gbw
File with points to evaluate V(r) ... 1_mep.xyz
Output file ... 1_mep.out

Reading the GBW file ... done
Reading the positions ... done (1728000 positions)
[file orca_tools/Tool-GTO-Integrals/gtovpot.cpp, line 1193]: Error: Prescreening matrix could not be read!

[file orca_tools/Tool-GTO-Integrals/gtovpot.cpp, line 1193]: Error: Prescreening matrix could not be read!

Traceback (most recent call last):
File "3.py", line 130, in
vpot = read_vpot(f"{basename}_mep.out")
File "3.py", line 38, in read_vpot
with open(vpot) as fp:
FileNotFoundError: [Errno 2] No such file or directory: '1_mep.out'

This is hard to debug because I don't know what your molecule is. Try the following:

1. Make sure the water example works.
2. Run your molecule using a reasonable number of points per side, 40, for example.

If this works, increase the number of points.

i'm trying on the same molecule as provided (water) with different in theory only (b3lyp) unfortunately no matter how many time i tried, it doesn't work still giving the same problem

I can reproduce the error. I will open a ticket on the forum.

thank you for testing looking forward to solve this problem

Apparently, the bug is known https://orcaforum.kofo.mpg.de/viewtopic.php?f=11&t=11657&p=47529&hilit=vpot#p47529.

You will have to wait for the bug fixing release.

oh i see thank you so much waiting for good things