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,
I am having the fallowing error
Error (ORCA_VPOT): Cannot read density matrix of dimension 1314 (directorytoscfpFile.scfp).
Chemcraft read the cube file, however, no mep (actually any surface) is read, only structure of the molecule is showed.
thanks for your help
Kamil

Hi,
I am having the fallowing error
Error (ORCA_VPOT): Cannot read density matrix of dimension 1314 (directorytoscfpFile.scfp).
Chemcraft read the cube file, however, no mep (actually any surface) is read, only structure of the molecule is showed.
thanks for your help
Kamil

It was the newby's error! I needed to generate scfb file and install the numpy with conda. Working like a charm right now. Many thanks! Only one suggestion, could you please provide data for the citation of your script !

When the generated cube file is opened in ChemCraft, the charge cloud appears to be of the same colour over the complete complex having different kinds of atoms like Ru, C, H, N, O, S. May I know the reasons for this?

Change to Jmol

…

On Fri, Aug 6, 2021 at 11:04 AM Angelina Mary ***@***.***> wrote: ***@***.**** commented on this gist. ------------------------------ When the generated cube file is opened in ChemCraft, the charge cloud appears to be of the same colour (orange) over the complete complex having different kinds of atoms like C,H,N,O,S. May I know the reasons for this? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <https://gist.github.com/5501553#gistcomment-3848555>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/APTUEDJHLCSPAD4JEVG3D6DT3OQQ5ANCNFSM4HK2EDNA> . 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&utm_campaign=notification-email> .

Really great script, thank you!!
I am not familiar with the cube file format, so the portion of the code that actually generates the file is opaque to me. Would you be able to add comments to help explain exactly what is being done?
Also, is there any way to simply parallelize this code?

Hi, I am using ORCA 5.0.4 and calculation does not create .scfp file, even if I use "keepdens" keyword. The name of my .gbw file is test. When i run the script with command :
python mep.py test 80

I get the following error:

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

Reading the GBW file ... done
Reading the positions ... done (512000 positions)
Error (ORCA_VPOT): Cannot read density matrix of dimension 24 (test.scfp)
Traceback (most recent call last):
File "C:\Users\ssari\Desktop\ORCA_cube_generator\mep.py", line 130, in
vpot = read_vpot(f"{basename}_mep.out")
File "C:\Users\ssari\Desktop\ORCA_cube_generator\mep.py", line 38, in read_vpot
with open(vpot) as fp:
FileNotFoundError: [Errno 2] No such file or directory: 'test_mep.out'

If I copy the densities file and change the name to test.scfp then I get the following error:

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

Reading the GBW file ... done
Reading the positions ... done (512000 positions)
[file orca_tools/qcmat1.cpp, line 103]: Stored DataType 4616 does not match expected 8 in OrcaObject_ReadHeader!

[file orca_tools/qcmat1.cpp, line 103]: Stored DataType 4616 does not match expected 8 in OrcaObject_ReadHeader!

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

What am I doing wrong ?

@ssarix You should not have to copy any density file. Start with the water example; the input is posted above. First, run the Orca calculation using orca water.inp; after this finishes, run python mep.py water 80. Once this works for you, debugging other calculations will be easier.

@ssarix You should not have to copy any density file. Start with the water example; the input is posted above. First, run the Orca calculation using orca water.inp; after this finishes, run python mep.py water 80. Once this works for you, debugging other calculations will be easier.

Thank you for your answer. I did what you suggested with the water and it worked. I did an ORCA calculation with my molecule with the line:
! B3LYP 6-31G** OPT FREQ PAL6 VeryTightSCF LARGEPRINT keepdens xyzfile
The calculation went smoothly, but after that when I tried to run the script with the outputs I get the following error:

---

                  ORCA ELECTROSTATIC POTENTIAL GENERATION

---

GBW file ... 10a.gbw
File with points to evaluate V(r) ... 10a_mep.xyz
Output file ... 10a_mep.out

Reading the GBW file ... done
Reading the positions ... done (512000 positions)
Error (ORCA_VPOT): Cannot read density matrix of dimension 504 (10a.scfp)
re-compute Schwartz prescreening matrix in partically contracted basis ...done ( 0.100 sec)
Traceback (most recent call last):
File "C:\Users\ssari\Desktop\Orca_work\mep.py", line 130, in
vpot = read_vpot(f"{basename}_mep.out")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\ssari\Desktop\Orca_work\mep.py", line 38, in read_vpot
with open(vpot) as fp:
^^^^^^^^^^
FileNotFoundError: [Errno 2] No such file or directory: '10a_mep.out'

Is there something wrong with my basis set or functional ?

I had the same issue.
In the python script the line "subprocess.check_call([orca_vpot, f"{basename}_mep.inp"])" was the problem for me.
orca_vpot should be a string. Just put it in quotationmarks and it should work if orca_vpot is in your path.

Hi,

i have some problems. First orca 5 do not generate am *.scfp file. And the resulting shape is monochom ( tested with chimerax). In jmol and avogadro i can see only the molecule while opening the cube file. Is there a solution? The image show this for Fluoroethane in chimeraX

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

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