How to visualize MPSH

Asked by zhikai on 2018-12-02

Dear siesta developers and users

I am trying to visualize MPSH for some energy values after calculating MPSH eigenvalues via sisl https://github.com/zerothi/sisl.

Is there any way to visualize these orbitals such as HOMO and LUMO.

Any help will be appreciated!

best regards
Zhikai

Question information

Language:
English Edit question
Status:
Solved
For:
Siesta Edit question
Assignee:
No assignee Edit question
Solved by:
Nick Papior
Solved:
2018-12-11
Last query:
2018-12-11
Last reply:
2018-12-05
Nick Papior (nickpapior) said : #1

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.)

zhikai (zhikaizhao) said : #2

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

Best Nick Papior (nickpapior) said : #3

No, what I wrote was correct.

Please try and print H using your method and using my method, then you should see a difference.

zhikai (zhikaizhao) said : #4

Thanks Nick Papior, that solved my question.