Let $c \in {\mathbb R}$ and $d \in {\mathbb R}$. Without loss of generality we can consider a slightly different function:
\begin{equation}
f(c,d):= \int\limits_{-\infty}^\infty \exp(-t^2) \mbox{erfc}(t-c) \mbox{erfc}(t-d) dt
\end{equation}
Now we differentiate with respect to $c$. We have:
\begin{eqnarray}
\partial_c f(c,d) &=& \lim\limits_{\epsilon \rightarrow 0} \frac{\sqrt{2}}{\sqrt{\pi}} \exp(-\frac{c^2}{2})
\left( \sqrt{\pi} + 2 \sqrt{\pi} T(\epsilon,\frac{1}{\sqrt{2}},\frac{2d-c}{\sqrt{2}}) - 2 \sqrt{\pi} T(\epsilon,\frac{1}{\sqrt{2}},\frac{-2d+c}{\sqrt{2}})\right)\\
&=&\sqrt{2} \exp(-\frac{c^2}{2}) \left( 1-erf(\frac{c-2 d}{\sqrt{6}})\right)
\end{eqnarray}
where $T(h,a,b)$ is the generalized Owen's T function Generalized Owen's T function .
Now, all we need to do is to integrate. Since $f(-\infty,d)=0$ we integrate over $c$ from minus infinity to $c$. We have:
\begin{eqnarray}
f(c,d) &=& \sqrt{2} \left( \sqrt{\frac{\pi}{2}}(1+erf(\frac{c}{\sqrt{2}}) + 2erf(\frac{d}{\sqrt{2}})) + \sqrt{2 \pi} 2 T(c,\frac{1}{\sqrt{3}},-\frac{2 d}{\sqrt{3}})\right) \\
&=& \frac{1}{3 \sqrt{\pi}} \left( 12 \pi T\left(c,\frac{c-2 d}{\sqrt{3} c}\right)+12 \pi T\left(d,\frac{d-2 c}{\sqrt{3} d}\right)-6 \arctan\left(\frac{c-2 d}{\sqrt{3} c}\right)-6 \arctan\left(\frac{d-2 c}{\sqrt{3} d}\right)+3
\pi \text{erf}\left(\frac{c}{\sqrt{2}}\right)+3 \pi \text{erf}\left(\frac{d}{\sqrt{2}}\right)+4 \pi\right)
\end{eqnarray}
where $T(h,a)$ is Owen's T function https://en.wikipedia.org/wiki/Owen%27s_T_function .
In[1060]:= {c, d} = RandomReal[{-2, 2}, 2, WorkingPrecision -> 50];
NIntegrate[Exp[-t^2] Erfc[t - c] Erfc[t - d], {t, -Infinity, Infinity}]
(4 \[Pi] - 6 ArcTan[(c - 2 d)/(Sqrt[3] c)] -
6 ArcTan[(-2 c + d)/(Sqrt[3] d)] + 3 \[Pi] Erf[c/Sqrt[2]] +
3 \[Pi] Erf[d/Sqrt[2]] + 12 \[Pi] OwenT[c, (c - 2 d)/(Sqrt[3] c)] +
12 \[Pi] OwenT[d, (-2 c + d)/(Sqrt[3] d)])/(3 Sqrt[\[Pi]])
Out[1061]= 0.483318
Out[1062]= 0.48331807775609703646923225386370751256977344715
Update: Now let us take four real numbers $a_1 \in {\mathbb R}$, $a_2 \in {\mathbb R}$, $c\in {\mathbb R}$ and $d \in {\mathbb R}$ and consider a more general integral:
\begin{equation}
f^{(a_1,a_2)}(c,d):= \int\limits_{-\infty}^\infty \exp(-t^2) \mbox{erfc}(a_1 t-c) \mbox{erfc}(a_2 t-d) dt
\end{equation}
Then by doing the same calculations as above we easily arrive at the following formula:
\begin{eqnarray}
&&f^{(a_1,a_2)}(c,d)= \frac{1}{\sqrt{\pi}} \left(\right.\\
&&
4 \pi T\left(\frac{\sqrt{2} c}{\sqrt{a_1^2+1}},\frac{a_1 a_2 c-\left(a_1^2+1\right) d}{c \sqrt{a_1^2+a_2^2+1}}\right)+4 \pi T\left(\frac{\sqrt{2}
\sqrt{a_1^2+1} d}{\sqrt{a_1^2 a_2^2+a_1^2+a_2^2+1}},\frac{a_1 a_2 d-\left(a_2^2+1\right) c}{d \sqrt{a_1^2+a_2^2+1}}\right) +\\
&&-2 \arctan\left(\frac{a_1
a_2 c-\left(a_1^2+1\right) d}{c \sqrt{a_1^2+a_2^2+1}}\right)-2 \arctan\left(\frac{a_1 a_2 d-\left(a_2^2+1\right) c}{d \sqrt{a_1^2+a_2^2+1}}\right)+\\
&&\pi
\text{erf}\left(\frac{\sqrt{a_1^2+1} d}{\sqrt{a_1^2 a_2^2+a_1^2+a_2^2+1}}\right)+\pi
\text{erf}\left(\frac{c}{\sqrt{a_1^2+1}}\right)+\\
&&\pi +2 \arctan\left(\frac{a_1 a_2}{\sqrt{a_1^2+a_2^2+1}}\right)\\
&&\left.\right)
\end{eqnarray}
For[count = 1, count <= 200, count++,
{a1, a2, c, d} = RandomReal[{-5, 5}, 4, WorkingPrecision -> 50];
I1 = NIntegrate[
Exp[-t^2] Erfc[a1 t - c] Erfc[a2 t - d], {t, -Infinity, Infinity}];
I2 = 1/Sqrt[\[Pi]] (\[Pi] +
2 ArcTan[(a1 a2)/Sqrt[(1 + a1^2 + a2^2)]] -
2 ArcTan[ (a1 a2 c - (1 + a1^2) d)/(
Sqrt[(1 + a1^2 + a2^2)] c)] -
2 ArcTan[ (a1 a2 d - (1 + a2^2) c)/(
Sqrt[(1 + a1^2 + a2^2)] d)] + \[Pi] Erf[c/Sqrt[
1 + a1^2]] + \[Pi] Erf[(Sqrt[(1 + a1^2)] d)/ Sqrt[
1 + a1^2 + a2^2 + a1^2 a2^2]] +
4 \[Pi] OwenT[(Sqrt[2] c)/Sqrt[
1 + a1^2], (a1 a2 c - (1 + a1^2) d)/(
Sqrt[(1 + a1^2 + a2^2)] c)] +
4 \[Pi] OwenT[(Sqrt[2] Sqrt[(1 + a1^2)] d)/ Sqrt[
1 + a1^2 + a2^2 + a1^2 a2^2], (a1 a2 d - (1 + a2^2) c)/(
Sqrt[(1 + a1^2 + a2^2)] d)]);
If[Abs[I2/I1 - 1] > 10^(-3),
Print["results do not match", {a1, a2, c, d, {I1, I2}}]; Break[]];
If[Mod[count, 10] == 0, PrintTemporary[count]];
];
Your integral admits a closed form.
Theorem. Let $a$ and $b$ be any real numbers. Then
$$\int_0^{\infty }\!\! e^{\large-\left(\frac{a-\log x}{b}\right)^2}\!\! \left(\text{erf}\left(\frac{a-\log x}{b}\right)+1\right) \mathrm{d}x=\sqrt{\pi}\:|b|\:e^{{\large {a+\frac{b^2}{4}}}}\left(1-\text{erf}\!\left(\! \frac{\sqrt{2}\:b}{4}\!\right)\! \right) \tag1$$
where $\displaystyle \text{erf}(z)=\frac{2}{\sqrt{\pi }}\int _0^ze^{\large-t^2}\mathrm{d}t.$
Proof. Let's denote the left hand side of $(1)$ by $I(a,b)$.
By the change of variable $$\displaystyle U:=\frac{a-\log x}{b}, \quad dx=-b\:e^a \:e^{-b U}dU,$$ we have
$$
\begin{align}
I(a,b)&=|b|\:e^a \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \left(\text{erf}\left(U\right)+1\right) \mathrm{d}U\\\\
&=|b|\:e^a \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \text{erf}(U)\: \mathrm{d}U+|b|\:e^a \int_{-\infty }^{+\infty } e^{\large-U^2-bU} \mathrm{d}U \tag2
\end{align}
$$ Clearly, by the gaussian integral:
$$
\int_{-\infty }^{+\infty } e^{\large-U^2-bU} \mathrm{d}U =e^{\large b^2/4}\int_{-\infty }^{+\infty } e^{\large-(U+b/2)^2} \mathrm{d}U =\sqrt{\pi}\:e^{ b^2/4} .\tag3
$$ Set
$$
f(b):=\int_{-\infty }^{+\infty } e^{\large-U^2-bU} \text{erf}(U)\: \mathrm{d}U \tag4
$$ and observe that $$ f(0)=\frac{\sqrt{\pi}}{4}\left[(\text{erf}(U))^2\right]_{-\infty }^{+\infty }=0. \tag5$$ Differentiating $(4)$ with respect to $b$ and performing an integration by parts gives
$$
\begin{align}
f'(b)&=-\int_{-\infty }^{+\infty } U\:e^{\large-U^2-bU} \text{erf}(U)\:\mathrm{d}U\\\\
f'(b)&=\left[\frac 12 e^{\large-U^2}\left(e^{\large-bU} \text{erf}(U)\right)\right]_{-\infty }^{+\infty }-\frac 12\int_{-\infty }^{+\infty } e^{\large-U^2} \left(-b\:e^{\large-bU} \text{erf}(U)+\frac{2}{\sqrt{\pi }}e^{\large-bU}e^{\large-U^2}\right)\mathrm{d}U\\\\
f'(b)&=\frac{b}{2}f(b)-\frac{1}{\sqrt{\pi }}\int_{-\infty }^{+\infty } e^{\large-2U^2-bU}\mathrm{d}U
\end{align}
$$ or
$$
f'(b)=\frac{b}{2}f(b)-\frac{\sqrt{2}}{2}\:e^{\large b^2/8}. \tag6
$$ We classically solve the ODE $(6)$, using $(5)$, to obtain
$$
f(b)=-\sqrt{\pi }\:e^{\large b^2/4}\text{erf}\!\left(\! \frac{\sqrt{2}\:b}{4}\!\right) \tag7
$$ then plugging $(7)$, $(4)$ and $(3)$ into $(2)$ gives $(1)$ as desired.
Best Answer
Let $\;I(b)=\displaystyle\int_{-\infty}^{+\infty}e^{-ax^2}\mathrm{erf}\!\left(\!\dfrac{x+b}{\sqrt2}\!\right)\!\mathrm dx\,.$
By differentiating with respect to $\,b\,,$ we get that
$I’(b)=\displaystyle\int_{-\infty }^{+\infty}e^{-ax^2}\frac{\partial}{\partial b}\left[\mathrm{erf}\!\left(\!\dfrac{x+b}{\sqrt2}\!\right)\right]\!\mathrm dx=$
$=\dfrac2{\sqrt{\pi}}\displaystyle\int_{-\infty }^{+\infty}e^{-ax^2}\frac{\partial}{\partial b}\left(\!\int_0^{\frac{x+b}{\sqrt2}}e^{-t^2}\mathrm dt\!\right)\mathrm dx=$
$=\dfrac2{\sqrt{\pi}}\cdot\dfrac1{\sqrt2}\displaystyle\int_{-\infty }^{+\infty}e^{-ax^2}e^{-\frac{(x+b)^2}2}\mathrm dx=$
$=\dfrac{\sqrt2}{\sqrt{\pi}}\displaystyle\int_{-\infty }^{+\infty}e^{-ax^2-\frac12x^2-bx-\frac{b^2}2}\,\mathrm dx=$
$=\dfrac{\sqrt2}{\sqrt{\pi}}\displaystyle\int_{-\infty }^{+\infty}e^{-\frac{2a+1}2x^2-bx-\frac{b^2}2}\,\mathrm dx=$
$=\dfrac{\sqrt2}{\sqrt{\pi}}\displaystyle\int_{-\infty }^{+\infty}e^{-\frac{2a+1}2\left[x^2+\frac{2bx}{2a+1}+\frac{b^2}{(2a+1)^2}\right]+\frac{b^2}{2(2a+1)}-\frac{b^2}2}\,\mathrm dx=$
$=\dfrac{\sqrt2}{\sqrt{\pi}}\displaystyle\int_{-\infty }^{+\infty}e^{-\frac{2a+1}2\left(x+\frac b{2a+1}\right)^2-\frac{ab^2}{2a+1}}\,\mathrm dx=$
$=\dfrac{\sqrt2}{\sqrt{\pi}}\,e^{-\frac{ab^2}{2a+1}}\!\displaystyle\int_{-\infty }^{+\infty}e^{-\frac{2a+1}2\left(x+\frac b{2a+1}\right)^2}\,\mathrm dx=$
$\underset{\overbrace{\text{by letting }\,u=\sqrt{\frac{2a+1}2}\left(x+\frac b{2a+1}\right)}}{=}\dfrac{\sqrt2}{\sqrt{\pi}}\sqrt{\dfrac2{2a+1}}\,e^{-\frac{ab^2}{2a+1}}\!\displaystyle\int_{-\infty }^{+\infty}e^{-u^2}\,\mathrm du=$
$=\dfrac{\sqrt2}{\sqrt{\pi}}\sqrt{\dfrac2{2a+1}}\,e^{-\frac{ab^2}{2a+1}}\sqrt{\pi}=$
$=\dfrac2{\sqrt{\pi}}\sqrt{\dfrac{\pi}a}\sqrt{\dfrac a{2a+1}}\,e^{-\frac{ab^2}{2a+1}}=$
$=\displaystyle\dfrac2{\sqrt{\pi}}\sqrt{\dfrac{\pi}a}\!\cdot\!\dfrac{\partial}{\partial b}\left(\int_0^{b\sqrt{\frac{a}{2a+1}}}e^{-t^2}\mathrm dt\right)=$
$=\dfrac{\partial}{\partial b}\left[\sqrt{\dfrac{\pi}a}\,\mathrm{erf}\!\left(\!b\sqrt{\dfrac{a}{2a+1}}\right)\right].$
Consequently,
$I(b)=\sqrt{\dfrac{\pi}a}\,\mathrm{erf}\!\left(\!b\sqrt{\dfrac{a}{2a+1}}\right)+\mathrm{constant}\,.$
Since $\;I(0)=\displaystyle\int_{-\infty}^{+\infty}\underbrace{e^{-ax^2}\mathrm{erf}\!\left(\!\dfrac x{\sqrt2}\!\right)}_{\text{it is an odd function}}\,\mathrm dx=0= \sqrt{\dfrac{\pi}a}\,\mathrm{erf}\!\left(0\right)\,,\,$
it follows that $\;\mathrm{constant}=0\,,\,$ hence ,
$I(b)=\sqrt{\dfrac{\pi}a}\,\mathrm{erf}\!\left(\!b\sqrt{\dfrac{a}{2a+1}}\right)$
that is ,
$\displaystyle\int_{-\infty}^{+\infty}e^{-ax^2}\mathrm{erf}\!\left(\!\dfrac{x+b}{\sqrt2}\!\right)\!\mathrm dx=\sqrt{\dfrac{\pi}a}\,\mathrm{erf}\!\left(\!b\sqrt{\dfrac{a}{2a+1}}\right).$