Scissors operator for LCAO mode

Warning

Work in progress

gpaw.lcao.scissors.non_self_consistent_scissors_shift(shifts, dft)

Apply non self-consistent scissors shift.

Return eigenvalues ase a:

(nspins, nibzkpts, nbands)

shaped ndarray in eV units.

The shifts are given as a sequence of tuples (energy shifts in eV):

[(<shift for occupied states>,
  <shift for unoccupied states>,
  <number of atoms>),
 ...]

Here we open a gap for states on atoms with indices 3, 4 and 5:

eig_skM = non_self_consistent_scissors_shift(
    [(0.0, 0.0, 3),
     (-0.5, 0.5, 3)],
    dft)

In LCAO Mode we solve the following generalized eigenvalue problem:

\[\sum_\nu (H + \Delta H)_{\mu\nu} C_{\nu n} = \sum_{\nu} S_{\mu\nu} C_{\nu n} \epsilon_n,\]

where \(\Delta H\) is a scissors operator.

Space is divided into regions \(\Omega_i\) and for each region we define desired shifts of the occupied and unoccupied bands: \(\Delta^i_{\text{occ}}\) and \(\Delta^i_{\text{unocc}}\).

For each region, we diagonalize the density-matrix

\[\rho_{\mu\nu} = \sum_n C_{\mu n} f_n C_{\nu n}^*\]

in the orbitals belonging to \(\Omega_i\):

\[\sum_{\nu\in i}(S^{1/2}\rho S^{1/2})_{\mu\nu} V^i_{\nu\alpha} = \lambda^i_\alpha V^i_{\mu\alpha}.\]

Here, the eigenvalues \(\lambda^i_\alpha\) will be close to either zero or one and the scissors operator is now given as:

\[\Delta H_{\mu\nu} = \sum_i \sum_{\mu'\in i,\nu'\in i} \sum_\alpha S^{1/2}_{\mu\mu'} V^i_{\mu'\alpha} (\lambda^i_\alpha \Delta^i_{\text{occ}} + (1 - \lambda^i_\alpha) \Delta^i_{\text{unocc}}) V^{i*}_{\nu'\alpha} S^{1/2*}_{\nu'\nu}.\]

WS2 layer on top of MoS2 layer

Band structures for:

• no shifts (shifts=[])

• MoS2 gap opened up by 1.0 eV (shifts=[(-0.5, 0.5, 3)])

• MoS2 shifted up by 0.5 eV and WS2 down by 0.5 eV (shifts=[(0.5, 0.5, 3), (-0.5, -0.5, 3)])

from ase.build import mx2
from gpaw.new.ase_interface import GPAW


def mos2wds(shifts: list[tuple[float, float, int]], tag: str) -> None:
    """WS2 layer on top of MoS2 layer."""
    atoms = mx2(formula='MoS2', kind='2H', a=3.184, thickness=3.13,
                size=(1, 1, 1))
    atoms += mx2(formula='WS2', kind='2H', a=3.184, thickness=3.15,
                 size=(1, 1, 1))
    atoms.positions[3:, 2] += 3.6 + (3.13 + 3.15) / 2
    atoms.positions[3:] += [-1 / 3, 1 / 3, 0] @ atoms.cell
    atoms.center(vacuum=6.0, axis=2)
    k = 6
    atoms.calc = GPAW(mode='lcao',
                      basis='dzp',
                      nbands='nao',
                      kpts=dict(size=(k, k, 1), gamma=True),
                      eigensolver={'name': 'scissors',
                                   'shifts': shifts},
                      txt=f'{tag}.txt')
    atoms.get_potential_energy()
    bp = atoms.cell.bandpath('GMKG', npoints=50)
    bs_calc = atoms.calc.fixed_density(kpts=bp, symmetry='off')
    bs_calc.write(f'{tag}.gpw')
    bs = bs_calc.band_structure()
    bs.write(f'{tag}.json')


if __name__ == '__main__':
    for i, shifts in enumerate([[],
                                [(-0.5, 0.5, 3)],
                                [(0.5, 0.5, 3), (-0.5, -0.5, 3)]]):
        mos2wds(shifts, f'mos2ws2-{i}')

import matplotlib.pyplot as plt
import numpy as np
from ase.units import Ha
from gpaw.new.ase_interface import GPAW
from matplotlib.collections import LineCollection


def line_segments(x_k: np.ndarray, y_nk: np.ndarray) -> np.ndarray:
    """Helper function for plotting colored bands.

    Converts (x,y) points to line segments.
    """
    N, K = y_nk.shape
    S_nksv = np.empty((N, K, 3, 2))
    S_nksv[:, 0, 0, 0] = 0.0
    S_nksv[:, 1:, 0, 0] = 0.5 * (x_k[:-1] + x_k[1:])
    S_nksv[:, :, 1, 0] = x_k
    S_nksv[:, :-1, 2, 0] = 0.5 * (x_k[:-1] + x_k[1:])
    S_nksv[:, -1, 2, 0] = x_k[-1]
    S_nksv[:, 0, 0, 1] = np.nan
    S_nksv[:, 1:, 0, 1] = 0.5 * (y_nk[:, :-1] + y_nk[:, 1:])
    S_nksv[:, :, 1, 1] = y_nk
    S_nksv[:, :-1, 2, 1] = 0.5 * (y_nk[:, :-1] + y_nk[:, 1:])
    S_nksv[:, -1, 2, 0] = np.nan
    return S_nksv.reshape((N * K, 3, 2))


def plot(ibzwfs, bp, ax):
    x_k, xlabel_K, label_K = bp.get_linear_kpoint_axis()
    label_K = [label.replace('G', r'$\Gamma$') for label in label_K]

    eig_kn = []
    color_kn = []
    for wfs in ibzwfs:
        c_an = [(abs(P_ni)**2).sum(1) for P_ni in wfs.P_ani.values()]
        c_n = sum(c_an[:3]) / sum(c_an)
        color_kn.append(c_n)
        eig_kn.append(wfs.eig_n * Ha)

    eigs = line_segments(x_k, np.array(eig_kn).T)
    colors = np.array(color_kn).T.copy().flatten()
    lc = LineCollection(eigs)
    lc.set_array(colors)
    lines = ax.add_collection(lc)
    ax.set_xlim(0, x_k[-1])
    ax.set_ylim(-10, 0)
    ax.set_xticks(xlabel_K)
    ax.set_xticklabels(label_K)
    return lines


if __name__ == '__main__':
    fig, axs = plt.subplots(1, 3, sharey=True)
    for i, ax in enumerate(axs):
        bs_calc = GPAW(f'mos2ws2-{i}.gpw')
        bp = bs_calc.atoms.cell.bandpath('GMKG', npoints=50)
        lines = plot(bs_calc.dft.ibzwfs, bp, ax)
    cbar = fig.colorbar(lines)
    cbar.set_ticks(ticks=[0, 1], labels=['W', 'Mo'])
    # plt.show()
    plt.savefig('mos2ws2.png')

Tip

You can plot the JSON band-structure files with the command: ase band-strucuture <name>.json.