Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Multiple DOS overlapping #154

Open
BlessedAkarat opened this issue Jul 15, 2024 · 17 comments
Open

Multiple DOS overlapping #154

BlessedAkarat opened this issue Jul 15, 2024 · 17 comments

Comments

@BlessedAkarat
Copy link

BlessedAkarat commented Jul 15, 2024 •
edited
Loading

Dear pyprocar users.

Hello.
I'm trying to draw some Density of States by using pyprocar from Quantum ESPRESSO data, and here I have some questions.

1) At first, I draw some DOS of primitive Ruthenium with plain mode.
And I want to draw with axis units together.
image

I already know the x-axis represent Energy in (eV) units, but I'm not sure about the units for the y-axis of DOS.
Generally, the units for DOS are likely to be (states/eV), but in some examples, the units are not specified or are written as (a.u.) (possibly denoting Bohr radius).

Especially, in "User Guid", plotting DOS with VASP data in 'plain' mode denotes (a.u.) for the y-axis unit.
(user guid: https://romerogroup.github.io/pyprocar/user-guide/dos.html#)
But with QE, in my case, the figure came out with out written any units for y-axis. (watch above the picture)
How can I check the units of y-axis?

2) Second, Can anyone please tell me how to normalize y-axis with slab size?
When I try to calculate slab DOS, y axis seems not to be normalized.
But, I wan't to analyze these data with normalized same scale.
Please tell me how.

3) Last, how can I compare normalized DOS for slabs of different sizes effectively in a single Figure?

When slab size differs, DOS should be differs too.
In fact, there should be some difference between each DOS of each slab size and I want to observe them.

Additional question.
When plotting DOS with QE datas in plain mode, should I set the Fermi level on execution .py file?
It seems there is nothing change whether it distinguished or not.
<Dos.py>
image

But in the "User Guid", I linked above, they set the Fermi level for the option.
(Moreover, that output figure is not similar form with this QE output form, even if I didn't consider the spin.)
<User Guid, VASP, Plain>
image
How to set up the DOS plotting codes with QE outputs?

  • DOS for 2 atom slab (Ru primitive + 15 Angstrom vacuum)
    image

  • DOS for 10 atom slab (x5 super cell with z direction of Ru (001) + 15 Angstrom vacuum)
    image

Best regards,
G. H.
@ahromero
Copy link
Collaborator

ahromero commented Jul 15, 2024 via email

@lllangWV
Copy link
Member

Hey just to add to the conversation

I already know the x-axis represent Energy in (eV) units, but I'm not sure about the units for the y-axis of DOS.
Generally, the units for DOS are likely to be (states/eV), but in some examples, the units are not specified or are written as (a.u.) (possibly denoting Bohr radius).

ahromero is correct we just keep the units of the underlying DFT code. If you want to change the label, you can, add the y_label

pyprocar.dosplot(
                    code='vasp',
                    mode='plain',
                    orientation='horizontal',
                    fermi=5.3017,
                    grid=True,
                    y_label='DOS [a.u.]'
                    dirname=data_dir)

Second, Can anyone please tell me how to normalize y-axis with slab size?

Currently, there is not an option to do this in the package. But, this sounds interesting we can add the option to choose the normalization, either by the total number electrons like ahromero mentioned or by the max value. I will update you when i have done this.

Last, how can I compare normalized DOS for slabs of different sizes effectively in a single Figure?

You can compare different density of sates in the following way. dosplot returns matplotlib fig and ax object, so you can pass the ax object to dosplot for another density of states.

fig,ax =pyprocar.dosplot(
                    code='vasp',
                    mode='parametric_line',
                    orbitals = [1],
                    atoms=[2,3,4],
                    spins=[0],
                    clim=[0,1],
                    orientation='horizontal',
                    fermi=5.3017,
                    grid=True,
                    show=False,  # make sure show is False to avoid producing 2 plots
                    dirname=data_dir)


pyprocar.dosplot(
                    code='vasp',
                    mode='overlay_species',
                    orbitals = [4,5,6,7,8],
                    clim=[0,1],
                    orientation='horizontal',
                    fermi=5.3017,
                    grid=True,
                    ax=ax,   # Pass the ax object to dosplot
                    dirname=data_dir,
                    plot_total=False, # Do not show total as the previous plot has the same total
                    )

When plotting DOS with QE datas in plain mode, should I set the Fermi level on execution .py file?
It seems there is nothing change whether it distinguished or not.

Yes, set the fermi energy. Otherwise the energy axis will not be shifted by the fermi energy.

(Moreover, that output figure is not similar form with this QE output form, even if I didn't consider the spin.)

I apologies about this, this image is out of data we need to update to the current configuration.

Logan Lang

@BlessedAkarat
Copy link
Author

Dear Ahromero and Logan.

Thanks for both your help, and sorry about my late answer.
Because of my lack of understanding, I have to take some time before I ask more.

Well I have to check out my self about what you replied me and answer again.

Before, I wan't to know something about the y-axis unit.
You've told me the unit of y-axis is defined what DFT output unit is.
And there I have some output which has LDOS and PDOS data.

<Ru.k.pdos_atm#10(Ru)_wfc#>
image

This seems like y-axis refers to (E), but I have no idea about this unit.
Typically, the DOS is calculated as a the function of energy, with outputs represented in (states/eV)... as I understand :(.

Last, there seems to be no change in the DOS plot after varying the Fermi level in the code. (I'm using qe option).
What kind of role does this play?

  • Dos.py, Fermi level for Slab Ru (1.3074)
    image
    image

  • Dos.py, Random Fermi level
    image
    image

@lllangWV
Copy link
Member

Hey,

This seems like y-axis refers to (E), but I have no idea about this unit.
Typically, the DOS is calculated as a the function of energy, with outputs represented in (states/eV)... as I understand

You are absolutely correct in your assumption. Here, from there documentation for projwfc.x

All DOS(E) are in states/eV plotted vs E in eV

Last, there seems to be no change in the DOS plot after varying the Fermi level in the code. (I'm using qe option).
What kind of role does this play?

Hmmmm, this does not happen in my version. I am pushing updates to pypi tonight. I will message you when I have. Install those and try this again.

It should shift dos plot such that the fermi level is at 0. Also it will change the x axis label.

image

Logan Lang

@BlessedAkarat
Copy link
Author

Yes, you are right.

If the Fermi level option works, the DOS plot should shift accordingly.

But as I show above, there is nothing different after Fermi changed.
I'm really happy that you made some patches for this.
So, does that mean all the DOS I've drawn so far automatically have the Fermi level set, or are there incorrect values drawn? Is there any way to know what criteria were used to calculate the Fermi level?

  • Even though I had already seen the link for projwfc.x, for finding the unit of DOS, I hadn’t checked out the section about the DOS unit. If I had just scrolled down a bit more, it would’ve explained that in the output, but stupidly, I thought the explanation for the input was all there was and didn’t scroll down further.
    Thanks for bringing it up again!

Always appreciate your help.
Best regards,
G. H.

@lllangWV
Copy link
Member

Hey,

I pushed the version 6.2.1 to pypi, so you should be able to [pip install pyprocar.

So, does that mean all the DOS I've drawn so far automatically have the Fermi level set, or are there incorrect values drawn? Is there any way to know what criteria were used to calculate the Fermi level?

In a previous versions, we automatically shifted the fermi energy in the code. We decided to move away from this because other DFT package do not retain the self-consistent Fermi energy in the same directory as the band calculation. Quantum Espresso is able to do this as the {Prefix}.xml keeps the self-consistent fermi energy, but for consistency across the other DFT codes we choose not to shift the energy unless fermi is specified by the user. Though, we are considering the automatic shift back to DFT codes that retain the self-consistent fermi energy. You can also find the fermi energy in the output of the self-consistent calculation.

The change in behavior to the fermi energy happened in version 6.1.9 Mar 28th, 2024. Quantum Espresso never really had an issue with its fermi energy, so I would say the dos you've drawn should be correct.

Yeah, I've looked at their documentation many times and I always seem to find new details I didn't see before. It's just how it goes.

I'm glad I could help!

Best,

Logan Lang

@BlessedAkarat
Copy link
Author

fermi=20.3074
image

fermi=1.2854
image

Yeah, it works perfectly.
However, let me ask something more.
In my output folder, my final output .xml is overwritten non-self consistent (nscf) or bands calculation data.

Currently, the calculations are being performed sequentially, with SCF and NSCF being carried out together as a single task.
Because I need nscf output data for the BoltzTraP2 interpolation.
The .xml file will contain NSCF data, and I should use "fermi" for the NSCF Fermi level. Is that correct?

grep Fermi *out

image

It's a bit disappointing that, like you said, some tools can not directly access the Fermi level in the output data. It's actually more convenient when the Fermi level is automatically set like before (nonetheless, thanks for your consideration). In that sense, considering going back to the original setup seems like a good idea, but I think it wouldn't hurt to have a way to shift the Fermi level when you do input values, and only set it to default if you don't.

Also, is there any updates for normalization you told me above?

Currently, there is not an option to do this in the package. But, this sounds interesting we can add the option to choose the normalization, either by the total number electrons like ahromero mentioned or by the max value. I will update you when i have done this.

Thanks for your help.

Best regards,
G. H.

@lllangWV
Copy link
Member

The .xml file will contain NSCF data, and I should use "fermi" for the NSCF Fermi level. Is that correct?

You are right that for nscf mode over writes the data, including the fermi energy. However, for bands mode the fermi energy level is left unchanged. You should use the self-consistency fermi energy, though the difference between scf and nscf will not be large since they both use dense kmeshes.

Either way I do not like the consistency, so I changed to parsing of the fermi energy to scf.out

It's a bit disappointing that, like you said, some tools can not directly access the Fermi level in the output data.

I agree, I changed it so qe the shift by the fermi energy will automatically be applied.

Also, is there any updates for normalization you told me above?

I just added this feature in. You can pass normalize_dos_mode to dosplot. You can either normalize by the normalize_dos_mode=max or normalize_dos_mode=integral.

Here is an example:

fig,ax = pyprocar.dosplot(
                    code=code,
                    mode='plain',
                    orientation='horizontal',
                    normalize_dos_mode='max',
                    grid=True,
                    show=True,
                    dirname=data_dir)

I pushed the changes to the main branch of github, so you will have to install from the github repository.

Best,

Logan Lang

@BlessedAkarat
Copy link
Author

Dear Logan,
You work really hard.

I'm very appreciate every work you did.
However, I don't mean to bother you , I have few more questions...

First, what about using only 15x15x15 for bulk SCF, 50x50x50 for bulk NSCF and 15x15x1 for slab SCF, 50x50x4 for Slab NSCF?
You mentioned that I should use the SCF Fermi energy rather than NSCF, but if the k-meshes are dense enough, NSCF would be acceptable since there is not quite big difference.
In my case, are those k-points dense enough?
They usually have 0.03~0.1 eV differences between SCF and NSCF.

And, please let me know where can I look up the options about normalize_dos_mode?
I can't find out any differences between 'max' and 'integral'.
Also, I want to check out what kind of normalization each are.

  • Ru Primitive
    image
    image

image
image

Last, I can't understand your last sentence,

I pushed the changes to the main branch of gitthub, so you will have to install from the github repository.

I apologize for my bad knowledge. I'm not very familiar with github and python... so I have no idea what I have to do.
Is that mean, you already update this job in the 6.2.0 version?

Sincery,
G. H.

@ahromero
Copy link
Collaborator

ahromero commented Jul 25, 2024 via email

@BlessedAkarat
Copy link
Author

Dear Ahromero.

I think you have to answer these questions by reviewing what is DFT and how it is done…. There is a difference between what you do in SCF and NSCF…. The point is not the K-mesh!!!

First of all, as I understand it, the SCF calculation is used to converge the electron density function through self-consistent iteration. Since a lot of computational cost is required for the iteration, I understand that specifying many k-points through SCF is a challenging task.
On the other hand, NSCF is the process of using the electron density derived from SCF to calculate electronic band structures with a denser k-mesh.

So, when you asks why I use the non-self-consistent (NSCF) if the self-consistent (SCF) is the correct approach, I think it means I'm misunderstanding the DFT system. I thought that using SCF first to calculate a rough structure and then get a more detailed structure with NSCF.
That's why I believed the YouTube examples using NSCF data for interpolation when running calculations with a program called BoltzTraP2 were showing the right method.
If I calculate NSCF with a sufficiently dense k-mesh, shouldn't it represent the results from the SCF calculations?

What do you mean ty “find out”… did you check the numbers? The proper way to normalize DOS is with the number of electrons, see Kittel or Cardona’s book… you can always normalize to one for particular comparisons but they have to be taken with care.. for example to compare supercell vs unit cell calculations

Well, when I run the DOS plot through that option, I get the graph data, right? But I never thought about changing it into numbers to take a look at it...
I’m saying this because you asked if I checked the numbers, Can I recheck these data (Post-processing data from the outputs obtained through projectwfc.x calculations) in to numbers?
And even if we quantify the data with numbers, if there's already no change in the figures that reflect that data, doesn't it make a difference?

I’m sorry if my question has made you feel a bit frustrated or angry.
But, I’m really trying my best to understand what I can, so I’d like to ask for your understanding on this part.

Thanks.
G. H.

@ahromero
Copy link
Collaborator

ahromero commented Jul 25, 2024 via email

@BlessedAkarat
Copy link
Author

Actually, the lab I belong isn't that big and it's a relatively new lab (just opened 2 years ago).
So there are very few students working on the same topic (just two of us, including me), and he also don't know much about DFT.
This is why I mainly gather information from places like YouTube or Google, but of course, I also try to find supporting materials like papers to verify their authenticity.

My advisor said that assigning a large k-point for SCF calculations uses too many computational resources, so it's better to use NSCF instead. I actually ran a few convergence tests with that approach. (If I remember correctly, for Ru, NSCF over 30x30x30 was sufficient for convergence.) That’s why I’m using NSCF for the current calculations, and the pyprocar homepage also suggests doing DOS plots in the order of SCF > NSCF > PDOS, so I was following that process.
However, there are opinions that using SCF is more accurate, while NSCF can give different values for some systems (especially semiconductors or semi-metals), so I think I should discuss that part again with my advisor. (Though all the systems I’m calculating are conductors...)

Well, thanks for your guidance.
I'll check this again.

And even if we quantify the data with numbers, if there's already no change in the figures that reflect that data, doesn't it make a difference?

This means that you asked me to check the numbers, but if there’s no difference between the two figures, then there shouldn’t be a noticeable difference in the numbers themselves, right?
I think I'm not really getting what you're trying to point out either.

Best regards,
G. H.

@ahromero
Copy link
Collaborator

ahromero commented Jul 27, 2024 via email

@BlessedAkarat
Copy link
Author

BlessedAkarat commented Aug 1, 2024 •
edited
Loading

I talked to my advisor again, and we concluded that using the Fermi level of SCF is indeed the right approach.

Lastly, I’m still pretty clueless about the normalize_dos_mode part. I think I need to think about it some more.

Anyway, thanks for your help!

@lllangWV
Copy link
Member

lllangWV commented Sep 9, 2024

Hey,

Apologies for the delayed response; my graduate school workload has been increasing.

To clarify, normalizing the density of states involves dividing the dataset by a specific reference value. There are two main approaches for choosing this reference: either by using the maximum value of the density of states (normalize_dos_mode=max), or by taking the integral of the total density of states with respect to energy (normalize_dos_mode=integral), which reflects the total number of states.

When you set the normalize_dos_mode, I adjust the dataset accordingly. Most people prefer normalizing by the maximum value, as it allows for easier comparison across systems. This approach also lends itself to a pseudo-probability interpretation, even though it’s not a true probability distribution since it doesn’t integrate to 1. On the other hand, using the integral mode makes the density of states more directly comparable to a probability distribution, as the integral is forced to equal 1.

Best regards,
Logan Lang

@BlessedAkarat
Copy link
Author

I really thought I had seen your comment and replied about it three weeks ago,
but I guess I didn’t hit the comment button after rewriting the post.

I apologize for unintentionally disregarding your sincere efforts.

Anyway, over the past few weeks, I've been so busy with urgent stuff that I haven't been able to properly check the bands and DOS yet.
I feel like my knowledge isn’t enough, so I still don’t fully understand it, but I'll think it over again and come back if there are still issues.

Thanks for all your help in many ways.

Sincery,
G. H.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@ahromero @lllangWV @BlessedAkarat