I am trying to perform the following integral.
$$I=\int d^{3} \vec{r_{1}} \int d^{3} \vec{r_{2}} \exp\Bigg(-\mid\vec{r_{1}}- \vec{R_{1}}\mid^{2}-\mid\vec{r_{1}}- \vec{R_{2}}\mid^{2}-\mid\vec{r_{1}}- \vec{R_{2}}\mid^{2}-\mid\vec{r_{2}} - \vec{R_{1}}\mid^{2}\Bigg)$$
Where, $$\vec{r_{1}}=x_{1} \vec{i}+y_{1} \vec{j}+z_{1} \vec{k}$$ and $$\vec{r_{2}}=x_{2} \vec{i}+y_{2} \vec{j}+z_{2} \vec{k}$$
and $\vec{R_{1}}$and $\vec {R_{2}}$ are constant vectors. $d^{3} \vec{r_{1}}$ and $d^{3} \vec{r_{2}}$ are volume integrals integrations is over full volume. (i.e) $$-\infty<x_{1}<\infty$$ $$-\infty<y_{1}<\infty$$ $$-\infty<z_{1}<\infty$$$$-\infty<x_{2}<\infty$$ $$-\infty<y_{2}<\infty$$ $$-\infty<z_{2}<\infty$$
I'm assuming, as @bob.sacamento mentioned, that your third term is supposed to have $\mathbf{r}_2$ in it instead of $\mathbf{r}_1$. Here are some simplifications: \begin{align*} &\int_{\mathbb{R}^3}\int_{\mathbb{R}^3}\exp\left(-|\mathbf{r}_1-\mathbf{R}_1|^2-|\mathbf{r}_1-\mathbf{R}_2|^2-|\mathbf{r}_2-\mathbf{R}_1|^2-|\mathbf{r}_2-\mathbf{R}_2|^2\right)\,d^3\mathbf{r}_1 \, d^3\mathbf{r}_2 \\ =&\int_{\mathbb{R}^3}\exp\left(-|\mathbf{r}_1-\mathbf{R}_1|^2-|\mathbf{r}_1-\mathbf{R}_2|^2\right)\,d^3\mathbf{r}_1 \cdot \int_{\mathbb{R}^3}\exp\left(-|\mathbf{r}_2-\mathbf{R}_1|^2-|\mathbf{r}_2-\mathbf{R}_2|^2\right)\,d^3\mathbf{r}_2 \\ =&\left[\int_{\mathbb{R}^3}\exp\left(-|\mathbf{r}_1-\mathbf{R}_1|^2-|\mathbf{r}_1-\mathbf{R}_2|^2\right)\,d^3\mathbf{r}_1\right]^2. \end{align*} Let $\mathbf{r}=\mathbf{r}_1-\mathbf{R}_1$ and $\mathbf{R}=\mathbf{R}_1-\mathbf{R}_2.$ Then $d^3\mathbf{r}_1=d^3\mathbf{r}$, and the integral becomes $$\left[\int_{\mathbb{R}^3}\exp\left(-|\mathbf{r}|^2-|\mathbf{r}+\mathbf{R}|^2\right)\,d^3\mathbf{r}\right]^2. $$ Now we suppose that $\mathbf{r}=x\mathbf{i}+y\mathbf{j}+z\mathbf{k},$ and that $\mathbf{R}=X\mathbf{i}+Y\mathbf{j}+Z\mathbf{k}.$ The integral becomes \begin{align*}&\left[\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\exp\left[-(x^2+y^2+z^2)-((x+X)^2+(y+Y)^2+(z+Z)^2)\right] \, dx \, dy \, dz\right]^{\,2}\\ =&\left[\int_{-\infty}^{\infty}e^{-x^2-(x+X)^2}\,dx\cdot\int_{-\infty}^{\infty}e^{-y^2-(y+Y)^2}\,dy\cdot\int_{-\infty}^{\infty}e^{-z^2-(z+Z)^2}\,dz\right]^{\,2} \end{align*} These integrals can be computed separately thus: $$\int_{-\infty}^{\infty}e^{-x^2-(x+X)^2} \, dx=\sqrt{\frac{\pi}{2}}\,e^{-X^2/2}. $$ You should be able to finish.