Phase information in EigenChannel scattering wavefunctions

[anshugaur]

Dear all,

Is there a way to visualize phase information of scattering wavefunctions using EigenChannel analysis? When I do EigenChannel analysis for a particular channel, I get "real", "imaginary", and "absolute" values of scattering wavefunctions in cube format. When I plot these using VESTA, I get "+ve" and "-ve" values of the wavefunctions (real, imaginary, or absolute), however, I do not get any phase information. Is there a way to find and visualize the phase of wavefunctions using these?

• imaginary wavefunctions are very small compared to real (4 order of magnitude) in my case.
  Thanks,
  Anshu

[tfrederiksen]

Hi @anshugaur, one possibility is to use sisl + py3Dmol in a jupyter-notebook, see https://github.com/tfrederiksen/inelastica/blob/master/tutorials/IN_02/run.ipynb. But maybe the sisl.viz module by @pfebrer can now also achieve this? It would be a great addition to Inelastica with some tools/scripts in this direction.

[pfebrer]

But maybe the sisl.viz module by @pfebrer can now also achieve this?

It does not interface with py3Dmol for now, but it is something I want to do! (because plotly 3D engine is somewhat bad for big systems and blender is just for a definitive figure). I will probably use your example here as reference :)

[tfrederiksen]

It does not interface with py3Dmol for now, but it is something I want to do! (because plotly 3D engine is somewhat bad for big systems and blender is just for a definitive figure). I will probably use your example here as reference :)

OK, but I was more thinking if your existing tools for showing isosurfaces with plotly could also handle complex wave functions (and thus could be used to also visualize eigenchannel scattering states)?

[pfebrer]

Yes, that's for sure! If you have a sisl Grid in grid, as in your notebook, then you can do

import sisl.viz

grid.plot(represent="deg_phase") # or rad_phase

But you could also store the values of the phase in a cube file, right?

import numpy as np

grid.apply(np.angle).write("Phase.cube")

PS: If you want to visualize it in 2D/1D I actually recommend using sisl visualization directly, otherwise for now I think it's a better user experience to just write it to a cube file and visualize it with some other tool.

[tfrederiksen]

Thanks @pfebrer, I was playing a bit around with this. I leave a few lines of code here that enabled me to generate some plots for a simple 1D Au chain:

import numpy as np
import sisl
import sisl.viz

# read real and imaginary part of the scattering state from cube-format
wf_re = sisl.get_sile("Chain.EC.1L.Re.cube").read_grid()
wf_im = sisl.get_sile("Chain.EC.1L.Im.cube").read_grid().copy(dtype=np.complex128)

# construction of complex-valued grid + smoothening for nicer plots
wf = (wf_re + wf_im * 1j).smooth()

# density of scattering state
dens = wf.apply(np.abs).apply(np.square)

# phase of scattering state
phase = wf.apply(np.angle)

# two plots:
dens.plot(isos=[{"frac": 0.1}, {"frac": 0.02}], axes='zx', zsmooth="best", plot_geom=True, reduce_method="sum")
phase.plot(axes='zx', colorscale="Edge", zsmooth="best", plot_geom=True, reduce_method="average")

In the latter plot we can roughly see how the traveling wave undergoes a phase shift of around pi/2 per site in the Au chain. This fits the simple picture of a half-filled s-band.

[pfebrer]

Aaah ok so inelastica produces a .Re and a .Im grid. Looks good to me!

[tfrederiksen]

It would be really cool if one could combine the above information into an isosurface contour of the density (top plot) with its surface colored according to the phase (bottom plot). Do you see a way to do this?

[pfebrer]

Yes! You need to get the isosurface vertices for the density and then find out what the phase is at those vertices.This code snippet should do it in plotly:

import sisl
import numpy as np
import plotly.graph_objects as go

# read real and imaginary part of the scattering state from cube-format
wf_re = sisl.get_sile("Chain.EC.1L.Re.cube").read_grid()
wf_im = sisl.get_sile("Chain.EC.1L.Im.cube").read_grid().copy(dtype=np.complex128)

# construction of complex-valued grid + smoothening for nicer plots
wf = (wf_re + wf_im * 1j).smooth()

# density of scattering state
dens = wf.apply(np.abs).apply(np.square)

# phase of scattering state
phase = wf.apply(np.angle)

# Get the vertices and faces of a given isosurface of the density
verts, faces, *_ = dens.isosurface(...) # Whatever arguments you want to pass

# These are the variables that plotly needs for plotting the isosurface
x, y, z = verts.T
i, j, k = faces.T

# Get the values of the phase at the vertices of the surface:
verts_indices = phase.index(verts)
color = phase.grid[tuple(verts_indices.T)]

# Build the plotly figure
fig = go.Figure(data=[
    go.Mesh3d(
        # Vertices and faces of the density isosurface
        x=x, y=y, z=z,
        i=i, j=j, k=k,
        
        # Intensity of each vertex, which will be interpolated and color-coded
        intensity = color,
        # Pass whatever colorscale that you wish
        colorscale="viridis",
    )
])

# Make it have the appropiate proportions
fig.update_layout(scene_aspectmode="data")

fig.show()

Looking forward to seeing what you produce with it! :)

[tfrederiksen]

Nice!!! This is what I get by setting dens.isosurface(0.0001) and colorscale="Edge" (to use a circular scale):

The phase discontinuity going from +pi to -pi annoys me, but I do not see an easy way to fix this.

[pfebrer]

That looks great! It looks like a worm hahaha

The phase discontinuity going from +pi to -pi annoys me, but I do not see an easy way to fix this.

Hmm that is tricky indeed.

The default is that the intensity is applied to the vertices and then the color is interpolated between them, but there is also the possibility to use intensitymode="cell". In that case you define the color of the whole face of the mesh. I guess you would then need to find the coordinates of the center of each face and get the phase for those coordinates:
https://plotly.com/python/3d-mesh/#intensity-values-defined-on-vertices-or-cells

A bit more involved but it may be worth it (?) Although for the same grid precision the plot will look less smooth.

[tfrederiksen]

With intensitymode="cell" I get an "interesting", but wrong, plot:

Maybe an issue with the indices?

[pfebrer]

Hahahah that looks beautiful indeed.

That's what I was referring to, you need to then calculate the color of each face (triangle) which is not as straightforward as calculating the color of each vertex.

Each face is defined by three indices (i,j,k as defined in the code), which contain the index of each vertex of the triangle. For each i,j,k, you would need to find the coordinates that better represent your face (I guess the center) and then you can retrieve the phase values as we did with the coordinates of the vertices.

What you are doing here is retreiving the phase at the vertices and assigning these values to faces, that's why it doesn't make sense.

[tfrederiksen]

That's what I was referring to, you need to then calculate the color of each face (triangle) which is not as straightforward as calculating the color of each vertex.

Thanks for your explanations, I went too fast here. I'll leave the calculation of face values for another time.

@anshugaur , I hope the above is useful for your purposes?

[pfebrer]

By the way, perhaps making the grid finer in the Z direction (before computing the phase) is a good patch to make the discontinuity less visible and therefore less annoying?

Something like:

new_wf_shape = list(wf.shape)
new_wf_shape[2] *= 4 # Or some other factor.

wf = wf.interp(new_wf_shape)

But I don't know how computationally expensive will things get in this particular case 😅.

If you were to provide tools like these, another approach would be to detect the discontinuity and add extra very small faces so that there is no interpolation happening between the edges of the discontinuity, but that may be too hard to do.

[tfrederiksen]

By the way, perhaps making the grid finer in the Z direction (before computing the phase) is a good patch to make the discontinuity less visible and therefore less annoying?

Excellent suggestion! Using a factor 20 along Z I now get this much more acceptable result:

[anshugaur]

Dear Thomas,
I followed the example and was able to plot both real and imaginary wave functions using sisl + py3Dmol for my system.
I need couple of clarifications:

1. En (energy), I believe it is relative to the Fermi energy (EF=0). Is this assumption correct?
2. In calculating wave functions on real-space grid for eigenchannels, "maxeig = 3 # just plot the first two channels". Does this mean the generated "real" and "imaginary" wave functions (in *.cube) files will have combined first 2 channels? Why is maxeig = 3 here? Also, when I keep maxeig > 1 (2 or 3), it only generates one set of wavefunctions (both real and imaginary, EC_0_Re/Im cube files). It shows an error "IndexError: Index is out of bounds for axis 1 with size 1)
3. Is there a way to export or save the generated wave function plots for further modification in a separate stand-alone software. We can use VESTA to plot the generated cube files, but coloring scheme for "real" and "imaginary" parts in combined plot has been a bottle neck. Any other recommended software?

Thanks a lot for your help. The exchanges between you and pfebrer have been enlightening and I will try to use some of it.

Anshu (anshugaur)

[pfebrer]

Excellent suggestion! Using a factor 20 along Z I now get this much more acceptable result:

Nice! I would probably use cubic interpolation (wf.interp(new_wf_shape, order=3)) instead of linear so that the isosurface looks smoother (if it's feasible).

[tfrederiksen]

Nice! I would probably use cubic interpolation (wf.interp(new_wf_shape, order=3)) instead of linear so that the isosurface looks smoother (if it's feasible).

Voilà!

[pfebrer]

This is some nice tutorial material :)

[tfrederiksen]

@anshugaur , based on the above solutions my recommendation would now be to first use the script EigenChannels -w cube to generate the Re/Im cube files and then to plot them with sisl.viz. The tutorial I referenced initially was meant as an intro to the eigenchannel algorithm in notebook format.

[zerothi]

The phase discontinuity going from +pi to -pi annoys me, but I do not see an easy way to fix this.

If you use a symmetric colorbar, the discontinuity would be gone, regardless of the grid.

[pfebrer]

Hmm I don't think so. Between two vertices, it seems to interpolate values, not colors. So it goes through the whole colorbar and even if it's symmetric you will see this weird feature. This is because intensity is thought for continuous scalar fields.

I've seen in the docs that there is also vertexcolor. I don't know if that does the same as intensity but interpolating colors. You could try. However that may not be what you want because you indeed want to interpolate values, it's just this particular case, where you would like to wrap around.

If this were to be a provided script to visualize phases, it would not be that difficult to detect faces where these discontinuities happens and subdivide them until it becomes practically invisible. The other option is of course to plot just the absolute value if you don't care about the sign...