-
-
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. 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 |
| # 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 | |
| * |
@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:
- Make sure the water example works.
- 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
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