Well you want to go from QFT to Classical mechanics.
Let's do this in three steps
1. QED to Dirac Equation
QED lagrangian with electric dipole is
$\mathcal{L} = \bar{\psi}\left(\gamma\cdot\Pi - \frac{\mathrm{i}d}{2}\sigma^{\mu\nu}\gamma^5 F_{\mu\nu}\right)\psi\\$
Where $\Pi\equiv \partial - \mathrm{i}eA$. This implies the hamiltonian $$\mathcal{H} = \gamma^0 m + \gamma^0 \vec{\gamma}\cdot \vec{\Pi} + e\phi -\mathrm{i}d\left(-\vec{\gamma}\gamma^0\cdot \vec{B} + \mathrm{i}\vec{\gamma}\gamma^0\gamma^5\cdot\vec{E}\right)$$
Now using the Dirac representation of gamma matrices, the eigenvalue equation $\mathcal{H}\Psi = \mathcal{E}\Psi$ becomes $$\left( \begin{array}{cc} e\phi + m - \mathcal{E} + \mathrm{i}d\vec{\sigma}\cdot\vec{E} & \vec{\sigma}\cdot\vec{\Pi} + d \vec{\sigma}\cdot\vec{B} \\ -\vec{\sigma}\cdot\vec{\Pi} + d \vec{\sigma}\cdot\vec{B} & -e\phi + m + \mathcal{E} + \mathrm{i}d\vec{\sigma}\cdot\vec{E} \end{array} \right) \left( \begin{array}{c} \psi_+ \\ \psi_-\end{array} \right) = 0$$
in units of $c=1$.
2. Dirac to Schrodinger-Pauli
In the weak field non relativistic limit eq2 shows that $\psi_- \ll\psi_+$, so that solving for the latter we get an equation for the former
$$\left( \mathcal{E} - m\right)\psi_+ = \left(\frac{1}{2m}\left(\vec{\sigma}\cdot\vec{\Pi}\right)^2 - \frac{d}{2m}\left[ \vec{\sigma}\cdot\vec{\Pi}, \vec{\sigma}\cdot\vec{B}\right] + d \vec{E}\cdot\vec{\sigma} + e\phi\right)\psi_+$$
Now using the pauli matrix identities the commutator can be rewritten as $\frac{\hbar d}{m}\vec{\sigma}\cdot\left(\vec{\nabla}\times\vec{B}\right)+ -d^2B^2$, now we used the weak field limit (only first order in $\vec{E}$ and $\vec{B}$) so that we can also approximate $\vec{\Pi}\approx \vec{p} \rightarrow -i\hbar\vec{\nabla} $.
Using maxwell's equations for stationary fields we finally get that the hamiltonian in a stationary EM field is $$H = \frac{\left(\vec{p}- e\vec{A}\right)^2}{2m} + \frac{e\hbar}{2m}\vec{\sigma}\cdot\vec{B} + e\phi + d\vec{\sigma}\cdot \vec{E}$$
3. to Classical Electrodynamics
we deduce from he last term that the electron has a dipole moment of $\vec{d} \propto \vec{S}d$, where $\vec{S}$ stands for spin. Because as you know Electrostatic energy is$$ E= \int dx^3\phi\rho = \phi_0\int dx^3\rho + \nabla \phi_0\cdot\int dx^3 \vec{x}\rho + \ldots\rightarrow q\phi_0 + \vec{d}\cdot\vec{E}_0+\ldots$$
Note: no guarantee there are no signs or i's mistakes, but this is irrelevant to demonstrating the point
Quite simply, the dipole moment of a charged system depends on the coordinate origin. There is nothing particularly surprising or unphysical about this, and there are plenty of other quantities (such as orbital angular momentum) with that property. The dipole moment of a charged system is not a quantity that is defined for the system itself; it is only defined for the system in relation to a given coordinate system.
Best Answer
The dipole moment of a continuos charge distribution is dependent (unless your total charge is 0) by the pole you choose. If we set our pole equal to the origin of our system then we have: $$\vec{p}=\int_V \rho (\vec{r}')\vec{r}' d\tau$$ where $\vec{r}'$ is the vector position from your origin to the point you consider during integration(if you have a discrete charge distribution, then $\vec{r_i}'$ is the position of the charge $q_i$ with respect to the origin, $i$ is the index of the sum that replace the integral), and remember that you are integrating upon $d\tau=dxdydz$ which is expressed with respect to the origin. Now let's say we choose another point $A=(x_A,y_A,z_A)$. Reasoning with discrete distribution we say that we are searching the contribution of the $q_i$ charge that is still at the point $P'=\vec{r_i}'$ (we didn't change the origin) to the moment but with respect to another point $A$ which is in turn in the position $\vec{r}_A$ (thus we have the term $(\vec{r_i}'-\vec{r}_A)$), but the volume upon we are integrating is the same $d\tau$ and the point $P'$ is expressed by the same $\vec{r_i}'$ because we didn't change the origin, but only the pole. It follows that the moment with respect to the point $A$ is: $$\vec{p}_A=\int_V \rho(P'=\vec{r}')(\vec{r}'-\vec{r}_A)d\tau$$ where $V$ and $d\tau=dxdydz$ are the same of before.