Source code for chemtools.utils.cube

# -*- coding: utf-8 -*-
# ChemTools is a collection of interpretive chemical tools for
# analyzing outputs of the quantum chemistry calculations.
#
# Copyright (C) 2016-2019 The ChemTools Development Team
#
# This file is part of ChemTools.
#
# ChemTools is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# ChemTools is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see <http://www.gnu.org/licenses/>
#
# --
"""The Cube Module."""


import logging
import numpy as np

from chemtools.wrappers.molecule import Molecule

try:
    from importlib_resources import path
except ImportError:
    from importlib.resources import path


__all__ = ['UniformGrid']


[docs]class UniformGrid(object): """Class for generating a cubic grid and writing cube files.""" def __init__(self, numbers, pseudo_numbers, coordinates, origin, axes, shape): """Initialize ``UniformGrid`` class based on the origin, axes and shape of the cube. Parameters ---------- numbers : np.ndarray, shape=(M,) Atomic number of `M` atoms in the molecule. pseudo_numbers : np.ndarray, shape=(M,) Pseudo-number of `M` atoms in the molecule. coordinates : np.ndarray, shape=(M, 3) Cartesian coordinates of `M` atoms in the molecule. origin : np.ndarray, shape=(3,) Cartesian coordinates of the cubic grid origin. axes : np.ndarray, shape=(3, 3) The three vectors, stored as rows of axes array, defining the Cartesian coordinate system used to build the cubic grid. shape : np.ndarray, shape=(3,) Number of grid points along `x`, `y`, and `z` axis. """ self._numbers = numbers self._pseudo_numbers = pseudo_numbers self._coordinates = coordinates if origin.shape[0] != 3: raise ValueError('Argument origin should be an np.ndarray with shape=(3,)') self._origin = origin if axes.shape[0] != 3 or axes.shape[1] != 3: raise ValueError('Argument axes should be an np.ndarray with shape=(3, 3)') self._axes = axes if shape.shape[0] != 3: raise ValueError('Argument shape should be an np.ndarray with shape=(3,)') self._shape = shape # # Make cubic grid # # Number of points along x, y and z axis npoints_x, npoints_y, npoints_z = self._shape # Total number of grid points self._npoints = npoints_x * npoints_y * npoints_z # Make an array to store coordinates of grid points self._points = np.zeros((self._npoints, 3)) coords = np.array( np.meshgrid(np.arange(npoints_x), np.arange(npoints_y), np.arange(npoints_z)) ) coords = np.swapaxes(coords, 1, 2) coords = coords.reshape(3, -1) coords = coords.T self._points = coords.dot(self._axes) # Compute coordinates of grid points relative to the origin self._points += self._origin # log information self._log_init()
[docs] @classmethod def from_molecule(cls, molecule, spacing=0.2, extension=5.0, rotate=True): """Initialize ``UniformGrid`` class from Molecule object. Parameters ---------- molecule: instance of `Molecule` Instance of Molecule class. spacing : float, optional Increment between grid points along `x`, `y` and `z` direction. extension : float, optional The extension of the cube on each side of the molecule. rotate : bool, optional When True, the molecule is rotated so the axes of the cube file are aligned with the principle axes of rotation of the molecule. """ numbers = molecule.numbers pseudo_numbers = molecule.pseudo_numbers coordinates = molecule.coordinates # calculate center of mass of the nuclear charges: totz = np.sum(pseudo_numbers) com = np.dot(pseudo_numbers, coordinates) / totz if rotate: # calculate moment of inertia tensor: itensor = np.zeros([3, 3]) for i in range(pseudo_numbers.shape[0]): xyz = coordinates[i] - com r = np.linalg.norm(xyz)**2.0 tempitens = np.diag([r, r, r]) tempitens -= np.outer(xyz.T, xyz) itensor += pseudo_numbers[i] * tempitens _, v = np.linalg.eigh(itensor) new_coordinates = np.dot((coordinates - com), v) axes = spacing * v else: # Just use the original coordinates new_coordinates = coordinates # Compute the unit vectors of the cubic grid's coordinate system axes = np.diag([spacing, spacing, spacing]) # maximum and minimum value of x, y and z coordinates max_coordinate = np.amax(new_coordinates, axis=0) min_coordinate = np.amin(new_coordinates, axis=0) # Compute the required number of points along x, y, and z axis shape = (max_coordinate - min_coordinate + 2.0 * extension) / spacing shape = np.ceil(shape) shape = np.array(shape, int) # Compute origin origin = com - np.dot((0.5 * shape), axes) return cls(numbers, pseudo_numbers, coordinates, origin, axes, shape)
[docs] @classmethod def from_cube(cls, fname): r"""Initialize ``UniformGrid`` class based on the grid specifications of a cube file. Parameters ---------- fname : str Cube file name with \*.cube extension. """ fname = str(fname) if not fname.endswith('.cube'): raise ValueError('Argument fname should be a cube file with *.cube extension!') # Extract the specifications of the cubic grid from cube file's header numbers, pseudo_numbers, coordinates, origin, axes, shape = cls._read_cube_header(fname) return cls(numbers, pseudo_numbers, coordinates, origin, axes, shape)
[docs] @classmethod def from_file(cls, fname, spacing=0.2, extension=5.0, rotate=True): """ Initialize ``UniformGrid`` class based on the grid specifications of a file. Parameters ---------- fname : str Path to molecule's file. spacing : float, optional Increment between grid points along `x`, `y` and `z` direction. extension : float, optional The extension of the cube on each side of the molecule. rotate : bool, optional When True, the molecule is rotated so the axes of the cube file are aligned with the principle axes of rotation of the molecule. """ # Load file logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') try: mol = Molecule.from_file(str(fname)) except IOError as _: try: with path('chemtools.data.examples', str(fname)) as fname: logging.info('Loading {0}'.format(str(fname))) mol = Molecule.from_file(str(fname)) except IOError as error: logging.info(error) return cls.from_molecule(mol, spacing, extension, rotate)
@property def numbers(self): """Atomic number of the atoms in the molecule.""" return self._numbers @property def pseudo_numbers(self): """Pseudo-number of the atoms in the molecule.""" return self._pseudo_numbers @property def coordinates(self): """Cartesian coordinates of the atoms in the molecule.""" return self._coordinates @property def centers(self): """Cartesian coordinates of the atoms in the molecule.""" return self._coordinates @property def origin(self): """Cartesian coordinate of the cubic grid origin.""" return self._origin @property def axes(self): """Array with axes of the cube. The three vectors, stored as rows of axes array, defining the Cartesian coordinate system used to build the cubic grid. """ return self._axes @property def shape(self): """Number of grid points along `x`, `y`, and `z` axis.""" return self._shape @property def npoints(self): """Total number of grid points.""" return self._npoints @property def points(self): """Cartesian coordinates of the cubic grid points.""" return self._points def _log_init(self): """Log an overview of the cube's properties.""" logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') logging.info("Initialized cube: {0}".format(self.__class__)) logging.info("Origin : {0}".format(self._origin)) logging.info("Axes 1 : {0}".format(self._axes[0])) logging.info("Axes 2 : {0}".format(self._axes[1])) logging.info("Axes 3 : {0}".format(self._axes[2])) logging.info("Shape : {0}".format(self._shape))
[docs] def generate_cube(self, fname, data): r"""Write the data evaluated on grid points into a cube file. Parameters ---------- fname : str Cube file name with \*.cube extension. data : np.ndarray, shape=(npoints,) An array containing the evaluated scalar property on the grid points. """ if not fname.endswith('.cube'): raise ValueError('Argument fname should be a cube file with `*.cube` extension!') if data.size != self._npoints: raise ValueError('Argument data should have the same size as the grid. ' + '{0}!={1}'.format(data.size, self._npoints)) # Write data into the cube file with open(fname, 'w') as f: # writing the cube header: f.write('Cubefile created with HORTON CHEMTOOLS\n') f.write('OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n') natom = len(self._numbers) x, y, z = self._origin f.write('{0:5d} {1:11.6f} {2:11.6f} {3:11.6f}\n'.format(natom, x, y, z)) rvecs = self._axes for i, (x, y, z) in zip(self._shape, rvecs): f.write('{0:5d} {1:11.6f} {2:11.6f} {3:11.6f}\n'.format(i, x, y, z)) for i, q, (x, y, z) in zip(self._numbers, self._pseudo_numbers, self._coordinates): f.write('{0:5d} {1:11.6f} {2:11.6f} {3:11.6f} {4:11.6f}\n'.format(i, q, x, y, z)) # writing the cube data: num_chunks = 6 for i in range(0, data.size, num_chunks): row_data = data.flat[i:i+num_chunks] f.write((row_data.size*' {:12.5E}').format(*row_data)) f.write('\n')
[docs] def weights(self, method='R'): """ Return integration weights at every point on the cubic grid. Parameters ---------- method : str, optional The method for computing the integration weights at every point on the grid. Options: - 'R' method perfors rectangle/trapezoidal rule, without assuming that the function is close to zero at the edges of the grid. - 'R0' method performing rectangle/trapezoidal rule, assuming that the function is very close to zero at the edges of the grid. """ if method == 'R': volume = np.linalg.norm(self._shape[0] * self._axes[0]) volume *= np.linalg.norm(self._shape[1] * self._axes[1]) volume *= np.linalg.norm(self._shape[2] * self._axes[2]) numpnt = 1.0 * self._npoints weights = np.full(self._npoints, volume / numpnt) elif method == 'R0': volume = np.linalg.norm((self._shape[0] + 1.0) * self._axes[0]) volume *= np.linalg.norm((self._shape[1] + 1.0) * self._axes[1]) volume *= np.linalg.norm((self._shape[2] + 1.0) * self._axes[2]) numpnt = (self._shape[0] + 1.0) * (self._shape[1] + 1.0) * (self._shape[2] + 1.0) weights = np.full(self._npoints, volume / numpnt) else: raise ValueError('Argument method {0} is not known.'.format(method)) return weights
[docs] def integrate(self, data, method='R0'): """ Integrate the data on a cubic grid. Parameters ---------- data : np.ndarray, shape=(npoints, m) Data at every point on the grid given as an array. The size of axis=0 of this array should equal the number of grid points. method : str, default='R0' The method for computing the integration weights at every point on the grid. Options: - 'R' method perfors rectangle/trapezoidal rule, without assuming that the function is close to zero at the edges of the grid. - 'R0' method performing rectangle/trapezoidal rule, assuming that the function is very close to zero at the edges of the grid. """ if data.shape[0] != self._npoints: raise ValueError('Argument data should have the same size as the grid for axis=0. ' + '{0}!={1}'.format(data.shape[0], self._npoints)) value = np.tensordot(self.weights(method=method), data, axes=(0, 0)) return value
@staticmethod def _read_cube_header(fname): """ Return specifications of the cubic grid from the given cube file. Parameters ---------- fname : str Cube file name with *.cube extension. """ with open(fname) as f: # skip the title f.readline() # skip the second line f.readline() def read_grid_line(line): """Read a number and (x, y, z) coordinate from the cube file line.""" words = line.split() return ( int(words[0]), np.array([float(words[1]), float(words[2]), float(words[3])], float) # all coordinates in a cube file are in atomic units ) # number of atoms and origin of the grid natom, origin = read_grid_line(f.readline()) # numer of grid points in A direction and step vector A, and so on shape0, axis0 = read_grid_line(f.readline()) shape1, axis1 = read_grid_line(f.readline()) shape2, axis2 = read_grid_line(f.readline()) shape = np.array([shape0, shape1, shape2], int) axes = np.array([axis0, axis1, axis2]) def read_coordinate_line(line): """Read atomic number and (x, y, z) coordinate from the cube file line.""" words = line.split() return ( int(words[0]), float(words[1]), np.array([float(words[2]), float(words[3]), float(words[4])], float) # all coordinates in a cube file are in atomic units ) numbers = np.zeros(natom, int) pseudo_numbers = np.zeros(natom, float) coordinates = np.zeros((natom, 3), float) for i in xrange(natom): numbers[i], pseudo_numbers[i], coordinates[i] = read_coordinate_line(f.readline()) # If the pseudo_number field is zero, we assume that no effective core # potentials were used. if pseudo_numbers[i] == 0.0: pseudo_numbers[i] = numbers[i] return numbers, pseudo_numbers, coordinates, origin, axes, shape