-
-
Save mretegan/5501553 to your computer and use it in GitHub Desktop.
| """ | |
| (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 |
| # 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 | |
| * |
Thank you very much
Hi,
I need to use the mep.py for create a .cube file of the electrostatic potential using ORCA. unfortunately it dos't work for me. I will be grateful if any one can help me to solve the problem.
Problem: At first I run the mep.py with " python mep.py basename npoints"
$python mep.py input 80
Traceback (most recent call last):
File "./mep.py", line 98, in
basename + "_mep.inp", basename + "_mep.out"])
File "/usr/lib64/python2.7/subprocess.py", line 506, in check_call
retcode = call(*popenargs, **kwargs)
File "/usr/lib64/python2.7/subprocess.py", line 493, in call
return Popen(*popenargs, **kwargs).wait()
File "/usr/lib64/python2.7/subprocess.py", line 679, in init
errread, errwrite)
File "/usr/lib64/python2.7/subprocess.py", line 1249, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory
( that my python version was python2.7) so I installed the python 3 and run again the mep.py and get the below error.
$ python3 mep.py input 80
Traceback (most recent call last):
File "mep.py", line 78, in
atoms, x, y, z = read_xyz(basename + ".xyz")
File "mep.py", line 29, in read_xyz
f.next()
AttributeError: '_io.TextIOWrapper' object has no attribute 'next'
please let me know the solution.
please let me know the solution.
For me the following worked: add the full path to your orca_vpot program in line 92:
subprocess.check_call(["YOURPATH/orca_vpot", basename + ".gbw", basename + ".scfp", basename + "_mep.inp", basename + "_mep.out"])
Hello nice ppl! Really need some help here. I'm pretty ignorant w/ terminals and all...tried to follow the intructions, and here is what I get:
ina:3F-bzn ina$ python mep.py 3F-bzn 80
Traceback (most recent call last):
File "mep.py", line 74, in
atoms, x, y, z = read_xyz(basename + ".xyz")
File "mep.py", line 32, in read_xyz
z.append(float(data[3]))
IndexError: list index out of range
ina:3F-bzn ina$
Anyone!?
Hello nice ppl! Really need some help here. I'm pretty ignorant w/ terminals and all...tried to follow the intructions, and here is what I get:
ina:3F-bzn ina$ python mep.py 3F-bzn 80
Traceback (most recent call last):
File "mep.py", line 74, in
atoms, x, y, z = read_xyz(basename + ".xyz")
File "mep.py", line 32, in read_xyz
z.append(float(data[3]))
IndexError: list index out of range
ina:3F-bzn ina$Anyone!?
Something is wrong with your 3F-bzn.xyz file. Do you have it in the folder where you run the script?
With Retegan 2019 revision and Anaconda python3 error
ORCA ELECTROSTATIC POTENTIAL GENERATION
GBW file ... qmmm_0.input.gbw
File with points to evaluate V(r) ... qmmm_0.input_mep.inp
Output file ... qmmm_0.input_mep.out
Reading the GBW file ... done
Reading the positions ... done (512000 positions)
Traceback (most recent call last):
File "mep.py", line 112, in
basename + '_mep.inp', basename + '_mep.out'])
File "/opt/anaconda3/lib/python3.7/subprocess.py", line 363, in check_call
raise CalledProcessError(retcode, cmd)
subprocess.CalledProcessError: Command '['/home/francesco/tmp2/orca_vpot', 'qmmm_0.input.gbw', 'qmmm_0.input.scfp', 'qmmm_0.input_mep.inp', 'qmmm_0.input_mep.out']' died with <Signals.SIGSEGV: 11>.
francesco@vaio64:~/tmp2$
thanks for help
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?
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, runpython 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 ?
great stuff. ty