Using SIESTA, I want to plot the wavefunctions around the Fermi level (similar to HOMO/LUMO). Also, I want to do COOP analysis.

From the SIESTA manual, I can setup to write the wavefunctions associated to a set of bands or selecting them by energy:

The user can optionally request that the wavefunctions corresponding to the computed bands be written to file. They are written to the SystemLabel.bands.WFSX file.


The user can optionally request that specific wavefunctions are written to file. These wavefunctions are re-computed after the geometry loop (if any) finishes, using the last (presumably converged) density matrix produced during the last self-consistent field loop (after a final mixing). They are written to the SystemLabel.selected.WFSX file.

Note that the complete set of wavefunctions obtained during the last iteration of the SCF loop will be written to SystemLabel.fullBZ.WFSX if the COOP.Write option is in effect.

I am using the following keywords:

WFS.Write.For.Bands     T
COOP.Write              T
WriteEigenvalues        T
WriteKbands             T
WriteBands              T
WriteWaveFunctions      T

The results are two files a SystemLabel.fullBZ.WFSX with ~800MB and a SystemLabel.bands.WFSX with 14GB.

A doubt:
Which one should I use to plot the wavefunctions corresponding to eigenvalues around the Fermi energy?

Thanks in advance.

  • $\begingroup$ Is the question answered? If so, could you tick correct to make it answered? :) $\endgroup$
    – nickpapior
    Commented Feb 24, 2022 at 11:06

3 Answers 3


The two different files contain different things:

  • *.fullBZ.WFSX contains all wavefunctions for the SCF sampled Brillouin zone, i.e. for all k-points defined in the kgrid.MonkhorstPack block. The k-point order can manually be inferred from the KP file.
  • *bands.WFSX contains the wavefunctions for the bandstructure sampled Brillouin zone, i.e. from the blocks BandLines or BandPoints.

Since you are interested in the Fermi-level and that the band-structure contains the Fermi-level this would be easier to deal with, but both contains the relevant information in this case.

I would however suggest you used sisl (disclaimer, I am the author) for reading wavefunction coefficients and coop analysis.

Something like this would suffice:

import numpy as np
import sisl
H = sisl.get_sile("RUN.fdf").read_hamiltonian()
wfsx = sisl.get_sile("*.fullBZ.WFSX", parent=H)

E = np.linspace(-2, 2, 201)
for es in wfsx.yield_eigenstate():

which will print COOP curves for every eigenstate contained in the WFSX file. Then you can selectively decide which ones you want.

  • $\begingroup$ Running you code, returns the following error: Traceback (most recent call last): File "/home/icamps/Quantum_Chemistry/0_BN_new/SIESTA_outputs/BN/del.py", line 7, in <module> for es in wfsx.yield_eigenstate(): $\endgroup$
    – Camps
    Commented Nov 29, 2021 at 23:52
  • $\begingroup$ File "/home/icamps/.local/lib/python3.9/site-packages/sisl/io/siesta/binaries.py", line 1045, in yield_eigenstate _bin_check(self, 'yield_eigenstate', 'could not read sizes.') File "/home/icamps/.local/lib/python3.9/site-packages/sisl/io/siesta/binaries.py", line 46, in _bin_check raise SileError(f'{str(obj)}.{method} {message} (ierr={ierr})') sisl.io.SileError: wfsxSileSiesta(*.fullBZ.WFSX, base=/home/icamps/Quantum_Chemistry/0_BN_new/SIESTA_outputs/BN).yield_eigenstate could not read sizes. (ierr=2) $\endgroup$
    – Camps
    Commented Nov 29, 2021 at 23:53
  • $\begingroup$ I think you put in the literal '*.fullBZ.WFSX' it should be the full filename. If you still have problems, could you please open an issue here: www.github.com/zerothi/sisl/issues $\endgroup$
    – nickpapior
    Commented Nov 30, 2021 at 21:09

I believe both can technically be used. Let me first point out that your calculation must have been long and heavy ! Good luck with those. Your fullBZ.WFSX files are used for the COOP bonding analysis you want to do and .bands.WFSX will give you your fat bands which you can use for your wave-functions around the Fermi.

Here's a tutorial from SIESTA that might help.


An alternative way to explore the states around a given energy is to compute the "local density of states", which is a kind of "charge density" obtained from the states in a given energy range.

If you want wave-functions, you can select a more interesting set by using energy filters. From the manual:

The user can narrow the energy-range used (and save some file space) by using the WFS.Energy.Min and WFS.Energy.Max options (both take an energy (with units) as extra argument), and/or the WFS.Band.Min and WFS.Band.Max options. Care should be taken to make sure that the actual values of the options make sense.

Note that the band range options could also affect the output of wave-functions associated to bands (see section 6.15.2), and that the energy range options could also affect the output of user-selected wave-functions with the WaveFuncKPoints block (see section 6.16).


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .