The apparent contradiction is the following. The general non-relativistic Hamiltonian
$$
\begin{aligned}
{\cal{H}}&=
-\sum_{a=1}^{\cal{N}}\frac{1}{2m_a}\nabla_{\rho,a}^2+\sum_{\substack{a,b=1 \\ a<b}}^{\cal{N}}\frac{\alpha z_az_b}{|\vec{\rho}_a-\vec{\rho}_b|} \\
&=
-\sum_{i=1}^{N_{\text{e}}}\frac{1}{2m}\nabla_i^2+\sum_{\substack{i,j=1 \\ i<j}}^{N_{\text{e}}}\frac{\alpha}{|\vec{r}_i-\vec{r}_j|}
-\sum_{i}^{N_{\text{e}}}\sum_{A=1}^{N_{\text{nuc}}}\frac{\alpha Z_A}{|\vec{r}_i-\vec{R}_A|}
-\sum_{A=1}^{N_{\text{nuc}}}\frac{1}{2M_A}\nabla_A^2+\sum_{\substack{A,B=1 \\ A<B}}^{N_{\text{nuc}}}\frac{\alpha Z_AZ_B}{|\vec{R}_A-\vec{R}_B|}
\end{aligned}
$$
is invariant under inversions $\vec{\rho}_a\rightarrow-\vec{\rho}_a$ (that is, $(\vec{r}_i,\vec{R}_A)\rightarrow(-\vec{r}_i,-\vec{R}_A)$), meaning its non-degenerate energy eigenstates are also parity eigenstates which cannot carry a permanent dipole moment:
$$
\vec{\mu}=\langle\Psi|\sum_{A=1}^{N_{\text{nuc}}}Z_A\vec{R}_A-e\sum_{i=1}^{N_{\text{e}}}\vec{r}_i|\Psi\rangle=0 \ .
$$
However, at the same time, the clamped nucleus approximation of $\cal{H}$ (esentially the first half of the Born-Oppenheimer separation: letting $M_A\rightarrow\infty$ and treating $\vec{R}$ as parameters) leads to
$$
\begin{aligned}
H(\vec{R})=
-\sum_{i=1}^{N_{\text{e}}}\frac{1}{2m}\nabla_i^2+\sum_{\substack{i,j=1 \\ i<j}}^{N_{\text{e}}}\frac{\alpha}{|\vec{r}_i-\vec{r}_j|}
-\sum_{i=1}^{N_{\text{e}}}\sum_{A=1}^{N_{\text{nuc}}}\frac{\alpha Z_A}{|\vec{r}_i-\vec{R}_A|}
+\sum_{\substack{A,B=1 \\ A<B}}^{N_{\text{nuc}}}\frac{\alpha Z_AZ_B}{|\vec{R}_A-\vec{R}_B|} \ ,
\end{aligned}
$$
which does give rise to a dipole moment:
$$
\vec{\mu}=\sum_{A=1}^{N_{\text{nuc}}}Z_A\vec{R}_A-e\langle\Phi(\vec{R})|\sum_{i=1}^{N_{\text{e}}}\vec{r}_i|\Phi(\vec{R})\rangle=\sum_{A=1}^{N_{\text{nuc}}}Z_A\vec{R}_A-e
\int\mathrm{d}^3r\vec{r}\rho(\vec{r};\vec{R}) \ .
$$
In the above, $\vec{R}$ refers to the (parametric) dependence on all the nuclear coordinates, and $\rho$ is just the electronic density that can be calculated from $\Phi$. Neglecting any vibrational effects and simply using the equilibrium (minimal energy) nuclear configurations gives us permanent dipole moments that are in rather good agreement with experiment. For example[1][2]:
$$
\begin{aligned}
\mu_{\text{calc}}(\text{HCl})\approx 1.084 \, \text{D} \ \ \ &\leftrightarrow \mu_{\text{exp}}(\text{HCl})\approx 1.093 \, \text{D} \ , \\
\mu_{\text{calc}}(\text{H$_2$O})\approx 1.840 \, \text{D} \ \ \ &\leftrightarrow \mu_{\text{exp}}(\text{H$_2$O})\approx 1.857 \, \text{D} \ .
\end{aligned}
$$
How could we generate a non-zero dipole moment from the clamped nucleus approximation, and how can such a value be extracted from the general formalism? The punchline is that dipole moment is only zero in a laboratory-fixed frame due to free rotations/inversions. Clamping the nuclei, however, picks out a molecule-fixed frame (position and orientation being fixed by the point charges of the nuclei), in which a non-zero dipole moment can be found; this corresponds to the fact that a molecule of $N_\text{nuc}$ nuclei has only $3N_\text{nuc}-6$ internal degrees of freedom ($3N_\text{nuc}-5$ when linear). The question is whether such a molecule-fixed frame can be found without clamping the nuclei.
The first thing to realize is that the kinetic energy associated with center-of-mass motion must be excluded from ${\cal{H}}$, otherwise the spectrum would be continuous, and no internal motion could be described:
$$
{\cal{H}}=-\frac{1}{2M_{\text{tot}}}\nabla^2_{\text{COM}}+H_\text{in} \ .
$$
There are, however, an infinite number of coordinate transformations that can achieve this[3][4]! Any invertible linear transformation
$$
\vec{\rho}'_a=\sum_{b=1}^{\cal{N}}t_{ab}\vec{\rho}_b
$$
satisfying the additional conditions
$$
t_{1b}=\frac{m_b}{M_\text{tot}} \ \ \ , \ \ \
\sum_{b=1}^{\cal{N}}t_{ab}=\delta_{a1}
$$
is sufficient to completely decouple the motion of $\vec{\rho}'_1=\vec{\cal{R}}_{\text{COM}}$.
Let us find the appropriate molecule-fixed frame for diatomic molecules; the rest of my answer is limited to this case. The internal coordinates are specifically
$$
\begin{aligned}
\vec{\rho}'_1&=\vec{R}'_1=\vec{\cal{R}}_{\text{COM}} \ , \\
\vec{\rho}'_2&=\vec{R}'_2=\vec{R}_2-\vec{R}_1 \ , \\
\vec{\rho}'_{i+2}&=\vec{r}'_i=\vec{r}_i-\frac{1}{2}(\vec{R}_1+\vec{R}_2) \ ,
\end{aligned}
$$
so that an internuclear vector $\vec{R}=\vec{R}'_2-\vec{R}'_1$ emerges, and all electron coordinates are measured from the geometric center of the molecule (it is easy to check that this is a special case of the above general transformation).
This transformation leads to
$$
H_\text{in}=-\frac{1}{2M_+}\nabla_R^2-\frac{1}{2m}\sum_{i=1}^N{\nabla'}_i^2
-\frac{1}{8M_+}\sum_{i,j=1}^N{\nabla'}_i\cdot{\nabla'}_j+\frac{1}{2M_-}\nabla_R\cdot\sum_{i=1}^N{\nabla'}_i+V \ ,
$$
where we introduced
$$
\frac{1}{M_\pm}=\frac{1}{M_2}\pm\frac{1}{M_1} \ .
$$
Note that this frame is space-fixed in the sense that its orientation is still tied to that of the original laboratory-fixed frame. To finally have the molecule-fixed frame, we write
$$
\vec{R}=R\cos(\phi)\sin(\theta)\vec{e}'_x+R\sin(\phi)\sin(\theta)\vec{e}'_y+R\cos(\theta)\vec{e}'_z \ ,
$$
and switch to a new coordinate system whose "$z$" axis coincides with $\vec{R}$ in the primed coordinate system:
$$
\begin{aligned}
\vec{e}''_x&=\frac{\partial_\theta\vec{R}}{|\partial_\theta\vec{R}|}=\cos(\phi)\cos(\theta)\vec{e}'_x+\sin(\phi)\cos(\theta)\vec{e}'_y-\sin(\theta)\vec{e}'_z \ , \\
\vec{e}''_y&=\frac{\partial_\phi\vec{R}}{|\partial_\phi\vec{R}|}=-\sin(\phi)\vec{e}'_x+\cos(\phi)\vec{e}'_y \ , \\
\vec{e}''_z&=\frac{\vec{R}}{R}=\cos(\phi)\sin(\theta)\vec{e}'_x+\sin(\phi)\sin(\theta)\vec{e}'_y+\cos(\theta)\vec{e}'_z \ .
\end{aligned}
$$
The notation is somewhat confusing, since these are actually the $\vec{e}_\theta,\vec{e}_\phi,\vec{e}_R$ unit vectors of the spherical coordinate system, but this is the convention set by Kolos and Wolniewicz[5] that everyone seems to follow. Using
$$
\vec{r}'_i=x'_{i}\vec{e}'_x+y'_{i}\vec{e}'_y+z'_{i}\vec{e}'_z=
x''_{i}\vec{e}''_x+y''_{i}\vec{e}''_y+z''_{i}\vec{e}''_z
$$
the components in the space-fixed (primed) and molecule-fixed (double-primed) systems can be seen to related:
$$
\begin{bmatrix}
x_i'' \\
y_i'' \\
z_i''
\end{bmatrix}
=
\begin{bmatrix}
\cos(\phi)\cos(\theta) & \sin(\phi)\cos(\theta) & -\sin(\theta) \\
-\sin(\phi) & \cos(\phi) & 0 \\
\cos(\phi)\sin(\theta) & \sin(\phi)\sin(\theta) & \cos(\theta)
\end{bmatrix}
\begin{bmatrix}
x_i' \\
y_i' \\
z_i'
\end{bmatrix}
\ ,
$$
and with this in hand, the Hamiltonian can be transformed to the molecule-fixed coordinates. The dipole moment operator is defined in this very frame:
$$
\hat{\vec{\mu}}=\frac{Z_2-Z_1}{2}R\vec{e}''_z-e\sum_{i=1}^n\vec{r}''_i \ ,
$$
with the most important property that it does not flip sign upon inversion in the space-fixed coordinate system. The angles change as $\theta\rightarrow\pi-\theta$, $\phi\rightarrow\pi+\phi$ upon inversion, which means that[5][6]
$$
(x_i',y_i',z_i')\rightarrow(-x_i',-y_i',-z_i') \ \ \ \Leftrightarrow \ \ \ (x_i'',y_i'',z_i'')\rightarrow(-x_i'',+y_i'',+z_i'') \ .
$$
The dipole moment thus does not vanish by symmetry in the molecule-fixed frame and it can in principle be calculated as an expectation value with the eigenstates of $H_\text{in}$ (one can show that it does vanish for the homonuclear case $Z_1=Z_2$, $M_1=M_2$, as it should).
The main (or basically only) application of this formalism was the perturbative calculation of the dipole moment for the hydrogen deuteride (HD) molecule. The clamped nucleus approximation could only predict $\vec{\mu}(\text{HD})=\vec{\mu}(\text{H}_2)=\vec{\emptyset}$, since it is insensitive to isotopic mass differences; in reality, HD does have a tiny dipole moment of roughly ${\mu}_\text{exp}(\text{HD})\approx8.78\cdot10^{-4} \, \text{D}$. Blinder, Kolos, Wolniewicz and many other workers of the 60s-80s calculated this to be roughly ${\mu}_\text{calc}(\text{HD})\approx8\cdot10^{-4} \, \text{D}$ with the above approach, although the results did not (and to my knowledge, still do not) agree about the second significant digit.
See Ref. [4] about the details of this technique and for the references to the aforementioned original works. See also Ref. [7] for current experimental and theoretical values, and for an alternative post-Born-Oppenheimer calculation.
This answer is incomplete in the sense that I only dealt with diatomic molecules. Molecule-fixed frames can be found for more than two nuclei as well, but I do not know of any dipole moment calculation for such systems (I doubt anyone actually tried it, if nothing else, due to the severe computational cost). Also, for polyatomic molecules, questions like this one are closely related to the question of whether molecular structure is a meaningful concept beyond the Born-Oppenheimer framework, which is not settled, to say the least.
References
$ $
[1]: NIST database
$ $
[2]: Damour, Quintero-Monsebaiz, Caffarel, Jacquemin, Kossoski, Scemama, Loos: J. Chem. Theory Comput. 19 1 221 (2022)
$ $
[3]: Sutcliffe, Woolley: Phys. Chem. Chem. Phys. 7 3664 (2005)
$ $
[4]: Fernandez, Echave: Chapter 6. of Computational Spectroscopy (editor: Jorg Grunnenberg; 2010) (arXiv link)
$ $
[5]: Kolos, Wolniewicz: Rev. Mod. Phys. 35 3 473 1963
$ $
[6]: Landau, Lifshitz: Quantum Mechanics $-$ Non-relativistic Theory; Sec. 86.
$ $
[7]: Hobson, Valeev, Csaszar, Stanton: Mol. Phys. 107 8 1153 2009