Inquiry about Denchar equivalence in sisl

[el-abed]

**Good afternoon, I recently asked about utilizing denchar for my calculations. Nick mentioned that it could be possible to do denchar calculations using sisl and I would like to test it to calculate the HOMO and LUMO of my systems. May I ask how can I properly use sgrid to find them? **

Version details For example in denchar we can start by using denchar -k 3 -w 4 file.fdf
which will plot only the wave-function with (original) index 4 of the third k-point in the (WFSX) file. How can we do the same with sisl? I assume it would be in sgrid in this case. How? Which nc file should I find? Thank you and looking forward to your thoughts EL-abed

[pfebrer]

Hey! While you wait for Nick's answer, which will have probably other ways to do it, I just want to let you know that a first version of the visualization module that I was developing has been already merged to sisl. For now, you can install it from github, as there has not been a release since then:

pip install "sisl[viz] @ git+https://github.com/zerothi/sisl.git"

(where [viz] is to add the extra dependencies that only the visualization module uses). You can also of course just git pull if you are already using it from a git clone, but then you will have to install the extra dependencies manually.

Then, you can use it to plot wavefunctions like this:

import sisl
import sisl.viz

sisl.get_sile("file.fdf").plot.wavefunction(i=4, k=(0,0,0))

For this, you will need to have the hamiltonian (in any of its forms; HSX, TSHS, nc file) and the basis files of each atom (i.e. .ion or .ion.xml files).

If you want to write the grid to a file, you can just store the plot in a variable, and find the grid under its grid attribute.

plot = sisl.get_sile("file.fdf").plot.wavefunction(i=4, k=(0,0,0))
# The grid where the wf is stored can be found at plot.grid.
wfgrid = plot.grid
# Write it
wfgrid.write("wf4.cube")

However, note that the visualization module is already very flexible on visualizing grids, so you can try to use the settings of the plot that can be tweaked. For example, you can quickly switch the visualization from 3d to 2d or 1d, toggle the geometry, apply some transforms....

# Check all the available settings
print(plot.settings)
# You can update the settings using update_settings
plot.update_settings(axes=[0,1], plot_geom=False, transforms=["square"])

You could also plot from the terminal:

splotly wavefunction -f file.fdf --i 4

But there is no direct way to store the grid to a file using this by now.

Documentation for all the visualization module is already written and I hope it can be quite illustrative of all the options. For some reason it was failing when Nick tried to build it, but I was just going to look into it when I saw your issue :)

Hope it helps!

[el-abed]

Thank you very much @pfebrer for the help. I am eager to test the new software. I successfully installed viz where i can see it clearly in my home directory. The problem is that unlike previously installing sgrid or sisl, i could not find the executable file splotly. Is there a reason for that ? Thank you and looking forward to your thoughts El-Abed

[pfebrer]

Great! Thanks for checking it out.

I see that Nick commented out the CLI when he merged. That is because we are still not sure how the CLI should look like for plotting stuff. So basically, for now you would have to do it using python as I mentioned here:

import sisl
import sisl.viz

sisl.get_sile("file.fdf").plot.wavefunction(i=4, k=(0,0,0))

Otherwise you could go to the place where sisl is installed and modify the file sisl/viz/plotly/splot by adding:

if __name__ == "__main__":
    splot()

then you could use

python -m sisl.viz.plotly.splot

as a replacement for splotly. But I think that's dirty and I encourage you to use it from python until we find a satisfactory CLI :) It is much more flexible if you run a jupyter notebook, because once you have the plot you can tweak it very easily.

[el-abed]

Thank you very much for the quick reply. I followed the more suitable advice which is working right now. I assume first that i=4 means the wave function number 4 right? My main feedback is asked as follows: in sgrid there is a nice command in order to use the relaxed structure found in the XV file such as the following: sgrid PG.DeltaRho.grid.nc --geometry zeroPG.XV --diff P.DeltaRho.grid.nc --diff G.DeltaRho.grid.nc diffn.cube

The reason i am asking is because if i do not do that I will get a figure such as: where the isosurfaces are far away. Any advice on such? Thank you again @pfebrer and hope to hearing from you soon! EL-abed

[pfebrer]

For now you'll have to set the geometry before writing the file. Hadn't thought about this!

wfgrid = plot.grid
wfgrid.geometry = sisl.Geometry.new("zeroPG.XV")
wfgrid.write("wf.cube")

Sorry there isn't a faster way for now :)

You can check what each setting means by getting the help of the plot that you just got. I.e.:

help(plot)

[pfebrer]

Sorry, there's a faster way: if you use the hamiltonian instead of the fdf file, it will use the last geometry of the siesta run I believe. That is:

import sisl
import sisl.viz

plot = sisl.get_sile("file.TSHS").plot.wavefunction(i=4) # I'm using the TSHS, but it could be HSX or the nc hamiltonian
plot.grid.write("wf4.cube")

If you want to set any other geometry then you'll have to do it as in the previous comment.

Let me know if it works for you!

[el-abed]

@pfebrer thank you very much for the reply. 1- Concerning your suggestion here wfgrid = plot.grid wfgrid.geometry = sisl.Geometry.new("zeroPG.XV") wfgrid.write("wf.cube")

The code did not show any error. However, the final result remained the same. Any reason why?

2- Concerning the last suggestion i got the following error when i used HSX file: ValueError: hsxSileSiesta xij(orb) -> xyz did not find same coordinates for different connections It does not make sense since it belongs to the same run.

Thank you and looking forward to your reply. EL-abed

[pfebrer]

Hmm, and when you do:

sisl.Geometry.new("file.XV").plot()

are the atoms inside the unit cell?

[pfebrer]

Can you send me the files (fdf, H, .ion and XV) so that I can check?

[el-abed]

Here are the files(fdf, H, .ion and XV)

sisl.zip

[pfebrer]

Doing sisl.get_sile("file.fdf").plot.wavefunction(i=4, k=(0,0,0)) I get the error that you were getting here:

ValueError: hsxSileSiesta xij(orb) -> xyz did not find same coordinates for different connections

Do you have the hamiltonian in some other form (.TSHS or .nc)? Because probably it was using a different one when it worked.

[el-abed]

Yes you will get that error because you need to use my fdf :) So you need to say sisl.get_sile(""0.158_P_G.fdf"").plot.wavefunction(i=4, k=(0,0,0))

[pfebrer]

Yes, of course I used your fdf, sorry :)

Hmm I don't know why I get this, maybe you can send the entire folder

[zerothi]

The hsx file implementation is not really good and may be buggy.

I would heavily suggest you to use tshs or nc file format as @pfebrer suggests. But we don't know if you are using 4.1?

[zerothi]

Also indices in python are 0 based, so i=2 means the 3rd wavefunction.

[el-abed]

1- @pfebrer I should be more clear, when i used sisl.get_sile(""0.158_P_G.fdf"").plot.wavefunction(i=4, k=(0,0,0)) I got no error whatsoever. 2- @zerothi Thank you for the reply. I am using Siesta Version : v4.1-b4. Was wondering which command allows me to get the Hamiltonian nc files? Because Im doing a siesta not transiesta calculation so tshs is not part of my calculation. 3- Oh i forgot about python being 0 based. Thanks for that. 4- I want to send the entire folder but even when zipped its massive (>10 GB) due to the WFSX file. Any other option?

[pfebrer]

sisl.get_sile(""0.158_P_G.fdf"").plot.wavefunction(i=4, k=(0,0,0))

That's what I did, but for some reason I get the error.

Was wondering which command allows me to get the Hamiltonian nc files? Because Im doing a siesta not transiesta calculation so tshs is not part of my calculation.

TS.HS.Save true

You can use it in regular siesta calculations.

I want to send the entire folder but even when zipped its massive (>10 GB) due to the WFSX file. Any other option?

No need to send it then :) You can try with the TSHS hamiltonian, should work.

[el-abed]

@zerothi @pfebrer Hello again! When i used the following command: plot = sisl.get_sile("158PG.TSHS").plot.wavefunction(i=4, k=(0,0,0)) I got the following error: The plot has been initialized correctly, but the current settings were not enough to generate the figure. I want to figure it out on my own but Honestly I do not know what went wrong? The job ran successfully where I only added TS.HS.Save true Any thoughts?

[pfebrer]

The plot has been initialized correctly, but the current settings were not enough to generate the figure.

Could you send the full error message? I think I know what it is though. I forgot the .TSHS doesn't has the basis stored. So you actually have to pass a geometry with a basis:

geometry = sisl.get_sile("file.fdf").read_geometry(output=True) # This will read it from your XV
plot = sisl.get_sile("158PG.TSHS").plot.wavefunction(i=4, k=(0,0,0), geometry=geometry)

but I think it should be equivalent to the

sisl.get_sile("file.fdf").plot.wavefunction(i=4, k=(0,0,0))

approach. Because it is actually using the XV file. At least this is how I meant it. It makes no sense to plot the wavefunction using the input structure :sweat_smile:. So, if you see different results, then I will need to check what's wrong with the code.

If not, are you sure the XV coordinates are inside the unit cell? When you do this:

sisl.Geometry.new("file.XV").plot()

do atoms show up inside the drawn unit cell? If not, what happens if you:

plot = sisl.get_sile("file.fdf").plot.wavefunction(i=4, k=(0,0,0))
plot.geometry.plot().show()
# and
plot.H.geometry.plot().show()

By the way, you are the first person that it's not me trying the WavefunctionPlot, so thanks and sorry :)

[el-abed]

The entire error is: The plot has been initialized correctly, but the current settings were not enough to generate the figure. (Error: wavefunction: Cannot create wavefunction since no atoms have an associated basis-orbital on a real-space grid)

Which i got the same when I used: geometry = sisl.get_sile("file.fdf").read_geometry(output=True) # This will read it from your XV plot = sisl.get_sile("158PG.TSHS").plot.wavefunction(i=4, k=(0,0,0), geometry=geometry) The geometry line worked fine. The second line however gave me the same error. Here is my XV file which i had to convert to txt file. I am not sure if the zero's in the end are the reason of the mistake?

158PGtxtfile.txt Because I am not sure what is the out put of the command: sisl.Geometry.new("file.XV").plot()? As for WavefunctionPlot and any code, I am happy to help man! Really feel proud to help out! So what ever code you have happy to help! Looking forward to your reply!

[pfebrer]

Aah ok depending on where you are running you need to show the plot for it to display:

sisl.Geometry.new("file.XV").plot().show()

[pfebrer]

You can send me the .TSHS and I can try again :)

[el-abed]

158PGtshsfile.zip Here you go the TSHS file in zip format. Let me know how it goes

[zerothi]

@el-abed have you seen these tutorials: http://zerothi.github.io/sisl/tutorials/tutorial_siesta_1.html http://zerothi.github.io/sisl/tutorials/tutorial_siesta_2.html

They both create wavefunction outputs, but also real-space charges. This is what @pfebrer's backend does behind the scenes :)

[el-abed]

Oh! So I will have to start using Jupiter to be able to see what is going on more clearly. @zerothi Thank you! will have a look when I can. But I just want to know if my job ran wrong or is there something missing in the code?

[pfebrer]

@zerothi I think the issue @el-abed is having might be because of some bug with the writing of the grid to a "cube" file. If I do:

plot = sisl.get_sile("0.158_P_G.fdf").plot.wavefunction(i=50)
plot.show()

I get:

which I don't know if makes sense -maybe there is some bug in plotting the wavefunction-. But the geometry is inside the unit cell (I don't do any modification to the coordinates, you can check at plot.grid.geometry).

Now if I write the grid:

plot.grid.write("wf.cube")

it shows up in VESTA like this:

[pfebrer]

Seeing the wf.cube file I don't know what could be happening. Seems to be fine. Maybe a problem with the second lattice vector having a negative component?

[zerothi]

Yeah, the wavefunction and density methods really needs some comparisons! It would be of great value is somebody took the time (do you sign up for this @el-abed, that would be great!!!!) to do some system calculations and compare with denchar

I would recommend 3 systems, and ideally both in a non-orthogonal and orthogonal configuration (so 6 systems, really ;))

1. hBN (because it has two species)
2. AB-stacked hBN
3. AB-stacked hBN with atomic coordinates outside of the unit-cell

Then start with trying to recreate the density from denchar and using sisl.

import sisl
DM = sisl.get_sile("RUN.fdf").read_density_matrix()
grid = sisl.Grid(0.02, geometry=DM.geometry)
DM.density(grid)
grid.write("same_uc.cube")
# now compare with denchar orthogonal cell
grid = sisl.Grid(0.02, sc=[10, 10, 100])
DM.density(grid)
grid.write("sq_uc.cube")

this should compare how sisl works for skewed and non-skewed unit-cells. :)

Once this is cleared, then it will be easier to do the wavefunction (they are more or less the same with some details regarding phases etc.

Seeing the wf.cube file I don't know what could be happening. Seems to be fine. Maybe a problem with the second lattice vector having a negative component?

Yeah, that would be the cause. Some viewers don't like this, but I really don't know vesta, could you try with vmd or xcrysden?

[pfebrer]

Yeah, that would be the cause.

Yep, tried to turn the negative value to positive and now the third axis is fine (ofc the modified one is not):

I tried in vmd but I don't know how to see the volume data (it shows me nothing :sweat_smile:). I read the cube specifications and it says this:

The second through fourth values on this line are the components of the vector X⃗ defining the voxel X-axis. They SHOULD all be non-negative; proper loading/interpretation/calculation behavior is not guaranteed if negative values are supplied.

[zerothi]

Yeah ok. So it is a general limitation.

I am not surprised your geometry is now outside the unit-cell ;)

You'll have to pass a grid with a supercell that is corrected in order to get the correct lattice vectors, perhaps just doing:

sc = H.geometry.sc.copy()
sc.cell[1, :] *= -1
grid = Grid(0.02, sc=sc)
eigenstate.wavefunction(grid)

or something similar, that should do things correctly. But still, the lowest plane seems off?