The above answer may be generalized to cover the case where rho is a user selected value and the values $E(x)=E(y)=\mu$ and $V(x)=V(y)=\sigma^2$ are equal and known. In this case:
$E(x)=a+b=E(y)=a+c=\mu$, therefore $b=c$ and $(a+b)*(a+c)=\mu^2$ (note $\mu<1$)
The $a,b,c$,and $d$ values are joint probability values such that:
$a+b+c+d=1$ by definition. So that $d= 1-(a+b+c)$
And by using the definition of the variance for Bernoulli random variables, by setting $V(x)=V(y)=\sigma^2$ the equation below follows:
$V(x)=(a+b)\cdot[1-(a+b)]=V(y)=(a+c)*[1-(a+c)]=\sigma^2=(\mu^2)\cdot(1-\mu)^2$
Then, using the originally derived equations for rho(xy)and Cov(x,y) above, it is possible to define rho in terms of the values of just a and mu, such that:
$\rho= \frac{a - \mu^2}{\mu(1-\mu)}$
Then solving for a we see that all values are known from rho, mu and sigma:
$$a= \mu^2 + [\rho(\mu(1-\mu))]\\
b= \mu - a\\
c= b\\
d= 1-(a+b+c)$$
Two correlated Bernoulli random variables can be simulated by first defining a discrete random variable $Z$ with four values: $3,2,1,0$ with the respective probabilities of $a,b,c,d$ corresponding to the four joint probability states. Then,
$x=1$ when $Z=3$ or $2$, $x=0$ otherwise
$y=1$ when $Z=3$ or $1$, $y=0$ otherwise
A Monte Carlo simulation with an Excel add-in then demonstrates that the $x$ and $y$ variables are correlated as specified by rho given the known values of $\mu$ and $\sigma$. The Excel CORREL function for $5000$ samples in a Monte Carlo simulation verifies this. Also, as a cross check when $\rho=0$ you may quickly show that the values of $a,b,c$, and $d$ are properly computed for the case where $x$ and $y$ are independent random variables.
Thanks to all of the original contributors for the work that prompted the above derivation.