In this Wikipedia article the polarization $$ \vec P = \rho_\mathrm{b} \vec d $$ is determined from the density of bound charges $\rho_\mathrm{b}$ and the displacement vector $\vec d$, which denotes the displacement of positive and negative charges in a small volume. If I recall correctly, the volume as a whole has neutral charge.
Now if I am not mistaken, the relation $\rho_\mathrm{b} = - \nabla \vec P$ should hold as well. But how can one verify this?
Imagine some volume filled with point-like charges at the location $\mathbf{R}$, the charge density is then $$\rho = \sum_{n=1}^N q_n \delta(\mathbf{r} - \mathbf{r}_n - \mathbf{R}) $$ where $\mathbf{r}_n$ is the vector that shows from the centre of the volume to the nth charge. The macroscopic charge density however is an average over this microscopic charge density, so $\rho_m = \langle \rho \rangle$. To average the microscopic charge density use some smearing function $g(\mathbf{r})$ and take the convolution with $\rho$ to get the averaged charge density at the location $\mathbf{r}$ \begin{align*} \langle \rho\rangle (\mathbf{r}) &= \int \mathrm{d}^3r^{\prime}\, g(\mathbf{r}^{\prime}-\mathbf{r}) \sum_{n=1}^{N}q_n \delta(\mathbf{r}^{\prime}-\mathbf{r}_n - \mathbf{R}) \\ &= \sum_{n=1}^{N} q_n \int \mathrm{d}^3r^{\prime}\, g(\mathbf{r}^{\prime}-\mathbf{r}) \delta(\mathbf{r}^{\prime}-\mathbf{r}_n -\mathbf{R}) \end{align*} If you now taylor expand the integrand, assuming that $\mathbf{r}<<\mathbf{R}$ you get the following terms \begin{align*} \langle \rho \rangle (\mathbf{r}) &= \sum_{n=1}^{N} q_n \int \mathrm{d}^3r^{\prime}\, g(\mathbf{r}^{\prime}-\mathbf{r}) \delta(\mathbf{r}^{\prime}-\mathbf{R}) \\ & + \sum_{n=1}^{N} q_n \int \mathrm{d}^3r^{\prime} g(\mathbf{r}^{\prime}-\mathbf{r}) \Big(-\mathbf{r}_n \cdot \nabla_{\mathbf{r}^{\prime}} \delta(\mathbf{r}^{\prime} - \mathbf{R}) \Big) \end{align*} Now the first term gives you a sum over the smeared out charges, which is simply 0 if the Volume has no net charge. The second term however can be further simplified with integration by parts \begin{align*} &\quad\ \sum_{n=1}^{N}q_n \mathbf{r}_n \cdot \int \mathrm{d}^3r^{\prime} \delta(\mathbf{r}^{\prime}-\mathbf{R}) \nabla_{\mathbf{r}^{\prime}} g(\mathbf{r}^{\prime}-\mathbf{r}) \\ &= \sum_{n=1}^{N}q_n \mathbf{r}_n \cdot \nabla_{\mathbf{R}}\,g(\mathbf{R}-\mathbf{r}) \\ &= \mathbf{d} \cdot \nabla_{\mathbf{R}}\, g(\mathbf{R}-\mathbf{r}) \end{align*} Where $\mathbf{d}$ is the microscopic dipole moment. We can now change the variables and get that the averaged bound charge density can be expressed as \begin{align*} \langle\rho\rangle(\mathbf{r}) &= -\mathbf{d} \cdot \nabla_{r}\, g(\mathbf{R}-\mathbf{r}) \\ &= - \nabla \cdot \mathbf{d}\, g(\mathbf{R}-\mathbf{r}) \end{align*} This smeared out dipole vector is the polarisation $\mathbf{P}$, so we get that $\langle \rho \rangle = - \nabla \cdot \mathbf{P}$. And recalling that $\mathbf{d} = \sum_{n=1}^{N} q_n \mathbf{r}_n$, so the polarisation can be written as \begin{align*} \mathbf{P}(\mathbf{r}) &= \sum_{n=1}^{N} q_n \mathbf{r}_n g(\mathbf{R}-\mathbf{r}) \end{align*} For simplicity let's look at two charges with opposite sign now so then \begin{align*} \mathbf{P}(\mathbf{r}) &= \Big( q \mathbf{r}_1 - q \mathbf{r}_2 \Big) g(\mathbf{R}-\mathbf{r}) \\ &= q \mathbf{D} g(\mathbf{R}-\mathbf{r}) \end{align*} where $\mathbf{D}$ is the displacement vector. Now for the units to work out $g$ has to have units of $1 / V$ so $q g$ is charge per Volume which is the charge density so $\mathbf{P} = \rho \mathbf{D}$
