Source code for chemtools.utils.mesh

# -*- 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/>
#
# --
"""Module for constructing different mesh (grid coordinates)."""


import numpy as np


[docs]def mesh_plane(points, spacing, extension): """Return the grid points on the plane spanned by the given three coordinates. Parameters ---------- points : np.ndarray(3, 3) Points that are on the plane. Rows corespond to the different points, columns correspond to the x, y, and z components. The first set of coordinate will be used to establish the new vertical (i.e. y) direction in the plane. spacing : float Upper bound to the spacing between adjacent grid points. This means that the spacing between adjacent points will be less than the given value. extension : float Distance from the center of the three points that will define the edges of the mesh. Returns ------- grid_points : np.ndarray(N_v, N_h, 3) Points on the plane spanned by the given three points. First index corresponds to the position along the vertical direction of the plane mesh. This direction is the same as the direction of the first coordinate wrt origin. `N_v` is the number of points vertically along the the mesh. Second index corresponds to the position along the horizontal direction of the plane mesh. `N_v` is the number of points horizontally along the the mesh. Raises ------ TypeError If `plane_points` is not two-dimensional numpy array of shape (3, 3). If `spacing` is not a float. If `extension` is not int or float. ValueError If `spacing` is less than or equal to zero. If `extension` is less than or equal to zero. If the three points on the plane are on a line. """ if not (isinstance(points, np.ndarray) and points.ndim == 2 and points.shape == (3, 3)): raise TypeError("Arguments points must be given as a 2D numpy array of shape (3, 3).") if not isinstance(spacing, float) or spacing <= 0: raise TypeError("Argument spacing must be a float greater than 0.") if not isinstance(extension, (int, float)) or extension < 0: raise TypeError("Argument extension must be int or float greater than 0.") center = np.average(points, axis=0) vec = points - center length_vec = (np.sum(vec ** 2, axis=1) ** 0.5)[:, None] if np.any(length_vec == 0): raise ValueError("Three points on the plane cannot be in a line.") unit_vec = vec / length_vec if np.unique(unit_vec, axis=0).shape[0] < 3: raise ValueError("Three points on the plane cannot be in a line.") # edges of the cube edges = unit_vec * extension + points # length of the edges from the center length_edges = (np.sum((edges - center) ** 2, axis=1) ** 0.5)[:, None] # use first point as a reference angle_01 = np.arccos(unit_vec[0].dot(unit_vec[1])) angle_02 = np.arccos(unit_vec[0].dot(unit_vec[2])) height_1 = length_edges[1] * np.sin(angle_01 - np.pi / 2) width_1 = length_edges[1] * np.cos(angle_01 - np.pi / 2) height_2 = length_edges[2] * np.sin(angle_02 - np.pi / 2) width_2 = length_edges[2] * np.cos(angle_02 - np.pi / 2) height = max(height_1, height_2) + length_edges[0] width = width_1 + width_2 unit_vertical = unit_vec[0] origin = center + unit_vertical * length_edges[0] if height_1 > height_2: vec_horizontal = edges[1] - (center - height_1 * unit_vertical) else: vec_horizontal = edges[2] - (center - height_2 * unit_vertical) origin += vec_horizontal unit_horizontal = vec_horizontal / np.sum(vec_horizontal ** 2) ** 0.5 output = np.linspace( origin, origin - unit_horizontal * width, num=int(width / spacing) + 2, endpoint=True ) output = np.linspace( output, output - unit_vertical * height, num=int(height / spacing) + 2, endpoint=True ) return output