How to visualize MPSH

This snippet should show the required information with some details, please search the documentation of sisl for further details.

from __future__ import print_function
import sisl as si
import numpy as np

# Read in geometry
C60 = si.get_sile('RUN.fdf').read_geometry()
# Retrieve all C60 atoms:
a_C60s = (C60.atoms.Z == 6).nonzero()[0]

# Read in Hamiltonian (with basis-information)
H = si.get_sile('RUN.fdf').read_hamiltonian()
# Reduce the Hamiltonian to a sub-region
# This is also known as Molecular Projected Self-consistent Hamiltonian (MPSH)
H = H.sub(a_C60s)

# Find bounding box for the C60's
xm = H.geom.xyz.min(0) - H.geom.maxR()
xM = H.geom.xyz.max(0) + H.geom.maxR()
print(xm, xM)
# Define the region of space to plot the eigenstates in
sc = si.SuperCell(xM - xm, origo=xm)

# Create a list of eigenstates to plot
plot_states = range(236, 248)

# First do the molecular orbital from a periodic picture
es = H.eigenstate()

# Specify the grid to plot the molecular orbitals in
# The first number defines the "voxel" side-lengths
# and is thus a single number that uniquely defines the same
# voxel size along each direction
g = si.Grid(0.1, sc=sc)
print(g)

print(es.eig[plot_states])
for state in plot_states:
    print('Creating: {}'.format(state))
    es.sub(state).wavefunction(g)
    print('Writing: {}'.format(state))
    si.get_sile('p_mo_{}.cube'.format(state), 'w').write_grid(g)
    # Reset grid such that we don't create a superposition
    g.fill(0.)

# Try and plot it again without periodicity
H.set_nsc([1, 1, 1])
es = H.eigenstate()

# Specify the grid to plot the molecular orbitals in
# In this case we know that we should only plot the non-periodic
# images, and thus we need to force the boundaries to be anything
# but periodic. By default a Grid is sisl IS periodic.
g = si.Grid(0.1, sc=sc, bc=si.Grid.OPEN)
print(g)

# Create a list of eigenstates to plot
print(es.eig[plot_states])

for state in plot_states:
    print('Creating: {}'.format(state))
    es.sub(state).wavefunction(g)
    print('Writing: {}'.format(state))
    si.get_sile('mo_{}.cube'.format(state), 'w').write_grid(g)
    # Reset grid such that we don't create a superposition
    g.fill(0.)

Dear Papior

Thank you for your reply! It helps a lot.

In your snippet,
# Read in Hamiltonian (with basis-information)
H = si.get_sile('RUN.fdf').read_hamiltonian()

do you mean " H = si.get_sile('xxx.TSHS').read_hamiltonian()" ?
the file 'xxx.TSHS' is generated after Transiesta calc under bias voltage

zhikai

Thanks Nick Papior, that solved my question.