

In [ ]:

    
from __future__ import print_function
import sisl
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In this analysis example we will show how to plot the wavefunction for a periodic system (the same scheme may be used to plot molecular orbitals).

The basic principle of plotting the real-space wave functions may be written as this (for a given $\mathbf k$-point): \begin{equation} \psi_i(\mathbf k, \mathbf r) = \sum_\nu e^{i\mathbf k \cdot \mathbf r_\nu}c_{i\nu}\phi_\nu(\mathbf r - \mathbf r_\nu), \end{equation} where $c_{i\nu}=\langle\tilde\phi_\nu|\psi_i\rangle$ is the coefficient for the $i$th eigenstate and $\nu$ basis orbital and $\phi_\nu(\mathbf r - \mathbf r_\nu)$ is the basis orbital $\nu$ centered at $\mathbf r_\nu$.

sisl will in most cases, automatically read the necessary information from a Siesta run to be able to construct the $\phi_\nu$ basis functions in real-space. If the basis-information is not available sisl will inform you when trying to calculate real-space quantities.

In this example we will calculate the eigenstates for a given $\mathbf k$-point for graphene and plot a cut through the real-space wavefunction at a fixed distance above the graphene plane.

Exercises

Run Siesta.
Read in the Hamiltonian using the RUN.fdf file, see e.g. S 1.
Calculate the eigenstate for the $\Gamma$-point (see. Hamiltonian.eigenstate)
Read the entry about the EigenstateElectron in the sisl.
Figure out which method you should use in order to calculate the real-space wavefunction.
Use the below method (plot_grid) to plot a cut through the wavefunction grid. HINT: this may be a useful grid grid = Grid(0.05, sc=H.geometry.sc)
If you don't see a warning like this:
```
 info:0: SislInfo: wavefunction: summing 18 different state coefficients, will continue silently!
```
then you have done it correctly! Look in the documentation and figure out how to take a subset of the eigenstates and only plot a single one of them on the grid. This is important since otherwise you are plotting the super-position of all eigenstates at the $\mathbf k$-point.
HINT: If you have VMD/XCrySDen on your laptop it may be fun to plot the real-space quantities using cube files, if you want, figure out how to save the Grid into a cube/xsf file.
Try and play with the supercell you pass to the Grid initialization and plot the wavefunction for different sizes of the grid.
What do you see for increasing size of grids? Are there any periodicites, if so, what is the periodicity and why?
Calculate the eigenstates such that the wavefunctions has a periodicity of $3$ along the first lattice vector and $2$ along the second lattice vector (only plot one of the eigenstates)
HINT: $e^{i\mathbf k \cdot \mathbf R}$



In [ ]:

    
def plot_grid(grid, plane_dist=1):
    """ Plot the grid in either 1 (real-only) or 2 (complex wavefunction) graphs
    
    A cut through the grid will be plotted corresponding to `plane_dist` above the z-coordinate
    """
    z_index = grid.index(plane_dist, 2)
    x, y = np.mgrid[:grid.shape[0], :grid.shape[1]]
    dcell = grid.dcell
    x, y = x * dcell[0, 0] + y * dcell[1, 0], x * dcell[0, 1] + y * dcell[1, 1]
    if grid.dtype in [np.complex64, np.complex128]:
        fig, axs = plt.subplots(1, 2, figsize=(15, 5))
        axs[0].contourf(x, y, grid.grid[:, :, z_index].real)
        im = axs[1].contourf(x, y, grid.grid[:, :, z_index].imag)
        for ax in axs:
            ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]')
        axs[0].set_title('Real part')
        axs[1].set_title('Imaginary part')
    else:
        fig, ax = plt.subplots(1, 1) ; axs = [ax]
        axs = [ax]
        im = ax.contourf(x, y, grid.grid[:, :, z_index])
        ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]')
    # Also plot the atomic coordinates
    try:
        xyz = grid.geometry.xyz
        for ax in axs:
            ax.scatter(xyz[:, 0], xyz[:, 1], 50, 'k', alpha=.6)
    except: pass
    fig.colorbar(im);



In [ ]:

    
# Add code here to 1) read in Hamiltonian with basis orbitals, 2) calculate the eigenstate for a given k-point,
# 3) create a grid to plot the wavefunction on and 4) plot the grid using the above method



In [ ]: