Contribution of vibrational degrees of freedom
You have to be careful to make a distinction between the coordinates needed to describe a molecule, and the degrees of freedom which the molecule can actually access. For this purpose, one might choose to define the number of degrees of freedom as the minimum number of spatial coordinates needed to fully describe the positions of atoms relative to a fixed reference frame. This will be $3N$ total coordinates which corroborates the other answer you link. We then recognize that motion for non-linear (linear) systems, motion in six (five) of these coordinates correspond to translation and rotation. Therefore, the other $3N-6$ ($3N-5$) coordinates describe the vibrations of the system.

Now, if one is concerned with the total energy of a system, then we should keep track of the Hamiltonian. In a classical context, the Hamiltonian is a function of all coordinates which are capable of holding energy in one form or another, $H(\mathbf{q},\mathbf{p})$. For molecules this means we must keep track of $3N$ spatial coordinates, $\mathbf{q}$, and $3N$ momentum coordinates, $\mathbf{p}$.

Now, if we want to know about the average kinetic energy of atom 1 in the $x$ direction, we perform an ensemble average

$$
\langle T_{x,1}\rangle=\frac{1}{2m}\frac{\int p_{x,1}^2e^{-\beta H(\mathbf{q},\mathbf{p})}d\mathbf{q}d\mathbf{p}}{\int e^{-\beta H(\mathbf{q},\mathbf{p})}d\mathbf{q}d\mathbf{p}}
$$

It is a nice exercise to show that if all the degrees of freedom in the Hamiltonian are quadratic, then this separates into lots of integrals over gaussians and the top and bottom will cancel in all coordinates except for the $x$ direction motion of atom 1. This integral will give you the average energy as $1/2kT$.

But, since the integral over the quadratic momentum and the quadratic position as in the harmonic potential are exactly the same, we have just shown that every type of energy which depends quadratically on a coordinate has an average energy of $1/2kT$.

Therefore, if we have a particle in 1-D subject to a harmonic potential, the Hamiltonian is,
$$
H=\frac{p^2}{2m}+\frac{1}{2}kx^2
$$

which means the average energy is $kT$ since there are two coordinates which the Hamiltonian depends on quadratically.

So, hopefully now it is clear why we double to count of vibrations and arrive at the answer,
$$
\gamma=1+\frac{2}{13}=1.154
$$

I know this was very long-winded and overly detailed, but I hope this gives you a little better picture into where the equipartition theorem comes from and how this convenient counting of coordinates we do is really a coincidence of the fact a harmonic potential and kinetic energy are each quadratic in their respective coordinates.



As an aside, the distinction I have made between coordinates and degrees of freedom is completely arbitrary. Indeed, people often refer to both a harmonic potential and the kinetic energy as quadratic degrees of freedom. I prefer this language, but then one must be careful to clarify between degrees of freedom of position and degrees of freedom of momentum because this is more coordinates than are needed to describe the system completely, so it seems there are more degrees of freedom than there really are.