What is the multivariate distribution of $(X_1, \ldots, X_n | X_1+ \dotsm + X_n = y)$, given $y$ and that $X_i$ is (otherwise; without being given $y$) unconditionally independently normally distributed with mean $\mu $ and standard deviation $\sigma$?
And, to the task at hand: How do I simulate random drawings of $X$ given $y,~ \mu,$ and $\sigma$?
As pointed out by @whuber, we can simply use the known formulas for conditional distributions in the multinormal, as given for instance at Deriving the conditional distributions of a multivariate normal distribution (we will use notation from there). At first I doubted, since the full covariance matrix $\Sigma$ in this case is singular, but a close reading of the proof by user macro shows that we only need that $\Sigma_{22}$ is non-singular. That $\Sigma$ itself is singular does not matter.
For this case, we can easily compute (details not given) that $$ \Sigma =\begin{pmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{pmatrix} =\sigma^2 \begin{pmatrix}I_n & 1_n \\ 1_n^T & n \end{pmatrix} $$ where $I_n$ is the $n \times n$ identity matrix, and $1_n$ is the column vector with $n$ $1$'s. Then using the formulas we find that the conditional distribution of $X$ given that $Y=y$, where $Y=\sum_1^n X_i$, is the n-dimensional multinormal distribution with mean $$ \mu_{x|y} = 1_n \frac{y}n $$ that is, all the components have the same conditional expectation $y/n$, and covariance matrix $$ \Sigma_{x|y} = \sigma^2 \begin{pmatrix} I_n - 1_n 1_n^T / n \end{pmatrix} $$ Note that $1_n 1_n^T$ is an $n \times n$-matrix with all components 1.
