Here's a generalization of gowers's construction that is practical for Bernoulli RVs with $p\ne1/2$. You want to generate $n$ Bernoulli RVs, each taking on value 0 or 1, each with mean $p<1/2$. (For $p>1/2$, do as described below for $1-p$, then complement the results.) Let $d=\sqrt{2}\text{ erfc}^{-1}\left(2p\right)$. (This is just the inverse survival function for the standard normal distribution.) Take unit vectors $v_1,...v_n$ as before. Generate a random $n$-vector $z$ whose components are IID standard normal RVs. Let $B_i=1$ iff $z\cdot v_i>d$. $z\cdot v_i$ is standard normal, so obviously gives the desired mean of $p$.

What about correlations? As in gowers's construction, these depend uniquely on the angle between vectors $v_i$ and $v_j$. Let $c_{ij}$ be the coincidence frequency between $B_i$ and $B_j$, i.e., the frequency with which both are 1, which is related to the correlation. If $\theta_{ij}=\cos^{-1}\left(v_i\cdot v_j\right)$, then

$$c_{ij}=\int_d^\infty \Phi\left(\frac{u\cos\theta_{ij}-d}{\sin\theta_{ij}}\right)\phi\left(u\right)du$$

where $\Phi(z)$ and $\phi(z)$ are the standard normal CDF and PDF, respectively. $c$ decreases monotonically from $p$ at $\theta=0$ to 0 at $\theta=\pi$. In a practical problem you'd probably want the inverse: you'd know $p$ and $c$ and want to get $\theta$. I doubt that can be done other than numerically, but $c$ is a single function of two bounded variables $p$ and $\theta$, so you can tabulate it numerically once and invert the interpolated function if you're going to be doing a lot of this.

Now you know what all dot products $v_i\cdot v_j=\cos{\theta_{ij}}$ need to be, it is simple to construct vectors at these angles. Let $v_1=\left(1,0,...,0\right)$. Then $v_2=\left(\cos\theta_{12},\sin\theta_{12},0,...,0\right)$. For $v_3$, solve

$$\pmatrix{v_{11}&v_{12}\cr v_{21}&v_{22}}\pmatrix{v_{31}\cr v_{32}}=\pmatrix{1&0\cr \cos\theta_{12}&\sin\theta_{12}}\pmatrix{v_{31}\cr v_{32}}=\pmatrix{\cos\theta_{13}\cr\cos\theta_{23}}$$

... then let $v_{33}=\sqrt{1-v_{31}^2-v_{32}^2}$. Continue to generate the rest of the $v_i$. Since the matrix at every stage is lower triangular, the solution is unique as long as the diagonal is positive. The construction fails only if the norm of the first $i-1$ components of $v_i$ is $\ge1$. I'm going to speculate that that occurs only if you give it a set of impossible coincidence frequencies (for instance, $c_{12}=c_{13}=p$, $c_{23}=0$), but I haven't attempted to show that.

Edit: Nope, I was too optimistic. For instance, if you have three mutually exclusive Bernoulli RVs with $p\le1/3$, which is clearly possible, this construction fails. Alas.