The Hamiltonian for fine structure (the atom with $\text{Z}$ protons and with electron interaction terms included) is $$H=\frac{\text{Z}^2}{ r}+\underbrace{\frac{p^2}{m}+\frac{p^4}{m^3}}_{\text{kinetic}}+\underbrace{\frac{\text{Z} \ L\cdot S}{r^3}}_{\text{spin-orbit}}+\underbrace{\frac{\text{Z}}{m^2}\delta(r)}_{\text{Darwin term}}$$ modulo constants in from of each summand.
Apparently there is a derivation of this using the Dirac equation. Could anyone give a link to this?
You can find the complete derivation in ref.1, by using the Dirac equation. You may want to complement it by having a look at ref.2. first, where the Dirac equation is derived from the principles of quantum electrodynamics (which is a more fundamental theory), thus obtaining a more complete picture.
Sketch: ref.2. takes the QED Lagrangian, with operator $\hat\psi(x)$, and derives the Dirac equation for the field $\psi_n(x):=\langle 0|\hat\psi(x)|n\rangle$, where $|n\rangle$ is a complete set of state-vectors and $|0\rangle$ is a vacuum state. Given the Dirac equation for $\psi_n(x)$, ref.1. performs the so-called Foldy-Wouthuysen transformation, and then takes the non-relativistic limit, thus obtaining the first-order corrected Schrödinger equation. All in all, refs.1,2 provide a first-principles derivation of the Hamiltonian OP is after (up to factors of $Z$, which are easy to reintroduce).
References.
Itzykson and Zuber, QFT, §2-2-4.
Weinberg, QFT, Vol.I, §14.1.
Either you are mixing something up or the book explains it poorly. What you call the electron interaction term (I suppose) is the usual Coulomb term describing the interaction between the proton and the electron,
$$H_{\text{Coulomb}} = -\frac{e^2}{4 \pi r}$$
The second term corresponds to the kinetic energy of the electron (note that we assume the proton to remain stationary). However, as you may know $ E = mv^2/2$ is a classical expression and the more complete expression is the relativistic one (which reduces to $ E = mv^2/2$ for $mc^2 >> pc$ which is exactly the first order in $p^2c^2/m^2$ in the Taylor expansion of E). So the third term that you wrote down is called the relativistic correction to the kinetic energy and has nothing to do with electron interactions.
The corrections that are jointly called fine-structure contain apart from the relativistic correction also the spin-orbit coupling and the Darwin term. They both correct in some sort the assumption that the proton remains stationary. All this is quite well explained on the Wikipedia page https://en.wikipedia.org/wiki/Fine_structure
We can obtain the Fine structure Hamiltonian from an expansion on the ratio between the Bohr Radius $a_0=\frac{4\pi\varepsilon_0 \hbar^2}{m e^2}$ and the Compton Wavelength $\lambda_0=\frac{\hbar}{mc}$, i.e. an expansion on $\alpha=\frac{\lambda_0}{a_0}$ the fine structure constant. The Dirac equation is given by:
$$ \left(\begin{matrix} i\frac{mc}{\hbar}+\frac{1}{c}\left(\frac{\partial}{\partial t}-i\frac{e}{\hbar}\phi\right) & \sigma_k\left(\frac{\partial}{\partial x^k}+i\frac{e}{\hbar c}A_k\right) \\ -\sigma_k\left(\frac{\partial}{\partial x^k}+i\frac{e}{\hbar c}A_k\right) & i\frac{mc}{\hbar}-\frac{1}{c}\left(\frac{\partial}{\partial t}-i\frac{e}{\hbar}\phi\right) \end{matrix}\right)\left(\begin{matrix} \psi^+\\ \psi^- \end{matrix}\right)=0 $$
in the limit where $\frac{mc}{\hbar}\gg \frac{e\phi}{\hbar c}$ it is convenient to take $\chi^{\pm}\equiv e^{i\frac{mc}{\hbar}t}\psi^{\pm}$, i.e. shifting the zero energy. Then we get:
$$ \left(\begin{matrix}\frac{1}{c}\left(\frac{\partial}{\partial t}-i\frac{e}{\hbar}\phi\right) & \sigma_k\left(\frac{\partial}{\partial x^k}+i\frac{e}{\hbar c}A_k\right) \\ -\sigma_k\left(\frac{\partial}{\partial x^k}+i\frac{e}{\hbar c}A_k\right) & 2i\frac{mc}{\hbar}-\frac{1}{c}\left(\frac{\partial}{\partial t}-i\frac{e}{\hbar}\phi\right) \end{matrix}\right)\left(\begin{matrix} \psi^+\\ \psi^- \end{matrix}\right)=0 $$
so
$$ \chi^{-}\sim-\frac{i\hbar}{2mc}\sigma^k\left(\frac{\partial}{\partial x^k}+i\frac{e}{\hbar c}A_k\right) \chi^{+}\ll\chi^{+} $$
since the dependence of $\chi^{+}$ over $x^k$ os dominated by the length $\frac{e\phi}{\hbar c}$. This means that at first order we can drop out the $\chi^{-}$ spinor and just work with the $\chi^+$. Alittle of algebra will give you
$$ \left[\frac{1}{2m}\left(\frac{\hbar}{i}\nabla-\frac{e}{c}\vec{A}\right)^2-\frac{e\hbar}{2mc}\vec{\sigma}\cdot\vec{B}+e\phi\right]\chi^{+}=i\hbar\frac{\partial}{\partial t}\chi^{+} $$
note that this formula predicts the gyromagnetic ratio $g=2$.
Now, for obtain the next order we cannot simply drop the $\chi^{-}$ out of the equation. For the first order expansion we use the fact that the spinors $\psi^{+}$ and $\psi^{-}$ are coupled and equally important but we can decouple them in this limit by the unitary transformation $\chi^{\pm}\equiv e^{i\frac{mc}{\hbar}t}\psi^{\pm}$. Now to obtain the next order we are going to generalize this unitary transformation to $\chi^{\pm}\equiv e^{iS}\psi^{\pm}$ such that the spinors $\chi^{+}$ and $\chi^-$ decouple up to terms of the next-next-order, i.e. $\sim \alpha ^{k+1}$ for order $k$ in $\alpha$ expansion. This unitary transformation are called Foldy–Wouthuysen transformation.
We can always organize the Dirac equation to looks like a Schrodinger equation, then under the Foldy–Wouthuysen transformation:
$$ H\left(\begin{matrix} \psi^+ \\ \psi^- \end{matrix}\right)=i\hbar\frac{\partial}{\partial t}\left(\begin{matrix} \psi^+ \\ \psi^- \end{matrix}\right)\rightarrow (e^{iS}He^{iS})\left(\begin{matrix} \chi^+ \\ \chi^- \end{matrix}\right)=i\hbar(e^{iS}\frac{\partial}{\partial t}e^{iS})\left(\begin{matrix} \chi^+ \\ \chi^- \end{matrix}\right) $$
Now, $H=\beta mc^2+\mathcal{T}+\mathcal{E}$, where just $\mathcal{T}$ is the responsible for coupling the spinors. The calculation is bit clumsy so I will skip that. The ideia is to find an $S$ for each order in $\alpha$, for example: keeping just $\alpha^0$ order we get
$$ (e^{iS}He^{iS})=\beta mc^2 +\mathcal{T}+\mathcal{E}+i[S,\beta]mc^2+\mathcal{O}(\alpha mc^2) $$
and then $S=-\frac{\beta\mathcal{T}}{2mc^2}$. Then we look at the next order terms $H'=\beta mc^2+\mathcal{T}'+\mathcal{E}'$ do a new Foldy–Wouthuysen transformation, obtaining $H''=\beta mc^2+\mathcal{T}''+\mathcal{E}''$.
The Hamiltonian for the $\alpha^4mc^2$ order is given by:
$$ H'''=mc^2+\frac{p^2}{2m}+e\phi-\frac{p^4}{8m^2c^3}+\frac{e}{2m^2c^2}\frac{1}{r}\frac{d\phi}{dr}L\cdot S + \frac{e\hbar^2}{8m^2c^2}\nabla^2\phi + \mathcal{O}(\alpha^5mc^2) $$
Turn out that the next order obtained from the Dirac action will not agree with the experiment since the next order get into the scale of QED effects ($\Delta E\sim\alpha^5 mc^2$) like the Lamb Shift.
You can see all this and more here
