Let $n:=d$. By rescalling, without loss of generality $\sigma=1$.
Let $X_1,\dots,X_n$ be iid random variables (r.v.'s) with the standard Laplace distribution, so that the joint pdf of $X:=(X_1,\dots,X_n)$ is
\begin{equation}
f_X(x)=\frac1{2^n}\,\exp\Big\{-\sum_1^n|x_i|\Big\}
\end{equation}
for $x=(x_1,\dots,x_n)$. Let $Y_1:=X_1+\dots+X_n$, $Y_2:=X_2,\dots,Y_n:=X_n$, so that $X_1=Y_1-Y_2-\dots-Y_n$, $X_2=Y_2,\dots,X_n=Y_n$. So, $X=(X_1,\dots,X_n)$ and $Y:=(Y_1,\dots,Y_n)$ are related by an invertible linear transformation of determinant $1$. Therefore, using transformation technique for systems of r.v.'s/change of variables in multi-fold integrals, we see that the joint pdf of $Y=(Y_1,\dots,Y_n)$ is
\begin{multline}
f_Y(y)=f_X(y_1-y_2-\dots-y_n, y_2,\dots,y_n) \\
=\frac1{2^n}\,\exp\Big\{-|y_1-y_2-\dots-y_n|-\sum_2^n|y_i|\Big\}
\end{multline}
for $y=(y_1,\dots,y_n)$. Therefore, the joint distribution of $X=(X_1,\dots,X_n)$ given that $X_1+\dots+X_n=0$ is determined by the conditional joint pdf of $(Y_2,\dots,Y_n)$ given $Y_1=0$, which is given by the expression
\begin{equation}
f_{Y_2,\dots,Y_n|Y_1=0}(y_2,\dots,y_n)=\frac{f_Y(0,y_2,\dots,y_n)}{f_{Y_1}(0)},
\end{equation}
where $f_{Y_1}$
is the pdf of $Y_1$.

So, to complete the calculation of $f_{Y_2,\dots,Y_n|Y_1=0}(y_2,\dots,y_n)$, it remains to compute $f_{Y_1}(0)$. This can be done quickly by using characteristic functions (c.f.'s). Indeed, the c.f. of $X_1$ is given by
\begin{equation}
\phi(t)=\int_{-\infty}^\infty e^{itx}\frac1{2}\,e^{-|x|}\,dx=\frac1{1+t^2}
\end{equation}
for real $t$. Hence, the c.f. of $Y_1=X_1+\dots+X_n$ is $t\mapsto\frac1{(1+t^2)^n}$ and
\begin{equation}
f_{Y_1}(0)=\frac1{2\pi}\int_{-\infty}^\infty e^{-it0}\frac{dt}{(1+t^2)^n}
=2^{-2 n-1} \binom{2 n}{n}.
\end{equation}

Thus, the conditional joint pdf of $(Y_2,\dots,Y_n)$ given $Y_1=0$ is
\begin{equation}
2^{n+1}\exp\Big\{-|y_2+\dots+y_n|-\sum_2^n|y_i|\Big\}\Big/\binom{2 n}{n}.
\end{equation}

To simulate this distribution, one can use Markov chain Monte Carlo methods.