For pedagogical reasons, I am looking for ways to calculate quantum-mechanical quantities such as $\langle m | \dot{m}\rangle, \langle m | \dot{n}\rangle, \langle m | \ddot{n}\rangle$ using wavefunctions $m,n$ output from DFT calculations (Quantum Espresso, specifically). The overhead dots denote derivatives. Ideally, I would be able to calculate the Berry phase for custom loops as well.

So far, I have tried using some Python tools (z2Pack, PythTB) to try and post-process wannier90 output, but as a beginner, it doesn't seem as if these tools allow us to calculate our own matrix elements. They seem to be just for intra-level quantities such as $\langle m | \dot{m}\rangle$.

Does anyone more-experienced have any advice on what existing software might be easiest to work with to calculate these custom-defined quantities? I would rather try messing with some software that is more likely to support these kinds of calculations. PythTB seems to be one of the best options I have, but Python may not be the best for large datasets (compared to Fortran, etc). Thank you for your time.

Edit: The issue also seems to be that several post-processing software tend to use the numerical method where one takes $\arg$ of some product of complex phases corresponding to each $k$ point. However to do something like $\langle m | \dot{n}\rangle$, it might be better to use a central difference method to carry out the derivative of the wavefunction (as opposed to an established discretized method for the Berry phase, as in section 4.5 in Ref [1]).

[1]: Tight-Binding Formalism in the Context of the PythTB Package, https://www.physics.rutgers.edu/pythtb/formalism.html

  • $\begingroup$ I gave my +1 to this a long time back! But @TribalChief I just wanted to let you know that your question has been mentioned here: mattermodeling.stackexchange.com/q/6422/5 Do you think you're able to help that (new) user with their first question on the entire Stack Exchange network? $\endgroup$ Commented Aug 9, 2021 at 16:59
  • 1
    $\begingroup$ @NikeDattani, thanks for bringing this to my attention. I was able to make a MATLAB implementation a while back. I am away for a few days but can get to the question sometime Thursday. The user will probably have to process things further from there using Python. $\endgroup$ Commented Aug 9, 2021 at 17:15
  • $\begingroup$ Beautiful! I think a MATLAB implementation would be good enough, especially if it works in Octave! I have more experience with MATLAB than Python, so perhaps I could help the user if they need further help after your answer. $\endgroup$ Commented Aug 9, 2021 at 17:20

3 Answers 3


You can use QE.6xxx with the support of the hdf5 library. To realize that purpose, you should add the following command when you compile QE:


or take a look at the official guide.

Then the saved wavefunction can be manipulated with the python package h5py.

import h5py
  • $\begingroup$ Thanks for the answer. Would you mind confirming that if I compile QE with that flag, QE will automatically save hdf5 files (with wf_collect = .true. in pw.x)? I ask because I see only wfc1.dat. I tried looking in the documentation for more on this, but all I see is how to enable hdf5. I just want to make sure that's all that is needed to be done. Additionally, do you know of any references/links that demonstrate manipulation of these QE hdf5 files? Perhaps I am looking at the wrong documentation... $\endgroup$ Commented Jan 25, 2021 at 2:01
  • $\begingroup$ @TribalChief You should add that tag in your configure file. Take a look section 2.3 in this link: quantum-espresso.org/Doc/user_guide.pdf $\endgroup$
    – Jack
    Commented Jan 25, 2021 at 7:14

In sisl one can do the derivatives if one has access to the Hamiltonian.

The sisl code is heavily dependent on the LCAO formalism. So if you can't work with Hamiltonians it might be more difficult (read implement your own derivative function).

One can e.g. do something like this:

import sisl as si
# Read in a Siesta Hamiltonian
H = si.get_sile("RUN.fdf").read_hamiltonian()
eig = H.eigenstate(k=[0.5, 0.25, 0])
print(eig.derivative(order=2, matrix=True))

the above will read in a Hamiltonian from a Siesta calculation, then diagonalize at the specified k-point, then calculate the 1st and 2nd order derivatives and return the off-diagonals as well. The derivative function is for instance used to calculate the velocities, and one can overlay these one the band-structure plot (plus numerical derivatives) to assert they yield the correct quantities.

  • $\begingroup$ Title changed after my answer, so probably this answer does not answer anything :) $\endgroup$
    – nickpapior
    Commented Mar 7 at 20:41

You may be interested in DFTK, which is a fully featured plane-wave code like Quantum Espresso but which fits in a mere 7k lines of code. Importantly, DFTK has automatic differentiation, which means you should be able to compute whatever derivative you wish for in a numerically accurate fashion.


You must log in to answer this question.

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