-
-
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 | |
* |
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 ?
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