Can I use mcu to calculate the matrix element <m|d/dk n> along a custom k-space path?

[tharindufernando]

For the cell-periodic parts |m>, |n> of Bloch wavefunctions, along some custom k-space path, I want to calculate the matrix element <m|d/dk n> = <m|dH/dk|n>/Enm (where Enm = En - Em for eigenvalues En,Em). I did this using finite differences for models with analytical Hamiltonians using MATLAB, Python, etc. In these models, I have that |m> and |n> are orthonormal instantaneous energy eigenstates.
But I wish to try it using VASP as well. I would ideally get at least the wavefunction for each band, and the derivative, so that I can do a finite difference to calculate the matrix element along the k-space path. However, it has been hard for me to achieve this. VASP itself offers LPEAD and LEPSILON, but they each give me an issue. For instance, LPEAD does not support a custom k-space path; while I do not know how to read the data from LEPSILON. I have limited experience with these techniques.

Is mcu capable of giving me an array of wavefunctions and derivatives? I noticed that cder gives a derivative.

Thanks for your help!