problem with FFT in PYTHON
In most implementations the FFT returns the following for the DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j 2\pi nk/N}$$
Which would result in the total output being scaled by $N$.
The FFT implementation may or may not be scaled by N so it is always good to check the documentation directly for the fft that is being used. The total power is the sum square of all the FFT bins for a full FFT of length N where the time series data is of length N, once proper scaling is accounted for an assuming no windowing beyond the implied rectangular window is used. (I refer to another post that covers how to properly handle the windowing for power estimation).  So for the result as given above, the total signal power would be (and assuming units of voltage or current in the time domain, with a normalized resistance of 1 then the units of the output would be watts):
$$P = \sum_k\bigg(\frac{X[k]}{N}\bigg)^2$$
Besides reading the documentation, a quick test is to do fft([1, 1, 1, 1]) where fft() refers to whatever fft function is used; if the function needs to be scaled by $N$ before computing the total power, the result for above would be [4,0,0,0]. If the result was [1,0,0,0] then the DFT would have already been scaled by $N$ or another option is the FFT is scaled by $1/\sqrt{N}$ so that the fft and ifft are symmetric operations and the result would be [2,0,0,0].
Tne numpy fft function the OP is using rfft only returns the positive half spectrum, which would be the complex conjugate of the negative half spectrum for real signals and handles odd and even $N$ differently in terms of how many samples are included. For a consistent result, I recommend using fft which will return $N$ samples scaled by $N$, in which case to estimate the power you would first divide the fft results by $N$ and then sum the square of all the bins (consistent with Parseval's Theorem).
Also for Python it is recommended to use the fft in scipy.fftpack and not numpy (which is now in both for reverse compatibility reasons but is expected to eventually be removed from numpy; numpy is to be used for the lower level vector processing and broadcasting operations that scipy depends on but not the higher level numerical processing algorithms that are now all in scipy).
If windowing is applied, then the power computation needs to consider the coherent and non-coherent gain of the window. This is straightforward for a single tone whose power spectral density only occupies one bin; it will be reduced by the coherent gain of the window which is given by the sum of the window coefficients and scaled by $N$ (which is always a loss). However for signals with spectral content that spans multiple bins, we must be careful to consider that the bandwidth of a windowed DFT bin will always span multiple bins, so after the decrease due to coherent gain, there will be an increase due to the double counting as all the bins are summed in the power computation previously given. The net result is a processing gain that is applicable to signals with spectral densities spanning multiple bins and described further in the links below.
Find the Equivalent Noise Bandwidth
Effect of windowing on noise