Check why PDOS sometimes skips the first semicore, and why PDOS and DOS are different

[AndresOrtegaGuerrero]

Update November 26

• We should set the default deltaE of DOS (and PDOS) to 0.01. This will reduce the probability to encounter the problem and "fix" the case of Au as tested by Michail Minotakis.

• We should expose in advanced settings deltaE in case a user would prefer a different default.

• We should still report the issue to QE developers

BELOW info from November 19

Two notes:

• the output of QE also has the integral of the DOS/PDOS: its derivative does not miss the first state.

• For the difference between DOS and PDOS, of course they are different in the high conduction bands, but can we check why in the valence they are different? Different method, just different smearing, ...?

Some additional comments from @cpignedoli and @giovannipizzi investigation:

Issue to better understand in QE:

• PDOS not matching total dos; fast protocol, maybe still PDOS k-mesh is too coarse?
• Dos.x bug, missing narrow semicore states. E.g. when doing Au with the app, the state at -106 eV wrt Fermi level (~-85.1 with QE zero) shows no peak at all in DOS and PDOS with fast protocol, and small but not “full” peak for moderate:
• Au, moderate protocol, default QE app energy range:

  • Restarted by hand, changing Emin in the energy range to start at lower E: it cures a little bit, but still bug:

  • Note from above: the integral actually “sees” something of the peak (even if the DOS and the derivative of the integral do not match exactly, which is more accurate? Should we show by default the derivative of the integral instead? Is it possible also for the PDOS?)

• Note that in the stdout, the code “sees” the contribution. E.g. (for ref, E fermi = 20.7534 eV and with these orbitals:

    state #   1: atom   1 (Au ), wfc  1 (l=0 m= 1)
    state #   2: atom   1 (Au ), wfc  2 (l=1 m= 1)
    state #   3: atom   1 (Au ), wfc  2 (l=1 m= 2)
    state #   4: atom   1 (Au ), wfc  2 (l=1 m= 3)
    state #   5: atom   1 (Au ), wfc  3 (l=0 m= 1)
    state #   6: atom   1 (Au ), wfc  4 (l=2 m= 1)
    state #   7: atom   1 (Au ), wfc  4 (l=2 m= 2)
    state #   8: atom   1 (Au ), wfc  4 (l=2 m= 3)
    state #   9: atom   1 (Au ), wfc  4 (l=2 m= 4)
    state #  10: atom   1 (Au ), wfc  4 (l=2 m= 5)
  The lowest energy state at a couple of energy points chosen arbitrarily shows:
    k =   0.0000000000  0.0000000000  0.0000000000
  ==== e(   1) =   -85.15582 eV ==== 
    psi = 0.973*[#   1]+0.027*[#   5]
   |psi|^2 = 1.000
  
  
   k =  -0.7071067812  0.2357022604  0.0000000000
  ==== e(   1) =   -85.13980 eV ==== 
    psi = 0.992*[#   1]+0.007*[#   5]
   |psi|^2 = 1.000
  

  Therefore, it knows it’s almost 100% the first pass at: hardcoded number of electrons, added slider #1 projectors (lowest s state).

[AndresOrtegaGuerrero]

@giovannipizzi @cpignedoli: Did you use the "fast" protocol for the calculations?

According to the protocol in the PdosWorkFlow, when using the "precise" protocol, the parameter DeltaE is set to 0.01, whereas for the "fast" protocol, DeltaE is set to 0.1.

Would it make sense to create a widget in the PDOS Settings tab that allows users to set the energy_grid_step for DOS and PROJWFC calculations? This would allow users to easily visualize and adjust the values as needed.

[cpignedoli]

Thanks a lot @AndresOrtegaGuerrero. I was not aware of the setting 0.01 for precise protocol. Thus I wonder if it still makes sense or not to have the default to 0.01. Maybe we should simply warn the user that with the fast protocol... and we can suggest setting the energy_grid_step to 0.01 in the Au exercise if we agree on adding the widget. @giovannipizzi what would you suggest?

[giovannipizzi]

Since calculations are fast I would put 0.01 as a default for all, possibly comparing the Dos with even smaller value for precise (compare runtime, file size, accuracy of results). If there is time ok to add a widget in the advanced settings (otherwise put as an issue for the next round, I would eventually do it, but there are other more urgent issues to fix first)

[giovannipizzi]

The issue of PDOS and DOS being different was not solved, and changing the energy grid to skip the first semi core was just a workaround. Therefore I'm reopening this issue, with additional investigation that might point to an actual solution.

1. There are 3 implementations of tetrahedra in QE, that are selected in the NSCF step via the occupations flag (see docs):

   • tetrahedra (Tetrahedron method, Bloechl's version)
   • 'tetrahedra_lin' (Original linear tetrahedron method)
   • 'tetrahedra_opt' (Optimized tetrahedron method)

2. Internally, these are mapped into the tetra_type variable (respectively 0, 1, 2). See e.g. the source code

3. We are getting the DOS from dos.x, while the PDOS from projwfc.x (as a note: we could also get the DOS directly from the output of projwfc.x, as the second column of the aiida.pdos_tot file. @mbercx @t-reents @AndresOrtegaGuerrero @cpignedoli do we need something else from dos.x, or can we just use that value, avoiding one calculation?)

4. In the workflow, we are running with the tetrahedra occupation, and the output file of dos.x confirms this

5. In the projwfc.x calculation, when tetrahedra is specified, instead for some reason the tetrahedra_lin is used! Not sure why, but see source code, and this is also confirmed by the string printed in the output.

This gives rise to the discrepancy between the DOS. This is shown in the plot below:

Above, the total DOS from dos.x and projwfc.x disagree (while the PDOS from projwfc.x is overlapping with the DOS for low energy, where the projectors are complete for the space).

If instead we use the tetrahedra_opt flag in the NSCF, the results are much better and in much better agreement:

Here we see that the DOS from dos.x and projwfc.x perfectly overlap.

I would guess that if we specify tetrahedra_lin in the input, also then we would not have a discrepancy since then both dos.x and projwfc.x would use the same algorithm. But in the QE docs they clearly say:

'tetrahedra_lin' : Original linear tetrahedron method. To be used only as a reference; the optimized tetrahedron method is more efficient.

Proposed solution

Therefore, here is my proposed solution:

1. we should always use the tetrahedra_opt as occupation flag instead of tetrahedra (@superstar54 @t-reents where is this set, in the aiidalab app or in the aiida-qe plugin?)
2. Possibly, we do not need to run dos.x but just get the DOS as well from the output of projwfc.x (I guess it makes sense to run dos.x only if you really do not want to compute the PDOS (or in the rare cases when the pseudos do not have projectors?)). But I do not know if this is ever the case in the AiiDAlab QE app?

I suggest to go ahead with 1 now, and possibly evaluate 2 with lower priority (priority P1 in the QE app demo project).

Note: I didn't check other details, e.g. the actual PDOS channels being lower than the total DOS etc, but the results above make me confident that my suggested solution will solve the problem.

Moreover, the tetrahedra_opt does see a peak for the semi core flat band, so this also seems to solve that issue (at least in my simple test)

[giovannipizzi]

A few notes:

• I think the main change to do is here, replacing tetrahedra with tetrahedra_opt in the YAML file of the PDOS protocol
• One has to fix the check here to at least support all tetrahedra flags (probably one can use 'startswith' instead?). Example change:

      if not parameters.get('SYSTEM', {}).get('occupations', '').startswith('tetrahedra'):
       return '`SYSTEM.occupations` in `nscf.pw.parameters` is not set to one of the `tetrahedra` options.'

• good to double check that the parser works correctly with the new flags, as probably this was never tested. See e.g. this comment in the legacy parser. Is this still used? Probably not, but good to check if the parsed information is correct

[giovannipizzi]

As a follow up, I quickly did the 2 changes above (YAML + fix validation) and tried the app again, with gold and fast protocol. Considering it's the fast protocol, I think results are very good and fix this issue.

Lowest semicore:

We now see the peak and its character. To improve, probably it would be good to start not the the min energy but maybe, say, 0.1 eV below. (Note: using more k-points does not help for this, I tried with a k-point linear density of 0.1 for the DOS as in the moderate protocol and the plot looks similar.

Good, Now the PDOS is always below the DOS and the sum of the PDOS channels gives the DOS (at least at low energy where orbitals can fully represent the states)

[cpignedoli]

great @giovannipizzi . For your point 2) in "proposed solution" I am not an expert in DOS, I always focused on bandstructures in my calculations. I would directly check with Nicola about what should be done. It would be great if we could skip the dos.x calculation.

[Minotakm]

Great @giovannipizzi, I will make the changes you suggested in the aiida-quantumespresso plugin and test it.

[AndresOrtegaGuerrero]

@giovannipizzi , @t-reents there is also this issue for projwfc.x in aiida-quantum espresso , aiidateam/aiida-quantumespresso#1041 , if dos.x calculation is removed we need to get the total_pdos for magnetic systems to be properly parsed

[giovannipizzi]

I commented there, if that's the only blocking step, let's just fix the parser

[AndresOrtegaGuerrero]

I was wondering @superstar54 , should we set tetrahedral_opt for projwfc.x and dos.x calculations in aiidalab-qe as a provisional solution while we wait for the changes in aiida-qe plugin ?

[edan-bainglass]

@AndresOrtegaGuerrero We're looking at this at the moment, but was this actually done? Can't remember why I closed it.

[AndresOrtegaGuerrero]

I didnt work on this, i was not around. I still have the same question of what is the decision to solve this @giovannipizzi. "should we set tetrahedral_opt for projwfc.x and dos.x calculations in aiidalab-qe as a provisional solution while we wait for the changes in aiida-qe plugin ?"

[giovannipizzi]

I think that, if there is a workaround in aiidalab qe for now, we should do it for this release, so yes, let's use the opt version for tetrahedra in the app. Let's also open another issue not to forget to do the actual fix, afterwards. At least we avoid getting Dos and pdos that don't match

[AndresOrtegaGuerrero]

ok , so to clarify, @giovannipizzi we will set 'tetrahedra_opt' in the nscf calculation occupations='tetrahedra_opt' , and in the dos.x bz_sum='tetrahedra_opt' , right?

[giovannipizzi]

Yes, but please double check that those two changes indeed produce the expected result (try gold with the fast protocol, i.e. the in app guide that edan prepared, since it will be the first thing that everybody will do and at least that should work)

[AndresOrtegaGuerrero]

@giovannipizzi , i cannot do this change at the aiidalab qe level, i tried but to set these options there must be PR in quantum espresso plugin , since

nscf doesnt allow occupations different than tetrahedra

def validate_nscf(value, _):
    """Validate the nscf parameters."""
    parameters = value['pw']['parameters'].get_dict()
    if parameters.get('CONTROL', {}).get('calculation', 'scf') != 'nscf':
        return '`CONTOL.calculation` in `nscf.pw.parameters` is not set to `nscf`.'
    if parameters.get('SYSTEM', {}).get('occupations', None) != 'tetrahedra':
        return '`SYSTEM.occupations` in `nscf.pw.parameters` is not set to `tetrahedra`.'

and for dos calculations it doesn allow bz_sum since the get_parameter_schema doesnt allow bz_sum

def get_parameter_schema():
    """Return the ``PdosWorkChain`` input parameter schema."""
    return {
        '$schema': 'http://json-schema.org/draft-07/schema',
        'type': 'object',
        'required': ['DeltaE'],
        'additionalProperties': False,
        'properties': {
            'Emin': {
                'description': 'min energy (eV) for DOS plot',
                'type': 'number'
            },
            'Emax': {
                'description': 'max energy (eV) for DOS plot',
                'type': 'number'
            },
            'DeltaE': {
                'description': 'energy grid step (eV)',
                'type': 'number',
                'minimum': 0
            },
            'ngauss': {
                'description': 'Type of gaussian broadening.',
                'type': 'integer',
                'enum': [0, 1, -1, -99]
            },
            'degauss': {
                'description': 'gaussian broadening, Ry (not eV!)',
                'type': 'number',
                'minimum': 0
            },
        }
    }

def validate_dos(value, _):
    """Validate DOS parameters.

    - shared: Emin | Emax | DeltaE
    - dos.x only: ngauss | degauss | bz_sum
    - projwfc.x only: ngauss | degauss | pawproj | n_proj_boxes | irmin(3,n_proj_boxes) | irmax(3,n_proj_boxes)

    """
    jsonschema.validate(value['parameters'].get_dict()['DOS'], get_parameter_schema())