Source code for pyprocar.scriptFermi3D

# -*- coding: utf-8 -*-

import numpy as np
import pyvista

from matplotlib import colors as mpcolors
from matplotlib import cm
from .core.surface import boolean_add
from .fermisurface3d import FermiSurface3D
from .splash import welcome
from .utilsprocar import UtilsProcar
from .procarparser import ProcarParser
from .procarselect import ProcarSelect
from .bxsfparser import BxsfParser
from .frmsfparser import FrmsfParser
from .qeparser import QEFermiParser
from .lobsterparser import LobsterFermiParser


[docs]def fermi3D( procar="PROCAR", outcar="OUTCAR", infile="in.bxsf", fermi=None, bands=None, interpolation_factor=1, mode="plain", supercell=[1, 1, 1], colors=None, background_color="white", save_colors=False, cmap="jet", atoms=None, orbitals=None, spin=None, spin_texture=False, arrow_color=None, arrow_size=0.015, only_spin=False, fermi_shift=0, projection_accuracy="normal", code="vasp", vmin=0, vmax=1, savegif=None, savemp4=None, save3d=None, perspective=True, save2d=False, camera_pos=[1, 1, 1], widget=False, show=True, repair=True, ): """ Parameters ---------- procar : str, optional (default ``'PROCAR'``) Path to the PROCAR file of the simulation e.g. ``procar='~/MgB2/fermi/PROCAR'`` outcar : str, optional (default ``'OUTCAR'``) Path to the OUTCAR file of the simulation e.g. ``outcar='~/MgB2/fermi/OUTCAR'`` infile : str, optional (default ``infile = in.bxsf'``) This is the path in the input bxsf file e.g. ``infile = ni_fs.bxsf'`` fermi : float, optional (default ``None``) Fermi energy at which the fermi surface is created. In other words fermi is the isovalue at which the fermi surface is created. If not defined it is read from the OUTCAR file. e.g. ``fermi=-5.49`` bands : list int, optional Which bands are going to be plotted in the fermi surface. The numbering is based on vasp outputs. If nothing is selected, this function will iterate over all the bands and plots the ones that cross fermi. e.g. ``bands=[14, 15, 16, 17]`` interpolation_factor : int, optional The kpoints grid will increase by this factor and interpolated at the new points using Fourier interpolation. e.g. If the kgrid is 5x5x5, ``interpolation_factor=4`` will lead to a kgrid of 20x20x20 mode : str, optional (default ``mode='plain'``) Defines If the fermi surface will have any projection using colormaps or is a plotted with a uniform plain color. e.g. ``mode='plain'``, ``mode='parametric'`` supercell : list int, optional (default ``[1, 1, 1]``) If one wants plot more than the 1st brillouin zone, this parameter can be used. e.g. ``supercell=[2, 2, 2]`` colors : list str, optional List of colors for each band. This argument does not work when a 3d file is saved. The colors for when ``save3d`` is used, we recomend using qualitative colormaps, as this function will automatically choose colors from the colormaps. e.g. ``colors=['red', 'blue', 'green']`` background_color : str, optional (default ``white``) Defines the background color. e.g. ``background_color='gray'`` save_colors : bool, optional (default ``False``) In case the plot is saved in 3D and some property of the material is projected on the fermi surface, this argument allows the projection to be stored in the 3D file. e.g. ``save_colors=True`` cmap : str, optional (default ``jet``) The color map used for color coding the projections. ``cmap`` is only relevant in ``mode='parametric'``. A full list of color maps in matplotlib are provided in this web page. `https://matplotlib.org/2.0.1/users/colormaps.html <https://matplotlib.org/2.0.1/users/colormaps.html>`_ atoms : list int, optional ``atoms`` define the projection of the atoms on the fermi surfcae . In other words it selects only the contribution of the atoms provided. Atoms has to be a python list(or numpy array) containing the atom indices. Atom indices has to be order of the input files of DFT package. ``atoms`` is only relevant in ``mode='parametric'``. keep in mind that python counting starts from zero. e.g. for SrVO\ :sub:`3`\ we are choosing only the oxygen atoms. ``atoms=[2, 3, 4]``, keep in mind that python counting starts from zero, for a **POSCAR** similar to following:: Sr1 V1 O3 1.0 3.900891 0.000000 0.000000 0.000000 3.900891 0.000000 0.000000 0.000000 3.900891 Sr V O 1 1 3 direct 0.500000 0.500000 0.500000 Sr atom 0 0.000000 0.000000 0.000000 V atom 1 0.000000 0.500000 0.000000 O atom 2 0.000000 0.000000 0.500000 O atom 3 0.500000 0.000000 0.000000 O atom 4 if nothing is specified this parameter will consider all the atoms present. orbitals : list int, optional ``orbitals`` define the projection of orbitals on the fermi surface. In other words it selects only the contribution of the orbitals provided. Orbitals has to be a python list(or numpy array) containing the Orbital indices. Orbitals indices has to be order of the input files of DFT package. The following table represents the indecies for different orbitals in **VASP**. +-----+-----+----+----+-----+-----+-----+-----+-------+ | s | py | pz | px | dxy | dyz | dz2 | dxz | x2-y2 | +-----+-----+----+----+-----+-----+-----+-----+-------+ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | +-----+-----+----+----+-----+-----+-----+-----+-------+ ``orbitals`` is only relavent in ``mode='parametric'`` e.g. ``orbitals=[1,2,3]`` will only select the p orbitals while ``orbitals=[4,5,6,7,8]`` will select the d orbitals. If nothing is specified pyprocar will select all the present orbitals. spin : list int, optional e.g. ``spin=[0]`` spin_texture : bool, optional (default False) In non collinear calculation one can choose to plot the spin texture on the fermi surface. e.g. ``spin_texture=True`` arrow_color : str, optional Defines the color of the arrows when ``spin_texture=True``. The default will select the colors based on the color map specified. If arrow_color is selected, all arrows will have the same color. e.g. ``arrow_color='red'`` arrow_size : int, optional As the name suggests defines the size of the arrows, when spin texture is selected. e.g. ``arrow_size=3`` only_spin : bool, optional If ``only_spin=True`` is selected, the fermi surface is not plotted and only the spins in the spin texture is plotted. fermi_shift : float, optional This parameter is useful when one wants to plot the iso-surface above or belove the fermi level. e.g. ``fermi_shift=0.6`` projection_accuracy : str, optional (default ``'normal'``) Selected the accuracy of projected properties. ``'normal'`` and ``'high'`` are the only two options. ``'normal'`` uses the fast but rather inaccurate nearest neighbor interpolation, while ``'high'`` uses the more accurate linear interpolation for the projection of the properties. e.g. ``projection_accuracy='high'`` code : str, optional (default ``'vasp'``) The DFT code in which the calculation is performed with. Also, if you want to read a .bxsf file set code ="bxsf" e.g. ``code='vasp'`` vmin : float, optional The maximum value in the color bar. cmap is only relevant in ``mode='parametric'``. e.g. vmin=-1.0 vmax : float, optional The maximum value in the color bar. cmap is only relevant in ``mode='parametric'``. e.g. vmax=1.0 savegif : str, optional pyprocar can save the fermi surface in a gif format. ``savegif`` is the path to which the gif is saved. e.g. ``savegif='fermi.gif'`` or ``savegif='~/MgB2/fermi.gif'`` savemp4 : str, optional pyprocar can save the fermi surface in a mp4 video format. ``savemp4`` is the path to which the video is saved. e.g. ``savegif='fermi.mp4'`` or ``savegif='~/MgB2/fermi.mp4'`` save3d : str, optional pyprocar can save the fermi surface in a 3d file format. pyprocar uses the `trimesh <https://github.com/mikedh/trimesh>`_ to save the 3d file. trimesh can export files with the following formats STL, binary PLY, ASCII OFF, OBJ, GLTF/GLB 2.0, COLLADA. ``save3d`` is the path to which the file is saved. e.g. ``save3d='fermi.glb'`` perspective : bool, optional To create the illusion of depth, perspective is used in 2d graphics. One can turn this feature off by ``perspective=False`` e.g. ``perspective=False`` save2d : str, optional The fermi surface can be saved as a 2D image. This parameter turns this feature on and selects the path at which the file is going to be saved. e.g. ``save2d='fermi.png'`` camera_pos : list float, optional (default ``[1, 1, 1]``) This parameter defines the position of the camera where it is looking at the fermi surface. This feature is important when one chooses to use the ``save2d``, ``savegif`` or ``savemp4`` option. e.g. ``camera_pos=[0.5, 1, -1]`` widget : , optional .. todo:: show : bool, optional (default ``True``) If set to ``False`` it will not show the 3D plot. Returns ------- s : pyprocar surface object The whole fermi surface added bands surfaces : list pyprocar surface objects list of fermi surface of each band """ welcome() if code == "vasp" or code == "abinit": if repair: repairhandle = UtilsProcar() repairhandle.ProcarRepair(procar, procar) print("PROCAR repaired. Run with repair=False next time.") if show: p = pyvista.Plotter() if code == "vasp": outcarparser = UtilsProcar() if fermi is None: e_fermi = outcarparser.FermiOutcar(outcar) else: e_fermi = fermi reciprocal_lattice = outcarparser.RecLatOutcar(outcar) procarFile = ProcarParser() procarFile.readFile(procar, False) data = ProcarSelect(procarFile, deepCopy=True) elif code == "qe": procarFile = QEFermiParser() reciprocal_lattice = procarFile.reclat data = ProcarSelect(procarFile, deepCopy=True) if fermi is None: e_fermi = procarFile.fermi else: e_fermi = fermi elif code == "lobster": procarFile = LobsterFermiParser() reciprocal_lattice = procarFile.reclat data = ProcarSelect(procarFile, deepCopy=True) if fermi is None: e_fermi = 0 else: e_fermi = fermi elif code == "bxsf": e_fermi = fermi procarFile = BxsfParser(infile= infile) reciprocal_lattice = procarFile.rec_lattice bands = np.arange(len(procarFile.bandEnergy[0, :])) elif code == "frmsf": e_fermi = fermi data = FrmsfParser(infile=infile) reciprocal_lattice = data.rec_lattice bands = np.arange(len(data.bands[0, :])) if bands is None: bands = np.arange(len(data.bands[0, :])) surfaces = [] spd = [] if mode == "parametric": if orbitals is None: orbitals = [-1] if atoms is None: atoms = [-1] if spin is None: spin = [0] data.selectIspin(spin) data.selectAtoms(atoms, fortran=False) data.selectOrbital(orbitals) for iband in bands: spd.append(data.spd[:, iband]) else: for iband in bands: spd.append(None) spd_spin = [] if spin_texture: dataX = ProcarSelect(procarFile, deepCopy=True) dataY = ProcarSelect(procarFile, deepCopy=True) dataZ = ProcarSelect(procarFile, deepCopy=True) dataX.kpoints = data.kpoints dataY.kpoints = data.kpoints dataZ.kpoints = data.kpoints dataX.spd = data.spd dataY.spd = data.spd dataZ.spd = data.spd dataX.bands = data.bands dataY.bands = data.bands dataZ.bands = data.bands dataX.selectIspin([1]) dataY.selectIspin([2]) dataZ.selectIspin([3]) dataX.selectAtoms(atoms, fortran=False) dataY.selectAtoms(atoms, fortran=False) dataZ.selectAtoms(atoms, fortran=False) dataX.selectOrbital(orbitals) dataY.selectOrbital(orbitals) dataZ.selectOrbital(orbitals) for iband in bands: spd_spin.append( [dataX.spd[:, iband], dataY.spd[:, iband], dataZ.spd[:, iband]] ) else: for iband in bands: spd_spin.append(None) counter = 0 for iband in bands: print("Trying to extract isosurface for band %d" % iband) if code == "bxsf": surface = FermiSurface3D( kpoints=procarFile.kpoints, band=procarFile.bandEnergy[:, iband], spd=spd[counter], spd_spin=spd_spin[counter], fermi=e_fermi + fermi_shift, reciprocal_lattice=reciprocal_lattice, interpolation_factor=interpolation_factor, projection_accuracy=projection_accuracy, supercell=supercell, file="bxsf", ) elif code == "qe": surface = FermiSurface3D( kpoints=procarFile.kpoints, band=procarFile.bands[:, iband], spd=spd[counter], spd_spin=spd_spin[counter], fermi=e_fermi + fermi_shift, reciprocal_lattice=reciprocal_lattice, interpolation_factor=interpolation_factor, projection_accuracy=projection_accuracy, supercell=supercell, file="qe", ) elif code == "lobster": surface = FermiSurface3D( kpoints=procarFile.kpoints, band=procarFile.bands[:, iband], spd=spd[counter], spd_spin=spd_spin[counter], fermi=e_fermi + fermi_shift, reciprocal_lattice=reciprocal_lattice, interpolation_factor=interpolation_factor, projection_accuracy=projection_accuracy, supercell=supercell, file="lobster", ) elif code == "frmsf": surface = FermiSurface3D( kpoints=data.kpoints, band=data.bands[:, iband], spd=spd[counter], spd_spin=spd_spin[counter], fermi=e_fermi + fermi_shift, reciprocal_lattice=reciprocal_lattice, interpolation_factor=interpolation_factor, projection_accuracy=projection_accuracy, supercell=supercell, file="bxsf", ) else: surface = FermiSurface3D( kpoints=data.kpoints, band=data.bands[:, iband], spd=spd[counter], spd_spin=spd_spin[counter], fermi=e_fermi + fermi_shift, reciprocal_lattice=reciprocal_lattice, interpolation_factor=interpolation_factor, projection_accuracy=projection_accuracy, supercell=supercell, ) if surface.verts is not None: surfaces.append(surface) counter += 1 nsurface = len(surfaces) norm = mpcolors.Normalize(vmin=vmin, vmax=vmax) cmap = cm.get_cmap(cmap) scalars = np.arange(nsurface + 1) / nsurface if save_colors is not None or save3d is not None: for i in range(nsurface): if surfaces[i].scalars is None: surfaces[i].set_scalars([scalars[i]] * surfaces[i].nfaces) if colors is None: colors = np.array([cmap(norm(x)) for x in (scalars)]).reshape(-1, 4) print(scalars) print(colors) print(cmap) if show or save2d: # sargs = dict(interactive=True) p.add_mesh( surfaces[0].brillouin_zone.pyvista_obj, style="wireframe", line_width=3.5, color="black", ) for isurface in range(nsurface): if not only_spin: if mode == "plain": p.add_mesh(surfaces[isurface].pyvista_obj, color=colors[isurface]) text = "Plain" elif mode == "parametric": p.add_mesh( surfaces[isurface].pyvista_obj, cmap=cmap, clim=[vmin, vmax] ) p.remove_scalar_bar() text = "Projection" else: text = "Spin Texture" if spin_texture: # Example dataset with normals # create a subset of arrows using the glyph filter arrows = surfaces[isurface].pyvista_obj.glyph( orient="vectors", factor=arrow_size ) if arrow_color is None: p.add_mesh(arrows, cmap=cmap, clim=[vmin, vmax]) p.remove_scalar_bar() else: p.add_mesh(arrows, color=arrow_color) if mode != "plain" or spin_texture: p.add_scalar_bar( title=text, n_labels=6, italic=False, bold=False, title_font_size=None, label_font_size=None, position_x=0.9, position_y=0.01, color="black", ) p.add_axes( xlabel="Kx", ylabel="Ky", zlabel="Kz", line_width=6, labels_off=False ) if not perspective: p.enable_parallel_projection() p.set_background(background_color) p.set_position(camera_pos) if not widget: p.show(cpos=camera_pos, screenshot=save2d) # p.screenshot('1.png') # p.save_graphic('1.pdf') if savegif is not None: path = p.generate_orbital_path(n_points=36) p.open_gif(savegif) p.orbit_on_path(path) # ,viewup=camera_pos) if savemp4: path = p.generate_orbital_path(n_points=36) p.open_movie(savemp4) p.orbit_on_path(path) # ,viewup=camera_pos) # p.close() s = boolean_add(surfaces) s.set_color_with_cmap(cmap=cmap, vmin=vmin, vmax=vmax) # s.pyvista_obj.plot() # s.trimesh_obj.show() if save3d is not None: extention = save3d.split(".")[-1] s.export(save3d, extention) return s, surfaces