Norm-conserving pseudopotentials are defined such that the energy derivative of the log derivative of the real and pseudo wavefunctions agree at $r>r_c$. I understand that although valence eigenenergies are discretised, they (and their associated valence wavefunctions) can be varied continuously by applying a perturbation to the potential, which allows us to define the energy derivative of the log derivative. However, I do not understand which perturbation we apply in this definition. I asked a related question here.

  • 1
    $\begingroup$ I gave my +1 long ago, but I wonder if you've had any luck with this since it's now been almost 14 months since you asked the question? Did you figure out the answer? Are you still working on it? $\endgroup$ Commented Jan 31, 2023 at 21:10
  • $\begingroup$ I think that the log derivative is defined on all solutions to the Schrodinger equation (with a continuous range of energy values), not only on valence eigenfunctions for which energy discretised (since $|\psi(r)|\rightarrow 0$ as $r\rightarrow \infty$). If I'm right about this, we don't consider a specific perturbation, rather the aforementioned continuous set of solutions to the unperturbed Schrodinger equation. Do you know whether I'm right about this? $\endgroup$ Commented Feb 3, 2023 at 15:04

1 Answer 1


You should not be thinking about a specific perturbation operator. The key thing here is that the orbitals are Kohn-Sham eigenfunctions, and that already reveals everything that is needed. Still, it took me a while to dig out the literature on how these orbital energy derivatives are calculated. However, eventually I landed on the paper Use of energy derivative of the radial solution in an augmented plane wave method: application to copper by Koelling and Arbman, which explains how.

Shortly put, most atomic codes use finite differences for solving the radial orbitals. (Mine use finite elements instead, see Int. J. Quantum Chem. 119, e25945 (2019) and the newest preprint on arXiv:2302.00440.) If you use finite differences, the radial functions $u_{nl}(r)$ are found from the equations

\begin{eqnarray} \hat{h}_{nl}u_{nl}(r) = & E_{nl} u_{nl}(r) \\ \hat{h}_{nl} =& -\frac{1}{r} \frac{{\rm d}^2} {{\rm d}r^2} r + \frac {l(l+1)} {r^2} + V(r) \end{eqnarray}

Koelling and Arbman's trick is to take the energy derivative $\partial / \partial E_{nl} = \dot{}$ from the first equation, yielding

\begin{eqnarray} (\hat{h}_{nl} -E_{nl}) \dot{u}_{nl}(r) = & u_{nl}(r) \end{eqnarray}

which is an inhomogeneous differential equation with $u_{nl}$ as the homogeneous solution. Note that you already know $u_{nl}(r)$ and $E_{nl}$ here! Koellig and Arbman therefore simply solve this differential equation for $\dot{u}_{nl}(r)$.

Next, the requirement that the wave functions are orthonormal

$$ \int_0^\infty {\rm d}r r^2 u_{nl}(r) u_{n'l}(r) = \delta_{nn'} $$

gives you the condition

$$ \int_0^\infty {\rm d}r r^2 [\dot{u}_{nl}(r) u_{n'l}(r) + u_{nl}(r) \dot{u}_{n'l}(r)] = 0 $$

which means that $u_{nl}$ and $\dot{u}_{nl}$ are supposed to be orthogonal. Koellig and Arbman simply orthogonalize the obtained numerical solution for $\dot{u}_{nl}(r)$ against $u_{nl}(r)$.


The finite difference schemes traditionally used in numerical atomic orbital calculations are not based on the spectral theorem, and instead rely on a number of tricks such as the shooting method to solve the self-consistent field equations. I think that the orbital energy derivative scheme is inherently tied to the solution method, which makes the scheme mathematically ill-defined.

The problem is that we know that the spectrum of the Hamiltonian is discrete, and that its solutions form a complete set of states. This means that one should be able to expand the energy derivative wave functions in terms of the original wave functions, $|\dot{\psi}_i\rangle = |\psi_k\rangle R_{ki}$. However, it is easy to see that the equation used to define the energy derivative has no solutions.

It is already clear from the definition of the Kohn-Sham-Fock matrix being diagonal in the canonical orbital basis $\hat{f} = \text{diag }{\boldsymbol{\epsilon}} $ that the only perturbation that can affect the energy of the $i$th orbital is a level shift $\Delta {\bf H}^{i} = s | i \rangle \langle i |$, which is an eigenoperator of the orbital $i$ and therefore does not result in a changed operator. I would be interested if anyone can prove me wrong on this point...


You must log in to answer this question.

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