Thank you to whoever was looking this closely at the code!

The term in Sheffield is

$\sum_j Z_j^2 n_j / \sum_k n_k Z_k$

if we divide both top and bottom by the total density $\sum_x n_x$ we get

$\sum_j Z_j^2 F_i / \sum_k F_k z_k$

where $F_i =n_i / \sum_x n_x$ is the fraction of each ion species.

Where, using the definition of zbar earlier in our code as
$\bar Z = \sum_k F_k Z_k$

we get
$\sum_j f_i Z_j^2 / \bar Z$

The $f_i$ is already included in our code, so I think the user's suggested code is correct.

I think we should change it as suggested. I will open a PR.

EDIT: math