"""
ePSproc orbPlot
Functions for calculating and plotting orbitals from quantum chemistry results.
Makes use of pyvista/itkwidgets, cclib and chemview/pyquante on the backend (see notes below).
13/05/20 v1, adapting from test notebook, cclib_chemlab_ITK_pyVista_dev_150420.ipynb
"""
# Required imports
import numpy as np
import sys
from pathlib import Path
import pyvista as pv
import cclib
from functools import reduce # For chemlab, qc.pgbf
# Import here doesn't fix issue, need to confirm scope.
# Adding "from functools import reduce" to chemlab.qc.pgbf.py does work.
# Set chemlab fn. import
[docs]def importChemLabQC(chemPath = None):
"""
Import chemlab. This is set as optional, and a local path can be defiend, since it is deprecated and may not install.
The chemlab.qc.wavefunction() method essentially repackages pyQuante functionality, for 3D volume calculation for orbitals defined by Gaussian basis sets.
Code: https://github.com/chemlab/chemlab/blob/master/chemlab/qc/wavefunction.py
Method: https://gabrielelanaro.github.io/blog/2015/01/14/visualizing-molecular-orbitals-ipython-notebook.html
To do: add more error parsing and logging, see, e.g., https://airbrake.io/blog/python-exception-handling/importerror-and-modulenotfounderror
"""
qcImport = False
try:
# import chemlab
from chemlab.qc import wavefunction
qcImport = True
except ImportError:
pass
if chemPath is not None:
try:
# Load dev scripts - sidesteps issues with chemlab installation
sys.path.append(chemPath)
# from chemlab.qc import wavefunction # Can't seem to get this version to work - likely due to init files?
from qc import wavefunction # This is OK with Base or subdir set
# from qc.utils import reduce # For additional function settings (or may need to patch?)
# This fix doesn't work - need to confirm scope here?
qcImport = True
except ImportError:
pass
if not qcImport:
print('Chemlab orbital functionality not found, orbital calculations not possible.')
return None
else:
print(f'Import OK: Chemlab module {wavefunction}')
return wavefunction
[docs]class molOrbPlotter():
"""
Class to define & calculate electronic structure properties & orbitals.
Parameters
----------
chemPath : str, optional
Set import path for chemlib. If not set, only standard python dir tree will be used.
fileIn : str, Path or jobInfo dict
Set file to use, either directly via str/Path, or use file from jobInfo dictionary (as returned by :py:func:`epsproc.headerFileParse`)
localPath : str, optional
Set local path to use for electronic structure file, instead of original path from jobInfo.
fileType : str, optional
Set to override electronic structure file type by suffix.
Usually this will be a .molden file for ePS, with a corresponding original file (e.g. Gamess .log).
CClib supports most formats, although doesn't read Molden. See https://cclib.github.io/data.html
Returns
-------
molOrbPlotter class object, with molecular data + plotting methods.
Notes
------
To do
-----
- File serach logic? Set default case to call :py:func:`epsproc.headerFileParse`? Other locations?
Electronic structure
--------------------
Read files and process with cclib.
https://cclib.github.io/how_to_parse.html
Orbitals
--------
Method from chemview.viewer.MoleculeViewer.add_isosurfaces
The chemlab.qc.wavefunction() method essentially repackages pyQuante functionality, for 3D volume calculation for orbitals defined by Gaussian basis sets.
"""
def __init__(self, chemPath = None, fileIn = None, localPath = None, fileType = None, geomInd = -1):
"""Init molOrbPlotter. Read electronic structure file & set-up grids."""
self.wavefunction = importChemLabQC(chemPath = chemPath)
if type(fileIn) is dict:
# Get file name from jobInfo
self.fileIn = Path(fileIn['Convert'][0].split()[1].strip("'"))
# Override path if local path specified
if localPath is not None:
self.fileIn = Path(localPath, self.fileIn.name)
else:
self.fileIn = Path(fileIn)
if fileType is not None:
self.fileIn.suffix
if self.fileIn.is_file():
print(f'Found electronic structure file: {self.fileIn}')
self.data = cclib.io.ccread(self.fileIn.as_posix())
# print(self.data.metadata)
print("Read %i atoms and %i MOs" % (self.data.natom, self.data.nmo))
# Set which molecular geometry to use.
self.geomInd = geomInd
# Set dict for storing calc fns.
self.f = {}
# Set grid
self.setGrid()
# Set dictionary to hold volume calcs.
# self.vol = {}
# Set pyVista object to hold orbital calculations. NO SET IN self.setGrid()
# Note ordering!
# self.pvObj = pv.StructuredGrid(self.X, self.Z, self.Y) # OK
print('*** Grids set OK')
print(self.vol)
else:
print('***Missing electronic structure file: {self.filleIn}')
[docs] def setGrid(self, extent = None, resolution = 50, padGrid = 1.0, sqGrid = True, geomInd = None):
"""Setup grids for orbital calculations.
Basic 3D grid setup, based on passed or atomic coords.
Currently limited to square grid for compatibility with ITK/VTK plotting - but should be able to make this more flexible.
Final grid set as PyVista object, self.vol
Parameters
----------
extent : list or np.array, [min, max], optional
Define limits of grid (in working units, Angs or Bohr).
If not specified will be defined from atomic coords (min and max values in any dimension).
resolution : int, optional, default = 50
Number of points per dimension, final grid size will be resolution^3
padGrid : float, optional, default = 1.0
Padding to add when using atomic coords to define grid (in working units angs/bohr).
Currently single valued, same for all dims.
sqGrid : bool, optional, default = True
Force square gridding? Currently only True is supported.
geomInd : int, optional, default = -1
Set to select a specific geometry from the electronic structure file (if present).
Defaults to the final geometry listed.
Overrides self.geomInd if set.
Notes
------
- Method originally adapted from chemview.viewer.MoleculeViewer.add_isosurfaces
- Q: do units matter here...?
"""
#*** Check global vals
if geomInd is not None:
self.geomInd = geomInd
#*** Check and set grid limits
if extent is None:
# This will generate a [2x3] array of (min,max) values for (x,y,z) dims.
extent = np.array([self.data.atomcoords[self.geomInd].min(axis=0), self.data.atomcoords[self.geomInd].max(axis=0)])
else:
extent = np.asarray(extent) # Ensure np.array type.
padGrid = False # Skip padding option if grid preset.
# Check for passed size - can be single or all dims.
# Set to all dims if required.
if extent.size < 3:
extent = np.c_[np.ones(3)*extent[0], np.ones(3)*extent[1]].T
# If set, force grid to square - same min and max for all dims
if sqGrid:
extent[0,:] = extent.min()
extent[1,:] = extent.max()
# Add grid padding
if padGrid:
extent[0,:] -= padGrid
extent[1,:] += padGrid
#*** Set axes
# TODO: better storage options here - Xarray or HV?
# See also efield.py for more.
axes = []
for col in extent.T:
# print(col)
axes.append(np.linspace(col[0], col[1], resolution))
self.axes = np.array(axes) # Set as axes here to avoid confusion with pyVista coords later.
# self.x = axes[0]
# self.y = axes[1]
# self.z = axes[2]
#*** Set grid
# self.X, self.Y, self.Z = np.meshgrid(axes[0], axes[1], axes[2]) # Store to self
X, Y, Z = np.meshgrid(axes[0], axes[1], axes[2]) # Local values only, used to set pyVista object
# spacing = np.array((area_max - area_min)/resolution)
# Set pyVista object to hold orbital calculations.
# Note ordering!
self.vol = pv.StructuredGrid(X, Z, Y)
[docs] def calcOrb(self, orbN = None):
"""
Calculate orbital volume.
Parameters
----------
orbN : int, optional, default = None
Set orbital number to plot.
This will default to the homo (as defined by cclib data.homos[-1]) if not set.
Note there is a -1 indexing offset for zero-indexed python arrays, which is **applied** in the function.
"""
if orbN is None:
orbN = self.data.homos[-1]
else:
orbN -= 1 # Set to python indexing
# Set orbital coeffs.
# For unrestricted spin case, try stacking, athough this may be incorrect!
if len(self.data.mocoeffs) == 1:
coeffs = self.data.mocoeffs[0][orbN]
else:
print('***Warning: Urestricted orbitals not yet tested!')
coeffs = np.array([self.data.mocoeffs[0][orbN], self.data.mocoeffs[1][orbN]])
# Set function & calculate volume
self.f[orbN+1] = self.wavefunction.molecular_orbital(self.data.atomcoords[self.geomInd], coeffs, self.data.gbasis) # Note first input is atom positions, also in data.atomcoords
# Calculate and assign as pyVista array, named as orbital N (note - must be str)
# self.vol[orbN+1] = self.f(self.x, self.y, self.z).flatten(order='K')
self.vol[str(orbN+1)] = self.f[orbN+1](self.vol.x, self.vol.y, self.vol.z).flatten(order='K')
# Store in dictionary (or Xarray, or other...?)
# AH - should just set pvObject here, and add each orbital as a data array.
# Now done at init, and set by orbN above.
[docs] def plotOrb(self, orbN = None, isoLevels = 6, isoValsAbs = None, isoValsPC = None, showAtoms = True, interactive = True, opacity = 0.5, subplots = False, rc = None, verbose = 0):
"""
Plot orbital(s) with pyVista/ITK
orbN : str or list of strs, optional, default = None
Orbital(s) to plot, by name in self.vol PV object.
By default these are set to orbital numbers, and exist for all calculated orbitals.
If not supplied, plot all available surfaces, i.e. items in self.vol.array_names
isoLevels : int, optional, default = 6
Number of isosurfs to compute & render.
isoValsAbs : list or array, default = None
Isovals to use (absolute values).
If both isoValsAbs and isoValsPC are passed, only PC values will be used.
isoValsPC : list or array, default = None
Isovals to use, percentage values. Plot values are set by +/-isoValsPC * np.abs(self.vol[item]).mean().
If both isoValsAbs and isoValsPC are passed, only PC values will be used.
showAtoms : bool, optional, default = True
If True, plot atoms (molecular structure).
interactive : bool, optional, default = True
Plot inteactively using itkwidgets if true.
Otherwise use pyVista.Plotter(), which produces static output.
opacity : float, optional, default = 0.5
Set isosurface opacity.
subplots : bool, optional, default = False
Set to one subplot per volume if True.
rc : array, optional, default = None
If set, and subplots is True, use to define layout grid [rows, columns].
If not set, this will be set for nCols = 3
To do
-----
- Fix naming & colourmapping for multiple objects in ITK plotter.
- Consolidate atoms > molecular geometry pyVista object.
- orbN as list of ints or np.array? Integrate with calcOrb()?
- Opacity mapping for isosurfs? For pv.Plotter() can do this via transfer fns, e.g. 'linear', https://docs.pyvista.org/examples/02-plot/opacity.html#sphx-glr-examples-02-plot-opacity-py
"""
# Add meshes per orbital, or use supplied list?
# NEEDS DEBUG!
if orbN is None:
orbN = self.vol.array_names
elif type(orbN) is str:
orbN = list(orbN)
# Set ITK widgets or pv.Plotter
# May also want to add other methods here, e.g. pv.BackgroundPlotter() # Runs plotter in separate window process.
if interactive:
pl = pv.PlotterITK()
else:
# Subplot code from wfPlot.plotWf()
# See also sphPlot.sphPlotPL()
if subplots:
# Set subplots if required.
# if subplots:
if rc is None:
pN = len(orbN)
rP = np.rint(np.sqrt(pN)).astype(np.int32)
cP = np.ceil(pN/rP).astype(np.int32)
sPlots=(rP,cP)
else:
sPlots=(rc[0],rc[1])
currPlot = [0, 0]
pl = pv.Plotter(shape = sPlots)
else:
pl = pv.Plotter()
# Add atoms
if showAtoms:
for n, item in enumerate(self.data.atomcoords[self.geomInd]):
# pl.add_mesh(pv.Sphere(center=item, radius=data.atomnos[n]/100)) # Relative scaling, small atoms
# pl.add_mesh(pv.Sphere(center=item, radius=np.sqrt(data.atomnos[n]/10))) # sqrt scaling, space filling atoms
pl.add_mesh(pv.Sphere(center=item, radius=np.sqrt(self.data.atomnos[n])/10), name = f'N{n}, Z={self.data.atomnos[n]}') # sqrt scaling, mid-sized atoms
# NOTE 18/07/20 - should be orbN list here? Didn't test yet.
# for item in self.vol.array_names:
# 26/10/20 - fixing list selection + added subplotting option
for item in orbN:
if item in self.vol.array_names:
if subplots:
pl.subplot(currPlot[0], currPlot[1])
# Increment
currPlot[1] += 1
if currPlot[1] > (sPlots[1]-1): # Wrap rows - probably a neater way to do this...!
currPlot[1] = 0
currPlot[0] += 1
limitVals = [self.vol[item].min(), self.vol[item].max(), np.abs(self.vol[item]).mean()]
# Set default case
isoValsOrb = np.linspace(-limitVals[2], limitVals[2], isoLevels) # Again set from mean.
# Override if alternative limits set (logic may be dodgy here)
if isoValsAbs is not None:
isoValsOrb = np.array(isoValsAbs)
if isoValsPC is not None:
isoValsPC = np.array(isoValsPC)
isoValsOrb = np.r_[isoValsPC, -isoValsPC] * limitVals[2] # Set PC vals from mean?
# Ensure ordering, may be required if using for clims
# isoValsOrb.sort()
# print(isoValsOrb)
# Add contours for currently selected scalars (orbital)
pl.add_mesh(self.vol.contour(isosurfaces = isoValsOrb, scalars = item), smooth_shading=True,
opacity=opacity,
clim=[isoValsOrb.min(), isoValsOrb.max()],
stitle=item)
# clim=[isoValsOrb[0], isoValsOrb[-1]]) # Requires sorted order
if subplots:
pl.add_scalar_bar(item) # For subplot case, this ensures cbar per subplot, and correct label, but THIS ALWAYS SEEMS TO set a shared cmap, WHY????
# Not the case in demos, e.g. https://docs.pyvista.org/examples/02-plot/multi-window.html
if verbose:
print(f"Volume: {item}, limitVals: {limitVals}, isoValsOrb: {isoValsOrb}")
else:
print(f"*** Volume {item} not found, skipping plot")
# Render plot
# In notebook tests this doesn't reneder unless called again?
self.pl = pl
# Additional options (for some plotters)
if hasattr(self.pl, 'add_axes'):
self.pl.add_axes()
self.pl.show()