What's the standard deviation in bond angle for water molecules at a particular temperature? Suppose I have a beaker of water at 20 °C. I want to know the standard deviation of the bond angles between the hydrogen atoms in the water molecules. How do I figure this out?

If I knew that the bond angle was the only possible degree of freedom of the molecules, then I'd be able to use the Boltzmann distribution to calculate this, if I knew the energy of the molecule as a function of bond angle. But I don't know the energy of the molecule as a function of bond angle, and also there are a lot of other degrees of freedom (e.g. the molecules can have overall kinetic energy, or they can be spinning rapidly, or the other vibrational modes of the water molecule can contain some of the energy.

Is there an easy way for me to look up the energy of the water molecule as a function of the positions of the hydrogen atoms? (I guess I could just calculate this myself using my favorite quantum chemistry software.)
I tried to answer by proposing a simple model. The model didn't work as well as I hoped, but the general procedure might serve as an answer.

A harmonic model

As a first approximation, you might consider a simple model for the energy of a single water molecule as a function of its angle. One possibility might be a harmonic function

$$E(\theta) = \frac{k \theta^2}{2}$$

where $k$ is an unknown parameter and $\theta$ is the angle. This is very simple and leaves out many complexities, of course.



If you had $E(\theta)$, the energy of a single water molecule at a particular angle, and if you assume an independent particle model (which may break down for a liquid), you could apply Boltzmann weighting

$$p(\theta) = \frac{1}{U}\exp\left(-\frac{E(\theta)}{R T}\right)$$

where $p(\theta)$ is the probability of finding a water molecule at a particular angle and $U = \int_0^\pi d\theta \exp\left(-\frac{E(\theta)}{R T}\right)$ is a normalization factor.



If you had a probability distribution $p(\theta)$ describing the spread of angles at a certain temperature, getting the mean and standard deviation is easy:

$$\langle\theta\rangle = \int_0^\pi d\theta\ \theta\ p(\theta)\\
\langle\theta^2\rangle = \int_0^\pi d\theta\ \theta^2\ p(\theta)\\
\sigma_\theta^2 = \int_0^\pi d\theta\ \left(\theta - \langle\theta\rangle\right)^2\ p(\theta) = \langle\theta^2\rangle - \langle\theta\rangle^2$$

So this would be how to get $\sigma_\theta^2$.
Below is my numerical attempt.
Beware: this model is too simple.



Numerical results

I did some calculations (mostly with SymPy) with the above and I came up with the following:

$$p(\theta) = \frac{1}{U} \exp\left(-\frac{k \theta^2}{2 R T}\right)\\
U = \sqrt{\frac{\pi R T}{2 k}} erf \left(\frac{\pi}{2} \sqrt{\frac{2 k}{R T}} \right)$$

$$\langle\theta\rangle = \frac{k}{R T U}\left(\exp\left( \frac{\pi^2 k}{2 R T} \right) -
 1 \right) \exp\left( -\frac{\pi^2 k}{2 R T} \right)\\
\sigma_\theta^2 = \frac{R T}{k} - \left[ \frac{\pi}{U} \exp\left( -\frac{\pi^2 k}{2 R T} \right) +
\left(\frac{RT}{kU}\right)^2 \left(\exp\left( \frac{\pi^2 k}{2 R T} \right) -
 1 \right)^2 \exp\left( -\frac{\pi^2 k}{R T} \right)
\right]$$

By substituting constants and temperature ($\pi = 3.14159265$, $R = 8.314459848$ J K$^{-1}$mol$^{-1}$ and $T = 293.15$ K) and setting the average above to as close as possible to 104.45° and solving numerically the complete least squares problem, I found an approximate solution for $k$: $384.05105591$ J/rad².

Sadly, this value is associated with an average angle of around 96°, which means the model is very poor.

Using this best $k$ gives us a $\sqrt{\sigma_\theta^2}$ of 

$$\sqrt{0.508770928037803} = 0.713281801280393 \text{ rad}$$

which is around 41°. Again, a too simplistic model: the best prediction gives angles as 96°$\pm$41°!



A bad model

By expanding the average in series, I found the following

$$\langle\theta\rangle \approx 1.570796325 + 0.000530047067361739 k + O(k^2)$$

The second order term above was $- 7.86976333117302\times 10^{-7} k^2$, so I dropped it.
Observe the first term is pretty much 90° in radians. All other terms decay too quickly.

By doing the same expansion for the standard deviation the following comes out:

$$\sigma_\theta^2 \approx - \frac{4.54747350886464 \times 10^{-13}}{k} + 0.411233515772252 + 0.000388544793228128 k + O(k^2)$$

The next term in the series above was $- 4.34803406790969 \times 10^{-7} k^2$ so I dropped it. Again, a (wrong) zeroth order term dominates.