Quasi-eigenstates ================= The quasi-eigenstates are defined as .. math:: |\Psi(E) \rangle &=\frac{1}{2\pi} \int_{-\infty}^{\infty} \mathrm{e}^{\mathrm{i}Et} |\psi(t) \rangle \mathrm{d}t \nonumber \\ &=\frac{1}{2\pi}\sum_i a_i \int_{-\infty}^{\infty} \mathrm{e}^{\mathrm{i}(E-E_i)t} |\phi_i \rangle \mathrm{d}t \nonumber \\ &=\sum_i a_i \delta(E-E_i)|\phi_i \rangle If the energy :math:`E` equals to an eigenvalue :math:`E_i`, then the quasi-eigenstate is the exact eigenstate corresponding to :math:`E_i`, or a superposition of degenerate eigenstates if :math:`E_i` is degenerate. In this tutorial, we will show how to calculate the quasi-eigenstates and compare them with the exact eigenstates. The script can be found at ``examples/advanced/quasi_eigen.py``. We begin with importing the necessary packages: .. code-block:: python import numpy as np import tbplas as tb And define the following functions to generate primitive cell and sample with a single vacancy: .. code-block:: python :linenos: def make_prim_cell() -> tb.PrimitiveCell: """ Make graphene primitive cell with a single vacancy. :return: primitive cell with vacancy """ prim_cell = tb.make_graphene_diamond() # Get the index of vacancy dim = (17, 17, 1) orb_map = tb.RegularOrbitalMap(dim, prim_cell.num_orb) vac_id = orb_map.pc2sc(8, 8, 0, 0) # Extend cell and remove orbital prim_cell = tb.extend_prim_cell(prim_cell, dim=dim) prim_cell.remove_orbitals([vac_id]) return prim_cell def make_sample() -> tb.Sample: """ Make graphene sample with a single vacancy. :return: sample with vacancy """ prim_cell = tb.make_graphene_diamond() # Get the index of vacancy dim = (17, 17, 1) orb_map = tb.RegularOrbitalMap(dim, prim_cell.num_orb) vac_id = orb_map.pc2sc(8, 8, 0, 0) # Create supercell and sample super_cell = tb.SuperCell(prim_cell, dim=dim, pbc=(True, True, False), vacancies=[vac_id]) sample = tb.Sample(super_cell) return sample We are going to create a :math:`17\times17\times1` graphene sample with a vacancy in the center. We utilize the :class:`.RegularOrbitalMap` class to get the index of the 0th orbital in the :math:`(8, 8, 0)` cell. For the :class:`.PrimitiveCell` class, vacancies are added by calling the ``remove_orbitals`` method, while for the :class:`.SuperCell` class they are added via the ``vacancies`` argument. Exact eigenstates ----------------- We define the following function to calculate the exact eigenstates: .. code-block:: python :linenos: def wfc_diag(prim_cell: tb.PrimitiveCell) -> None: """ Calculate wave function using exact diagonalization. :param prim_cell: primitive cell :return: None """ solver = tb.DiagSolver(prim_cell) solver.config.prefix = "diag" solver.config.k_points = np.array([[0.0, 0.0, 0.0]]) bands, states = solver.calc_states() i_b = prim_cell.num_orb // 2 wfc = np.abs(states[0, i_b])**2 wfc /= wfc.max() vis = tb.Visualizer() vis.plot_wfc(prim_cell, wfc, scatter=True, site_size=wfc * 100, site_color="r", with_model=True, model_style={"alpha": 0.25, "color": "gray"}, fig_name="diag.png", fig_size=(6, 4.5), fig_dpi=100) In line 8 we create the solver from :class:`.DiagSolver` class. Then we define the prefix of data files and the k-points on which the eigenstates will be calculated in line 9-10, where we consider the :math:`\Gamma` point. In line 11 we call the ``calc_states`` method to get the eigenvalues and eigenstates, as returned in ``bands`` and ``states``. After that, we extract the eigenstate at the Fermi level and normalize it in line 12-14. To visualize the eigenstate, we create a visualizer from the :class:`.Visualizer` class and call its ``plot_wfc`` method. Note that we must use scatter plot by setting the ``scatter`` argument to true. The eigenstate should be passed as the ``site_size`` argument, i.e., sizes of scatters will indicate the projection of eigenstate on the sites. We also need to show the model alongside the eigenstate through the ``with_model`` argument and define its plotting style with the ``model_style`` argument. Quasi-eigenstates ----------------- The quasi-eigenstate is evaluated and plotted in a similar approach: .. code-block:: python :linenos: def wfc_tbpm(sample: tb.Sample) -> None: """ Calculate wave function using TBPM. :param sample: graphene sample :return: None """ solver = tb.TBPMSolver(sample) solver.config.prefix = "tbpm" solver.config.rescale = 9.0 solver.config.num_random_samples = 1 solver.config.num_time_steps = 1024 solver.config.qe_energies = np.array([0.0]) qs = solver.calc_quasi_eigenstates() wfc = qs[0] # qs is already the squared norm wfc /= wfc.max() vis = tb.Visualizer() vis.plot_wfc(sample, wfc, scatter=True, site_size=wfc*100, site_color="b", with_model=True, model_style={"alpha": 0.25, "color": "gray"}, fig_name="tbpm.png", fig_size=(6, 4.5), fig_dpi=100) Evaluating quasi-eigenstates is a kind of TBPM calculation, so it follows the common procedure of TBPM. We firstly create the solver from :class:`.TBPMSolver` class in line 8, then set the calculation parameters in line 9-13. Especially, the ``qe_energies`` argument defines the energies for which the quasi-eigenstates will be evaluated. Then we call the ``calc_quasi_eigenstates`` method of :class:`.TBPMSolver` to get the eigenstates, normalize it, and visualize it in line 14-20. Finally, we define the driving function: .. code-block:: python :linenos: def main(): wfc_diag(make_prim_cell()) wfc_tbpm(make_sample()) if __name__ == "__main__": main() The output is shown in Fig. 1. It is clear that the eigenstate and quasi-eigenstate are localized around the vacancy with a 3-fold rotational symmetry. .. figure:: images/quasi_eigen/quasi_eigen.png :align: center :scale: 50% Spatial distribution of (a) exact eigenstate and (b) quasi-eigenstate of graphene sample with a vacancy in the center. The ``X`` marks indicate the position of the vacancy.