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?

[pfebrer]

Yes, but the correction should be done when writing to the cube file, no? I don't know if there's a general way of setting all lattice vectors to positive. Would shifting the geometry make sense? Like:

Original cell: [ [1,0,0], [0,-1,0], [0,0,1] ]

Final cell: [ [1,0,0], [0,1,0], [0,0,1] ]

And, given that coordinates are not fractional:

geometry.move([0,1,0])

Because mirroring (keeping coordinates fractional and just changing the sign) would change the geometry, right?

[pfebrer]

But still, the lowest plane seems off?

Yes, definitely

[pfebrer]

Would shifting the geometry make sense?

Ah no, but the grid also needs to be modified. In this case it's just a flip, but is it possible to generalize? Is this why you proposed to change the geometry before calculating wavefunctions? :smile:

[zerothi]

Hmm. Perhaps doing this would be ok

grid = ...
eigenstate.wavefunction(grid) ...

# Now correct origo
sc = grid.sc
sc.origo = sc.origo + sc.cell[1, :]
sc.cell[1, :] *= -1

grid.write("test.cube")

Hmm.. I have to think about this, if it actually just mirrors the data, or if one can do this.

[el-abed]

1. @zerothi I do not mind testing the software to be honest.
2. I was wondering how can I get the latest denchar software? Because it was not a long time ago when I downloaded the latest siesta version from github.
3. For the 6 systems you mentioned how many atoms do you want me to build? How big? How far away the atoms should be far away from the unit cell?

[zerothi]

@el-abed many thanks!!

Sorry for the late, return. I haven't forgotten about it... :)

So I created a small script that should create some structures that would be interesting to check.

import sisl as si

write = False

# Create hBN
B = si.Atom('B')
N = si.Atom('N')

itest = 0

dist = 1.44
for orthogonal in (False, True):
    # simple hBN
    hBN = si.geom.graphene(dist, [B, N], orthogonal=orthogonal)
    # Add some vacuum to be sure
    hBN = hBN.append(si.SuperCell(10), axis=2)

    if write:
        hBN.write(f"STRUCT{itest}.fdf")
    itest += 1

    # Then also created AB stacked
    hBN_B = hBN.translate([dist, 0., 3.35])

    hBN_AB = hBN.add(hBN_B)
    hBN_AB = hBN.add(hBN_B)
    if write:
        hBN_AB.write(f"STRUCT{itest}.fdf")
    itest += 1

    # now move top layer
    dxy = hBN_AB.cell[[0, 1]].sum(0)
    hBN_AB = hBN_AB.translate(dxy, atoms=range(hBN_AB.na // 2, hBN_AB.na))

    if write:
        hBN_AB.write(f"STRUCT{itest}.fdf")
    itest += 1

So there are quite a few of these structures. Probably the easiest to do would be to first test one of them. Then get a feeling for the procedure.

Also, it may be easier to start with the orthogonal=True cases since that is what Denchar intrinsically produces (it doesn't produce skewed lattices.

Please use this siesta/denchar version.

The accuracy of the calculations are not important since here we are extracting the same information from the same data files. Just be sure to add TS.HS.Save to the fdf file to be sure to get the information. :)

If you have any questions, just pop them in here :) THANKS!

[el-abed]

Thank you very much for the message @zerothi ! Really eager to see how that goes! I do have a couple of questions: 1- Concerning denchar, I assume i should reinstall siesta and then install all utilities is that correct? 2- Concerning k pts, cut off energy and such do you need me to do some convergence tests just to be as realistic as possible? I know you said the final answer doesnt have to be accurate but at least we could aim for good numbers right? 3- I love how the code gave me the 6 fdf files we need instead of having to put true then false to get 3 by 3. LOVE IT :) Thank you very much again Nick! Happy to be part of this mini project! Eager to see how it goes.

[zerothi]

Thank you very much for the message @zerothi ! Really eager to see how that goes! I do have a couple of questions: 1- Concerning denchar, I assume i should reinstall siesta and then install all utilities is that correct?

Yeah, use the link provided. And all utilities you use should be from that tar.gz :)

2- Concerning k pts, cut off energy and such do you need me to do some convergence tests just to be as realistic as possible? I know you said the final answer doesnt have to be accurate but at least we could aim for good numbers right?

Maybe use 300 Ry, and 15x15x1 kpoints. This should be ok. Use Diag.ParallelOverK true in fdf to speed up things. If it is too demanding, just go lower in k-points (but retain at least a couple).

3- I love how the code gave me the 6 fdf files we need instead of having to put true then false to get 3 by 3. LOVE IT :) Thank you very much again Nick! Happy to be part of this mini project! Eager to see how it goes.

Great!! :) I am too :)

[el-abed]

Hi @zerothi The latest siesta installation has been a problem so I raised it to the supercomputer administration to fix the issue. In the mean time, would you mind if i do the calculation using the current siesta i have 4.1b4? Then once the issue is resolved I can use denchar in the latest version?

[zerothi]

I would prefer if that version is used, there have been some fixes. What was the problem?

[el-abed]

No problem! With such supercomputer its usually a mixture of finding the right modules and netcdf issues. So i have been told to wait for someone from administration to assist me unfortunately. :(

[pfebrer]

Hey @el-abed you can try to install siesta using spack (thanks to @vdikan). It knows what dependencies SIESTA has and installs them locally (you don't need anything from the supercomputer apart from the C and fortran compilers and python).

However it sometimes has trouble to work in clusters due to restrictions. You can give it a try and if it doessn't work you could also install it in your workstation, since the systems that Nick proposed are very small.

The rel-4.1 branch (which Nick was proposing) is not an option as of now, but I already asked Vladimir to include it. Meanwhile you can try if it works by installing siesta@master. If it does, it is very straightforward to install the one that you need.

Here's a guide that I assembled for a student that I'm supervising:

Installing SIESTA in parallel using spack (insultingly easy)

Make sure you have git:

sudo apt-get update
sudo apt-get install git

Clone the spack repo from Vladimir:

git clone -b siesta-develop https://github.com/vdikan/spack.git

Then make sure the spack environment is loaded when you open the terminal (just so that you can run the spack command).

echo "source /path/to/spack/share/spack/setup-env.sh" >> ~/.bashrc
echo "source /path/to/spack/share/spack/spack-completion.bash" >> ~/.bashrc

where you should replace /path/to/spack with the path to the spack folder in your computer.

Close and open the terminal. You should be able to run the spack command now.

Do

spack info siesta

to check that SIESTA is there and see all the options that you have. Then just:

spack install siesta

When you looked at the info of the siesta package you could see that there is a psml version available. To specify the version to install just use "@". E.g.:

spack install siesta@psml

When it finishes (it will take time, just go do something else in the meantime), you should be able to load siesta by doing:

spack load siesta

Now you have the siesta command, as well as mpirun. So you can try to run it in 2 processors:

mpirun -np 2 siesta

to check that it works.

Enjoy :)

[el-abed]

@pfebrer Thank you so much! It is because of the restrictions that I can not do so. But will try and see how it goes. Really laughed out so hard when you said insultingly easy ( even though I would have NEVER figured how to do that :( ) Still very funny! Thank you again!

[pfebrer]

It is because of the restrictions that I can not do so.

Have you tried? Maybe it works. In principle there should be no problems. Spack only runs python, compiles libraries and stores them in its folder.

Really laughed out so hard when you said insultingly easy

Hahahaha I meant vs having to compile and link libraries manually.

[el-abed]

Hello everyone, I can finally say with the help of our supercomputer administration that siesta 4.1.5 is installed. I have done the siesta calculations but was wondering for denchar options if those would be enough? STRUCT5textfile.txt Or is mandatory to add something like this? Denchar.MinX 0 Ang Denchar.MaxX 20 Ang Denchar.MinY 0 Ang Denchar.MaxY 10 Ang Denchar.MinZ 0 Ang Denchar.MaxZ 20 Ang

[el-abed]

Really looking forward to see if you are ok with the fdf file

[zerothi]

Hello everyone, I can finally say with the help of our supercomputer administration that siesta 4.1.5 is installed. I have done the siesta calculations but was wondering for denchar options if those would be enough? STRUCT5textfile.txt Or is mandatory to add something like this? Denchar.MinX 0 Ang Denchar.MaxX 20 Ang Denchar.MinY 0 Ang Denchar.MaxY 10 Ang Denchar.MinZ 0 Ang Denchar.MaxZ 20 Ang

Great!

I think it would be favourable to do two tests. For the orthogonal unit cells it would be useful to check that the unit-cell parameters gives the same result (i.e. Min[XYZ] 0. Ang, and Max[XYZ] equal to the respective lattice vectors). Then also try a displaced lattice. Min[XYZ] -2.45 Ang Max[XYZ] 5.6 Ang

or just some values that are not an integer multiple of the lattice vectors. In sisl this is easy since it works for skewed cells.

Hope this makes sense. But again, start with the simplest cases (orthogonal lattice vectors). If they are "wrong" then the rest are also ;)

[el-abed]

I copied the full.WFSX to make it the main WFSX file and thus Denchar is still running since i used 10 10 1 k points.

1- When you say wrong which Kx.WF.y is right or wrong? I got thousands of files? Do you want me to check the X k point or Y k point?? 2- since it is a orthogonal system i assume X Y or S or Gamma Points are what we are interested in right?

[zerothi]

I copied the full.WFSX to make it the main WFSX file and thus Denchar is still running since i used 10 10 1 k points.

1- When you say wrong which Kx.WF.y is right or wrong? I got thousands of files? Do you want me to check the X k point or Y k point?? 2- since it is a orthogonal system i assume X Y or S or Gamma Points are what we are interested in right?

I would suggest you start with the density, that seems easier. But again, we need to assert both density and wavefunction.

For k-points, any points will do. Perhaps Gamma and some other k-point is fine, a combination of X and Y would be ideal. You'll note here that the k-point you select should have a periodicity depending on the values. So select one with some integer periodicity as that is easiest.
To limit the number of files you can manually change the WaveFuncKPoints fdf block.

[el-abed]

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.