Quantization SNR of sine wave doesn't match 1.761 + 6.02 * Q
Some issues here:

Your SNR formula only applies to full scale sine waves, your sine wave has -6dB amplitude so your SNR will be 6 dB lower
The formula also implies rounding, not truncation, that's another 6 dB
You use a frequency that's a small integer divider of the sample rate, that means you are just repeating the same samples over and over again and don't get enough sample coverage to get a statistically meaningful result.
Your SNR analysis in the frequency domain is needlessly complicated and prone do imprecision and masking errors. Just do it directly in the time domain.

Here is how this would look in Matlab
%% quantization noise of a 16-bit sine wave
fr = 975.3; % something odd
n = 8192;
quant = 2^15;
% make the sine wave
y0 = sin(2*pi*(0:n-1)'/n*fr);
% quantize
yq = round(quant*y0)/quant;
% noise
yNoise = yq-y0;
% SNR
fprintf('SNR = %6.2fdB\n', 10*log10(mean(y0.^2)/mean(yNoise.^2)));

Technically, you would also have to work around the fact that the positive max amplitude of the sine might clip, but for large quantization this makes no meaningful difference.
EDIT
It may be good to review where the formula comes from and what it actually means. It all starts with the quantization noise. If we quantize and round than the quantitation noise is uniformly distributed between $[-0.5 \delta ,-0.5 \delta ]$, where $\delta$ is the quantization step. For truncation it would be uniformly distributed $[0,\delta ]$. For rounding, the resulting noise powers are
$$ P_{round} = \frac{\delta ^2}{12}, P_{trunc} = \frac{\delta ^2}{3}$$
If $B$ is the number of bits, then for a signed signal, we simply have $\delta = 2^{B-1}$ and hence we get the noise level for 16 bit rounding as
$$P_{16} = \frac{2^{-30}}{12} \Rightarrow -101.1 dB $$
That's true for ANY signal that has a reasonable wide distribution of samples. Since a full scale sine wave has a power of $-3dB$ the resulting SNR would be $-98.1 dB$