SPL values from FFT of microphone signal
Let's assume that signal you are analysing is sinusoid with amplitude $A$:

$x=a\sin{2\pi f_{0} t}$

Its RMS value of amplitude is then: $\dfrac{a}{\sqrt{2}}$ as you noticed in your code. Before performing DFT it is good practice to window signal as it decreases leakage, etc. 

After transforming this signal into the frequency domain, you will obtain a two-sided spectrum. Then we have to take the magnitude of these complex values. Probably you would notice that all of them are very high. That's because energy is not normalised. Normalisation can be done via dividing values of magnitude by energy of your window. As people tend to take signal "as it is" (which is equal to applying rectangular window) the spectrum is then divided by a number of samples $N$. For different types of windows this factor that is equal to the sum of all window samples.

Additionally, we are only interested in one half of a spectrum, therefore the amplitude of all samples must be multiplied by $2$ to compensate the loss of energy, except DC component (which appears only once) and Nyquist bin (if it's present).

Lastly, as you want to analyse RMS of the signal, you must understand the following. All values of single-sided spectrum are of height $\dfrac{A_{i}^{2}}{2}$, which is equal to $ \left( \dfrac{A_{i}}{\sqrt{2}}\right)^{2}$, where $\dfrac{A_{i}}{\sqrt{2}}$ is the RMS magnitude amplitude for the i-th frequency component. Translating this into Python code you obtain:

import numpy as np
import matplotlib.pyplot as plt

fs = 5e5       # Sampling frequency [Hz]
duration = 10  # Duration of the test signal [s]

N = int(fs * duration)
t = np.arange(N, dtype=float) / fs

f = 1000  # Frequency of the test signal [Hz]
sensitivity = 0.004  # Sensitivity of the microphone [mV/Pa]
p_ref = 2e-5 # Reference acoustic pressure [Pa]
amp = sensitivity * np.sqrt(2)  # Amplitude of sinusoidal signal with RMS of 4 mV (94 dB SPL)
signal = amp * np.sin(2 * np.pi * f * duration * t)  # Signal [V]

# Calculate the level from time domain signal
rms_time = np.sqrt(np.mean(signal**2))
db_time = 20 * np.log10(rms_time / sensitivity / p_ref)

# Apply window to the signal
win = np.hamming(N)
signal = signal * win

# Get the spectrum and shift it so that DC is in the middle
spectrum = np.fft.fftshift( np.fft.fft(signal) )
freq = np.fft.fftshift( np.fft.fftfreq(N, 1 / fs) )

# Take only the positive frequencies
spectrum = spectrum[N//2:]
freq = freq[N//2:]

# Since we just removed the energy in negative frequencies, account for that
spectrum *= 2
# If there is even number of samples, do not normalize the Nyquist bin
if N % 2 == 0:
    spectrum[-1] /= 2

# Scale the magnitude of FFT by window energy
spectrum_mag = np.abs(spectrum) / np.sum(win)

# To obtain RMS values, divide by sqrt(2)
spectrum_rms = spectrum_mag / np.sqrt(2)
# Do not scale the DC component
spectrum_rms[0] *= np.sqrt(2) / 2

# Convert to decibel scale
spectrum_db = 20 * np.log10(spectrum_rms / sensitivity / p_ref)

# Compare the outputs
print("Difference in levels: {} dB".format(db_time - spectrum_db.max()))

plt.semilogx(freq, spectrum_db)
plt.xlim( (1, fs/2)  )
plt.ylim( (0, 120) )
plt.grid('on')
plt.xlabel("Frequency [Hz]")
plt.ylabel("SPL [dB]")

plt.show()





  Difference in levels: -1.4210854715202004e-14 dB


As you can see, this is producing consistent results for both time and frequency domain. Good luck!