Source code for pyprocar.core.isosurface

from .surface import Surface
import numpy as np
from scipy import interpolate
from skimage import measure

__author__ = "Pedram Tavadze"
__maintainer__ = "Pedram Tavadze"
__email__ = "petavazohi@mail.wvu.edu"
__date__ = "March 31, 2020"


[docs]class Isosurface(Surface): def __init__( self, XYZ=None, V=None, isovalue=None, V_matrix=None, algorithm='lewiner', interpolation_factor=1, padding=None, transform_matrix=None, boundaries=None, file = None ): """ This class contains a surface that finds all the poins correcponding to the following equation V(X,Y,Z) = f Parameters ---------- XYZ : TYPE, list of lists of floats, (n,3) DESCRIPTION. a list of coordinates [[x1,y1,z1],[x2,y2,z2],...] corresponding V V : TYPE, list of floats, (n,) DESCRIPTION. a list of values [V1,V2,...] corresponding to XYZ XYZ[0] >>> V[0] XYZ[1] >>> V[1] isovalue : TYPE, float DESCRIPTION. The constant value of the surface (f) V_matrix : TYPE, float (nx,ny,nz) DESCRIPTION. one can present V_matrix instead of XYZ and V. V_matrix is a matrix representation of XYZ and V together. This matrix is generated if XYZ and V are provided. algorithm : TYPE, string DESCRIPTION. The default is 'lewiner'. The algorithm used to find the isosurface, This function used scikit-image to find the isosurface. possibilities ['classic','lewiner'] interpolation_factor : TYPE, int DESCRIPTION. The default is 1. This module uses Fourier Transform interpolation. interpolation factor will increase the grid points in each direction by a this factor, the dafault is set to 1 padding : TYPE, list of float (3,) DESCRIPTION. padding is used for periodic datasets such as bands in a solid state calculation. e.g. The 1st BZ is not covered fully so one might want to pad the matrix with wrap(look at padding in numpy for wrap), afterwards one has to clip the surface to the first BZ. easily doable using pyvista of trimesh padding goes as follows np.pad(self.eigen_matrix, ((padding[0]/2, padding[0]/2), (padding[1]/2, padding[1]/2) (padding[2]/2, padding[2])), "wrap") In other words it creates a super cell withpadding transform_matrix : TYPE, (3,3) float DESCRIPTION. applies an transformation to the vertices VERTS_prime=T*VERTS boundaries : TYPE, pyprocar surface DESCRIPTION. The default is None. The boundaries in which the isosurface will be clipped with for example the first brillouin zone """ self.XYZ = np.array(XYZ) self.V = V self.isovalue = isovalue self.V_matrix = V_matrix self.algorithm = algorithm self.padding = padding self.interpolation_factor = interpolation_factor self.transform_matrix = transform_matrix self.boundaries = boundaries self.file = file if self.algorithm not in ['classic', 'lewiner']: print( "The algorithm chose has to be from ['classic','lewiner'], automtically choosing 'lewiner'" ) self.algorithm = 'lewiner' if self.V_matrix is None: self.V_matrix = map2matrix(self.XYZ, self.V) if self.padding is None: self.padding = [self.nX*2 // 2, self.nY*2 // 2, self.nZ*2 // 2] else: self.padding = [ self.nX // 2 * padding[0], self.nY // 2 * padding[1], self.nZ // 2 * padding[2] ] verts, faces, normals, values = self._get_isosurface( interpolation_factor) if verts is not None and faces is not None: if transform_matrix is not None: verts = np.dot(verts, transform_matrix) """ Python, unlike statically typed languages such as Java, allows complete freedom when calling methods during object initialization. However, standard object-oriented principles apply to Python classes using deep inheritance hierarchies. Therefore the developer has responsibility for ensuring that objects are properly initialized when there are multiple __init__ methods that need to be called. For this reason I will make one temporary surface and from there I will using the other surface provided. """ if boundaries is not None: suprecell_surface = Surface(verts=verts, faces=faces, face_normals=normals) if not np.isnan(suprecell_surface.pyvista_obj.points[0,0]): verts, faces = self.clip(suprecell_surface, boundaries) #verts, faces = self.clip(suprecell_surface, boundaries) Surface.__init__(self, verts=verts, faces=faces, face_normals=normals)
[docs] def clip(self, S1, S2): """ This function clips S1 using the boundaries of S2 and returns Parameters ---------- S1 : TYPE pyprocar surface DESCRIPTION. S2 : TYPE pyprocar surface DESCRIPTION. Returns ------- verts,faces """ for iface in range(len(S2.faces)): normal = S2.face_normals[iface] center = np.average(S2.verts[S2.faces[iface]], axis=0) S1.pyvista_obj.clip(origin=center, normal=normal, inplace=True) faces = [] courser = 0 for i in range(S1.pyvista_obj.n_faces): npoints = S1.pyvista_obj.faces[courser] face = [] courser += 1 for ipoint in range(npoints): face.append(S1.pyvista_obj.faces[courser]) courser += 1 faces.append(face) return S1.pyvista_obj.points, faces
@property def X(self): """ Returns ------- TYPE numpy array DESCRIPTION. list of grids in X direction """ return np.unique(self.XYZ[:, 0]) @property def Y(self): """ Returns ------- TYPE numpy array DESCRIPTION. list of grids in Y direction """ return np.unique(self.XYZ[:, 1]) @property def Z(self): """ Returns ------- TYPE numpy array DESCRIPTION. list of grids in Z direction """ return np.unique(self.XYZ[:, 2]) @property def dxyz(self): """ Returns ------- list DESCRIPTION. length between points in each direction """ dx = np.abs(self.X[-1] - self.X[-2]) dy = np.abs(self.Y[-1] - self.Y[-2]) dz = np.abs(self.Z[-1] - self.Z[-2]) return [dx, dy, dz] @property def nX(self): """ Returns ------- TYPE int DESCRIPTION. number of points in the grid in X direction """ return len(self.X) @property def nY(self): """ Returns ------- TYPE int DESCRIPTION. number of points in the grid in Y direction """ return len(self.Y) @property def nZ(self): """ Returns ------- TYPE int DESCRIPTION. number of points in the grid in Z direction """ return len(self.Z) @property def surface_boundaries(self): """ This function tries to find the isosurface using no interpolation to find the correct positions of the surface to be able to shift to the interpolated one to the correct position Returns ------- list of tuples DESCRIPTION. [(mins[0],maxs[0]),(mins[1],maxs[1]),(mins[2],maxs[2])] """ padding_x = self.padding[0] padding_y = self.padding[1] padding_z = self.padding[2] eigen_matrix = np.pad(self.V_matrix, ((padding_x, padding_x), (padding_y, padding_y), (padding_z, padding_z)), "wrap") try: verts, faces, normals, values = measure.marching_cubes_lewiner( eigen_matrix, self.isovalue) for ix in range(3): verts[:, ix] -= verts[:, ix].min() verts[:, ix] -= (verts[:, ix].max() - verts[:, ix].min()) / 2 #+self.origin[ix] verts[:, ix] *= self.dxyz[ix] mins = verts.min(axis=0) maxs = verts.max(axis=0) return [(mins[0], maxs[0]), (mins[1], maxs[1]), (mins[2], maxs[2])] except: return None def _get_isosurface(self, interp_factor=1): """ Parameters ---------- interp_factor : TYPE, optional DESCRIPTION. The default is 1. Returns ------- TYPE DESCRIPTION. verts TYPE DESCRIPTION. faces TYPE DESCRIPTION. normals TYPE DESCRIPTION. vlues """ # Amount of kpoints needed to add on to fully sample 1st BZ padding_x = self.padding[0] padding_y = self.padding[1] padding_z = self.padding[2] eigen_matrix = np.pad(self.V_matrix, ((padding_x, padding_x), (padding_y, padding_y), (padding_z, padding_z)), "wrap") bnd = self.surface_boundaries if interp_factor != 1: # Fourier interpolate the mapped function E(x,y,z) eigen_matrix = fft_interpolate(eigen_matrix, interp_factor) # after the FFT we loose the center of the BZ, using numpy roll we # bring back the center of the BZ # eigen_matrix = np.roll(eigen_matrix,4 , # axis=[0, 1, 2]) try: verts, faces, normals, values = measure.marching_cubes_lewiner( eigen_matrix, self.isovalue) except BaseException: print("No isosurface for this band") return None, None, None, None # recenter for ix in range(3): if self.file == "bxsf" or self.file == "qe" or self.file == 'lobster': # verts[:, ix] -= verts[:, ix].min() verts[:, ix] *= self.dxyz[ix] / interp_factor verts[:, ix] -= self.supercell[ix] else: verts[:, ix] -= verts[:, ix].min() verts[:, ix] -= (verts[:, ix].max() - verts[:, ix].min()) / 2 verts[:, ix] *= self.dxyz[ix] / interp_factor if bnd is not None and interp_factor != 1: verts[:, ix] -= (verts[:, ix].min() - bnd[ix][0]) #+self.origin[ix] # verts[:, ix] *= self.dxyz[ix] / interp_factor # print(self.dxyz) # if self.file == "bxsf": # verts[:, ix] -= 0.5 # if bnd is not None and interp_factor != 1: # print((verts[:, ix].min() - bnd[ix][0])) # verts[:, ix] -= (verts[:, ix].min() - bnd[ix][0]) # if self.file == "bxsf": # verts[:, ix] -= 0.50 # x_shift = verts[:,0].min() - bnd[0] # y_shift = verts[:,1].min() - bnd[1] # z_shift = verts[:,2].min() - bnd[2] # transfare from fraction to cartesian # verts = np.dot(verts, self.reciprocal_) # new_faces = np.zeros(shape=(len(faces), 4)) # new_faces[:, 0] = 3 # new_faces[:, 1:] = faces # faces = new_faces return verts, faces, normals, values
[docs]def map2matrix(XYZ, V): """ mapps an Irregular grid to a regular grid Parameters ---------- XYZ : TYPE DESCRIPTION. V : TYPE DESCRIPTION. Returns ------- mapped_func : TYPE DESCRIPTION. """ XYZ = XYZ V = V X = np.unique(XYZ[:, 0]) Y = np.unique(XYZ[:, 1]) Z = np.unique(XYZ[:, 2]) mapped_func = np.zeros(shape=(len(X), len(Y), len(Z))) #kpoint_matrix = np.zeros(shape=(len(kx), len(ky), len(kz), 3)) This was added to check if the mesh grid is working count = 0 for ix in range(len(X)): condition1 = XYZ[:, 0] == X[ix] count += 1 for iy in range(len(Y)): condition2 = XYZ[:, 1] == Y[iy] #print(count) for iz in range(len(Z)): condition3 = XYZ[:, 2] == Z[iz] tot_cond = np.all([condition1, condition2, condition3], axis=0) if len(V[tot_cond]) != 0: mapped_func[ix, iy, iz] = V[tot_cond][0] # kpoint_matrix[ikx, iky, ikz] = [ # kx[ikx], ky[iky], kz[ikz]] else: mapped_func[ix, iy, iz] = np.nan # kpoint_matrix[ikx, iky, ikz] = [np.nan, np.nan, np.nan] return mapped_func
[docs]def fft_interpolate(function, interpolation_factor=2): """ if I = interpolation_factor This function withh recieve f(x,y,z) with dimensions of (nx,ny,nz) and returns f(x,y,z) with dimensions of (nx*I,ny*I,nz*I) """ eigen_fft = np.fft.fftn(function) shifted_fft = np.fft.fftshift(eigen_fft) nx, ny, nz = np.array(shifted_fft.shape) pad_x = nx * (interpolation_factor - 1) // 2 pad_y = ny * (interpolation_factor - 1) // 2 pad_z = nz * (interpolation_factor - 1) // 2 new_matrix = np.pad( shifted_fft, ((pad_x, pad_x), (pad_y, pad_y), (pad_z, pad_z)), "constant", constant_values=0, ) new_matrix = np.fft.ifftshift(new_matrix) interpolated = np.real(np.fft.ifftn(new_matrix)) * (interpolation_factor** 3) return interpolated