How can I create a histogram showing the RMS amplitude for each frequency bin of an audio file
If what you want is to determine the amount of energy your signal has in particular frequency bands you can perform a FFT of your signal. For that you need to make sure that you have amplitude data and not energy data. So you would want to have sound pressure Pa and not Pa^2. After having calculated a FFT with that the values have to be multiplied by $1/\sqrt{2}$ to get to the RMS-format.
Depending on your time data length and sampling time you get a specific frequency resolution of your signal in the frequency domain. With $N$ time samples and a sampling time of $\Delta t$ you get a frequency resolution of $1/(\Delta t \cdot N)$ Herz. Now if you want to know the energy contribution of a specific frequency range you can calculate the RMS-value of that frequency range using the earlier calculated linear spectra in RMS-format. According to Parseval's theorem one can use the following equation to calculate the RMS-value of a specific frequency range:

Therein the $A_k$s are the amplitudes of different frequencies in the spectra:

I used that approach to calculate the Energy of 10 Hz bins of a signal containing two sines with amplitude $2$ and Frequencies $5$ and $45$ Hz. Sampling frequency is 0.0098.
Here is the time Signal:

The resulting linear Spectra with Peak to Peak values is:

From that I calculated the RMS-Values of $\Delta 10$ Hz bins. So the first bar indicates the RMS-Value of the signal from $0 - 10$ Hz. The second from $10 - 20$ Hz and so on.

If some noise is added the spectra changes to:

The RMS-values of the $\Delta 10$ Hz bins become then:

I further suggest looking at Root Mean Square and Overall Level
Here is the Matlab Code I used to do the above (probably still containing errors though):
clear all
close all

t=linspace(0,10,1024); % time vector
y=2*sin(5*2*pi*t)+2*sin(45*2*pi*t)+randn(1,length(t)); % time signal
Fs=1/(t(2)-t(1));
Fn=Fs/2;
N=length(y);
yf=fft(y,N);
df = Fs/N; % Frequency Resolution

plot(t,y)
grid on
title('Time Signal')
ylabel('Amplitude')
xlabel(['Time Resolution: ',num2str(1/Fs),' s'])

amplH = abs(yf);
amplitudengang = fftshift(amplH/N);

x_fn = 0 : df : Fn-df;
x_fa = 0 : df : Fs-df;

figure
stem(x_fa-Fn, amplitudengang, 'b.-')
axis([-Fn Fn 0 max(amplitudengang)])
title('Two sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),' Hz'])
grid

amplitudengang=[amplitudengang(513) amplitudengang(514:end).*2];

figure
stem([0:df:(Fn-df)], amplitudengang, 'b.-')
axis([0 Fn 0 max(amplitudengang)])
title('Single Sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),'Hz'])
grid

amp10=amplitudengang(1:round((10-df)/df));
amp20=amplitudengang(round((10-df)/df):round((20-df)/df));
amp30=amplitudengang(round((20-df)/df):round((30-df)/df));
amp40=amplitudengang(round((30-df)/df):round((40-df)/df));
amp50=amplitudengang(round((40-df)/df):round((50-df)/df));

rms10=amp10./sqrt(2);
rms20=amp20./sqrt(2);
rms30=amp30./sqrt(2);
rms40=amp40./sqrt(2);
rms50=amp50./sqrt(2);

rms10=sqrt((rms10(1)/2)^2+sum(rms10(2:(end-1)).^2)+(rms10(end)/2)^2);
rms20=sqrt((rms20(1)/2)^2+sum(rms20(2:(end-1)).^2)+(rms20(end)/2)^2);
rms30=sqrt((rms30(1)/2)^2+sum(rms30(2:(end-1)).^2)+(rms30(end)/2)^2);
rms40=sqrt((rms40(1)/2)^2+sum(rms40(2:(end-1)).^2)+(rms40(end)/2)^2);
rms50=sqrt((rms50(1)/2)^2+sum(rms50(2:(end-1)).^2)+(rms50(end)/2)^2);

figure
stem([5 15 25 35 45],[rms10,rms20,rms30,rms40,rms50])
axis([0 Fn 0 max([rms10,rms20,rms30,rms40,rms50])])
title('RMS values per 10Hz')
ylabel('Amplitude')
xlabel(['Frequency in Hz'])
grid on

Edit
To reduce leakage the time data needs to be windowed. After that the resulting spectra needs to be multiplied by the matching correction factor: