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

[el-abed]

I tried both denchar and sisl and got the following shape:

Which is why i asked the question before in the previous post. I will try again soon otherwise I will send the entire folder and hope we can get to the bottom of this

[zerothi]

Hello again,
So I did specify that I want to plot the gamma using WaveFunckPoints block. But since we did 10 10 1 that means we will have more than one k point to worry about. My question is which one? It will be troublesome to plot them all using VESTA.

Like I said, just take one or two points. It doesn't really matter which ones.

Perhaps this:

WaveFuncKPointsScale ReciprocalLatticeVectors
%block WaveFuncKPoints
0.000 0.000 0.000 1 3 5
0.25 0.2 0.0 1 3 5
%endblock WaveFuncKPoints

then you could compare with the same k-points in sisl (k=[0,0,0] ; k = [0.25, 0.2, 0]).

When the grids have the same shape you can easily check by subtracting the two and finding the difference. For the skewed lattices it probably needs visual inspections. Thanks! :)

[el-abed]

then you could compare with the same k-points in sisl (k=[0,0,0] ; k = [0.25, 0.2, 0]).

When i did that, the following errors were shown:

plot = sisl.get_sile("siesta.TSHS").plot.wavefunction(k=[0,0,0] ; k = [0.25, 0.2, 0])
File "", line 1
plot = sisl.get_sile("siesta.TSHS").plot.wavefunction(k=[0,0,0] ; k = [0.25, 0.2, 0])
^
SyntaxError: invalid syntax

plot = sisl.get_sile("siesta.TSHS").plot.wavefunction(k=[0,0,0], k = [0.25, 0.2, 0])
File "", line 1
SyntaxError: keyword argument repeated
plot = sisl.get_sile("siesta.TSHS").plot.wavefunction(k=[0,0,0])
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)

[el-abed]

As for your suggestion using denchar the k points selected were empty files. Not sure what is going on

[zerothi]

Hmm...

First, I would suggest you write the grid's to a cube file and use an external tool to visualize. But @pfebrer should know what goes wrong?
I think it has to do with not reading from the fdf file. So you should change to:

plot = sisl.get_sile("RUN.fdf").plot.wavefunction(...)

or something @pfebrer may correct me here.

But again, if you follow this and store the grid to a cube file. Then you can plot it in vesta/xcrysden/vmd.

As for your suggestion using denchar the k points selected were empty files. Not sure what is going on

Could you attach your fdf files?

[pfebrer]

Yes, you need to do it from the fdf as Nick said. But it's probably better to do first the "raw" version as Nick said as well :)

[el-abed]

nick.zip
Here is the entire folder. I am really not sure why nothing is working out ?

[zerothi]

Which tool were you using for visualization?

I can see that denchar does not write correct atomic coordinates and this is what's causing some incompatibilities.

The simplest thing is to manually edit the cube files and change something like this:

$> head *.cube
 siesta.K2.WF.3.REAL.cube                                    
 siesta.K2.WF.3.REAL.cube                                    
    0    0.000000    0.000000    0.000000
   50    0.192829    0.000000    0.000000
   50    0.000000    0.192829    0.000000
   50    0.000000    0.000000    0.385659
 -0.10658E-01 -0.10681E-01 -0.10936E-01 -0.10378E-01 -0.83976E-02 -0.67459E-02
 -0.58670E-02 -0.59373E-02 -0.71887E-02 -0.10049E-01 -0.15923E-01 -0.28054E-01

into this:

$> head *.cube
 siesta.K2.WF.3.REAL.cube                                    
 siesta.K2.WF.3.REAL.cube                                    
    1    0.000000    0.000000    0.000000
   50    0.192829    0.000000    0.000000
   50    0.000000    0.192829    0.000000
   50    0.000000    0.000000    0.385659
8 0.0 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
 -0.10658E-01 -0.10681E-01 -0.10936E-01 -0.10378E-01 -0.83976E-02 -0.67459E-02
 -0.58670E-02 -0.59373E-02 -0.71887E-02 -0.10049E-01 -0.15923E-01 -0.28054E-01

note the change of the 0 to 1, and also an added line. This should make them readable with various plotters. I'll check with Siesta about this.

[el-abed]

1-Vesta is the only software i have for visualization. Any other suggestions?
2- Do you reckon it can be fixed?

[zerothi]

I think the same fix could be done for vesta, could you see whether it works by doing the above edits.

I guess it wont be fixed in denchar anytime soon. I am not an expert in that code, so it is a bit problematic for me...

[el-abed]

May I ask what the numbers 8 and 0 stand for though?
Also I changed it and still the same result

[zerothi]

it is formatted like this: Z 0. x y z (Z atomic number), (x, y, z are coordinates).

However, you are not going to use it for anything, so it doesn't really matter. ;)

[el-abed]

lol
fair enough! So what happens next? How can I help further?

[zerothi]

It would be nice to check that sisl can reproduce the same grids :)
I.e. do the same with sisl as with denchar and check that they actually match :)

[el-abed]

I did get the same issue when i used the commands for sisl. So the errors did match

[zerothi]

could you share the code you used for sisl

[el-abed]

so using python3
import sisl
import sisl.viz
plot = sisl.get_sile("siesta.TSHS").plot.wavefunction(k=(0,0,0))

Then i tried
sisl.get_sile("STRUCT5.fdf").plot.wavefunction(i=3,k=(0,0,0))

[zerothi]

Could you please try and use the tutorial I linked previously? I.e. this

Be sure to save the grid. And also note that you should have the same number of grid-points to compare numerical quantities.

[el-abed]

Hello again,
I really tried to understand the manual suggested last time and how to apply for example fat bands, PDOS, LDOS and band structure calculations using sisl to the systems, But I honestly am confused since the code seems to be applicable to only graphene.
1-I am not sure for example which Hamiltonian should i be using once I have done a SIESTA run.
2- If I have graphene/hbn or a perovskite and need to use sisl to find the band structure, is it feasible?
3- If i only want to use sisl to extract fatbands, PDOS, LDOS and so on, is it doable? Or do i have to use 1 core SIESTA job within sisl to get the Hamoltonian?

I really want to learn everything about sisl but i feel the manual needs more insight or more tutorials. Hope to hearing from you soon

[zerothi]

No, the code is not applicable to only graphene. Please try and have a go and try and understand how they function together.

1. Don't you only have 1 Hamiltonian after a succesful Siesta run?
2. yes.
3. yes it is doable in sisl if you have the Hamiltonian.

I really want to improve the manual, so any suggestion is very much welcome!

[el-abed]

1- There is the HS, HSX, TSHS files which are related to Hamiltonian.
To improve the manual, i really think a step by step procedure to be able to use sisl would be very helpful.
How to find fatbands for example? What to do exactly after a siesta run? What commands are needed to get the run right from the first time? I do not see that in the link. I see a code but not sure how to use it. I really want to see this succeed.
Hope that makes more sense.

[zerothi]

But http://zerothi.github.io/sisl/tutorials/tutorial_siesta_2.html is really a step by step procedure. I don't know what else to do? What is it about that tutorial that is not clear?

[el-abed]

So the way i understood so far is that we open python and start with the following code:
import os os.chdir('siesta_2') import numpy as np from sisl import * import matplotlib.pyplot as plt %matplotlib inline

1-For example where did siesta 2 come from? Why matplotlib is a comment? Even without comment it had an issue with inline

2- To be fair i was able to read the Hamiltonian. That was done but I was not sure was it reading from TSHS file or HSX? How many types do we have?

3-I want to focus on the band structure or PDOS example. After my siesta run is successful which had both PDOS and Band BLOCKS, I really did not understand how can I use my defined PDOS and Band structure blocks to get fat bands?

I really want to help because I even asked my colleagues if I misunderstood it. Same reaction :(

I hope you take this as a way we can make the tutorial better.

[zerothi]

The link is a jupyter notebook. I have to assume users are Python knowledgeable, I can't have everything in a tutorial. So if you are not familiar with python, and in this case jupyter notebooks, then this is the reason for your initial question.

1. The tutorial is based on the directory where I did the examples, if you know what chdir does, then this should be clear
2. If you can read the Hamiltonian, then you are ready. Either file you read is fine, they both contain the same information.
3. Please carefully go through the tutorial I linked and try and understand how you can change the lines to get different fat-bands. It is not clear to me whether you have issues understanding what a fatband plot is or what sisl does/can do. If you have the PDOS for all states + the bandstructure, then you can easily plot the fatbands...

Also, note that this original issue is not related at all with PDOS and fatbands. Denchar does an entirely different thing.

Regarding your colleagues, what question did you ask them, if it was sisl specific and they have never used the tool, then how could they know ;)
I really want to improve the tutorial, but I still don't know where you go wrong. From your questions it seems you at least lack some basic python knowledge.

[el-abed]

1- I think I understand the main issue and that is on me. If I collect all the bits and pieces in that page and put them into a python file like Papior.py Then use python Papior.py I will be able to get PDOS and band structure is that correct?
2- My other comment is if I already have my Hamiltonian within certain K pts and energy range, do I still have to define K and E?

[zerothi]

1. Yes, you could do that.
2. TSHS and HSX contains the Hamiltonian at no specific k-point. You can generate the Hamlitonian for any k using that information.

[el-abed]

Thank you a lot for your help and sorry for taking away the DENCHAR post.
But:
2- I do not understand. How can we extract the k pts and energy range we chose in our fdf file without defining them in the python code? Or should i just reintroduce E and K?

[zerothi]

2. I think you are misunderstanding what the Hamiltonian contains and what the E and k means? Please consider what each of them mean in the context of analysing DOS, how are k and E related to the DOS and examine the importance of the sampling of these as well.

[el-abed]

I guess my question is for large systems where we already introduced DOS and band blocks, why should we redefine them again in the python code? That is the only thing i do not understand.

[zerothi]

I guess my question is for large systems where we already introduced DOS and band blocks, why should we redefine them again in the python code? That is the only thing i do not understand.

But if you already calculated it with siesta, then there is no point in doing it with sisl, it is just calculating the same thing again?
However, in some cases it may be a lot easier to do convergence tests in sisl than in siesta since siesta requires a couple of SCF steps before analysis.