Generating a timeseries with an arbitrary power spectrum
You are reversing your phases okay.  You should be able to use an even value of N if you set the phase for N/2 to zero.

The problem you have is that your PSD spans the entire N.  Unless $ P_k = P_{N-k} $ you are not generating a DFT with complex conjugate symmetry which is what you need to produce a real valued signal.

Also, you should be using the "ifft" instead of the "fft" to do an inverse FFT.

Hope this helps,

Ced

Edit:  I didn't see your edit before I posted.  It seems you figured this out.

======================

Followup:

For even N:


def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2-1)*2*np.pi
    phase[1:N/2] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase