I am trying to understand this answer.
The answer is the following:
I'll assume that $V'$ is bounded and that $\mathbf{M'}$ is continuously differentiable. To your primary question:
Is [the] divergence theorem applicable to of equation (4)?
The answer is, strictly speaking, no. The hypotheses of the divergence theorem include that the vector field in question is continuously differentiable, and $\frac{\mathbf{M'(r)}}{r}$ is not (in general). However, one can still rigorously employ the divergence theorem to yield your expression by removing a ball $B_\epsilon(P)$ from $V'$. Writing $B_\epsilon '(P) \equiv B_\epsilon(P) \cap V' $,
$$\iiint_{V'} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' = \iiint_{V' \backslash B_\epsilon'(P)} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' + \iiint_{B_\epsilon'(P)} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' $$ The divergence Theorem may be applied to the first term for any $\epsilon>0$, and one can show that the second term goes to $0$ in the limit $\epsilon \to 0$ via your chain of equalities in equation $(4)$, which hold $a.e.$ in $B_\epsilon'(P)$, since $\mathbf{M'}$, being continuous on the compact set $\overline{V'}$, is bounded, say by $M>0$.
With some computation, one can show that the integral over $\partial(V' \backslash B_\epsilon '(P))$ arising from the first term approaches the integral over $S' = \partial V'$ in the same limit, assuming $S'$ is smooth in a neighborhood of $P$. To be explicit,
$$ \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' = \left( \iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)} - \iint_{\partial B_\epsilon '(P) \backslash \partial V'} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' $$
That is, the integrals are the same up to the difference of the integrals over the piece of $\partial B_\epsilon '(P)$ from $\partial V'$ and that from $\partial B_\epsilon(P)$. Since $S'$ is smooth near $P$, it is locally the graph of a smooth function $f(x,y)$ on the tangent plane at $P$ with standard linear coordinates $(x,y)$. So, for $\epsilon$ sufficiently small, we may compute the first integral on the RHS above as an integral in the tangent plane, with area form $\sqrt{1 + f_x^2 + f_y^2}dxdy=\sqrt{1 + f_x^2 + f_y^2}sdsd\theta$ in polar coordinates $(s,\theta)$. The term $\sqrt{1 + f_x^2 + f_y^2}$ is continuous and therefore bounded, say by $C$ (independent of $\epsilon$ for $\epsilon$ sufficiently small), and further $r'(s,\theta)=\sqrt{s^2+f(s,\theta)^2} \geq s$. Hence
$$\left|\iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)}\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'\right| \leq C \int_0^\epsilon \int_0^{2\pi}|\mathbf{M}(s,\theta)| \frac{s}{\sqrt{s^2+f^2}}dsd\theta \leq 2\pi MC \epsilon$$
Note that the exact integration here may not range over all of $0 \leq s \leq \epsilon$, but the bound still holds. Finally, in the second term on the RHS above, $r'=\epsilon$, so that integral is bounded by $4\pi M\epsilon$, and we have
$$\left| \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' \right| \leq 2\pi M(2+C)\epsilon $$
which goes to $0$ as $\epsilon \to 0$. Since my first equation holds for all $\epsilon >0$, we may take the limit of both sides as $\epsilon \to 0$ to obtain
$$\iiint_{V'} \left[\nabla' \cdot \left( \frac{\mathbf{M(r')}}{r'} \right) \right]dV' = \unicode{x222F}_{\partial V'} \left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'$$ As a naive application of the divergence theorem would have suggested.
Note that that "singularity" in the RHS integration does not cause the integral to diverge for essentially the same reason that it didn't in the volume integral: by passing to an appropriately defined set of coordinates in a neighborhood of $P$, we see that the area form gets small fast enough in this neighborhood for the integral to converge, just as the volume form did in your computations.
My misunderstandings are the following:
(1) What does $$ \left( \unicode{x222F}_{\partial V'} - \unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' = \left( \iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)} - \iint_{\partial B_\epsilon '(P) \backslash \partial V'} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS' $$ mean? Why are they equal?
My interpretation:
$\unicode{x222F}_{\partial V'}$ means the integral over the whole of the surface of $V'$.
$\text{ }$
$\unicode{x222F}_{\partial(V' \backslash B_\epsilon '(P))}$ means the integral over the surface of "$B_\epsilon '(P)$ removed from $V'$.
I have no idea about $\displaystyle\iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)}$ and $\displaystyle\iint_{\partial B_\epsilon '(P) \backslash \partial V'}$ Please explain.
Original answerer here. You are correct regarding what the two integrations mean individually. The notation on the LHS of the equation under consideration simply means the difference of the surface integrals over the two domains you've identified (of the same integrand $\frac{\mathbf{M}'(\mathbf{r}') \cdot \mathbf{\hat{n}}}{r'}$). Clearly the result will be an integral of this same integrand; the question is the domain.
The integrations over the intersection of the two original domains cancel, so our resulting integration won't contain any of the boundary of your rectangular prism except the bit above the removed ball, and the only other surface we need to include is the portion of the boundary of $B_\epsilon '(P)$ inside the rectangular prism. Note that $B_\epsilon'(P)$ in your drawn scenario is a half-ball, the boundary of which is a hemisphere (part of $\partial B_\epsilon(P)$, which is the full sphere of radius $\epsilon$ centered at $P$) plus a circular disc (part of $\partial V'$), and these two pieces are exactly what we've identified as what we want to integrate over. Note each of these pieces came from a different term on our LHS, so they'll be integrated over on the RHS with different signs. In particular, the disc should be integrated over with a positive sign, while the hemisphere should be integrated over with a negative sign. That is, our integral should be $$ \left( \iint_{\text{disc}} - \iint_{\text{hemisphere}} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'$$
All that remains is to check that my original notation agrees with this. But $\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)$ just means "the hemisphere and disc, remove anything contained in the full sphere", which is just the disc, while $\partial B_\epsilon '(P) \backslash \partial V'$ just means "the hemisphere and disc, remove anything contained in the rectangular prism's boundary", which is just the hemisphere, so the above integral is the same as $$ \left( \iint_{\partial B_\epsilon '(P) \backslash \partial B_\epsilon(P)} - \iint_{\partial B_\epsilon '(P) \backslash \partial V'} \right)\left[ \frac{\mathbf{M'(r')} \cdot \hat{\mathbf{n}}}{r'} \right]dS'$$ As our equation claims. The advantage of using this notation is that it represents the correct domain no matter what the shape of $V'$ is, and it emphasizes that the difference only depends on a region within distance $\epsilon$ of $P$.