I guess for starters, try finding a torsion angle on a 3 point molecule.
Is it fast enough, accurate enough?
The real issue is why do some force-fields have intra-molecular bonds and others do not.
Including bond and angle forces in water models make their computation more expensive. There are classes of flexible forcefields for water. There are even more indepth forcefields that are polarizable.
Why don't we use the most sophisticated forcefield models for water?
The more sophisticated the model, the longer our simulations take. Often, rigid water models are good enough for certain properties. So, we use the simplest model that gives us the right answer for the property we are after. This is just a general good rule.
Parameters
Force-fields use parameters, and proteins use ALOT of parameters. These parameters are often fit for a specific water model. Hence, you are kind of stuck using that water model if you are using that Force-Field for the protein/solute molecule.
Complexity
As our interests move towards more involved calculations that involve water it does become important to include more effects such as flexibility and polarizability. This is increasingly the case in calculating Free Energies, whether it is the Free Energy of hydration or the Free Energy of binding a ligand to a protein. Another example is modelling supercritical water - only polarizable models have been successful. Rigid models such as SPC/E and TIP3P simply are not good enough for many properties under many conditions. You can of course fudge parameters to make anything fit a particular problem, but we want to use generalizable force-fields as much as possible. In this case, models with more parameters (flexible, polarisable) give us additional accuracy for the difficult problems.
Bandwagon
Another reason many people use one model more than another is because their colleagues use it, so they use it because other people use it and it likely won't be questioned since other people use it. This is not a good reason for using a certain model! usually the first people to use a certain model had a reason, it solved their problem. However often the following bandwagon users don't read the strength/weaknesses of the model but instead just use it because they see it used lots for seemingly similar areas of research. Welcome to molecular modelling.
Mixing Rules
Generally atoms have dispersion and electrostatic terms. In the case of dispersion, when the force or energy between two atoms is being calculated a single parameter for each term is required for the calculation. Take the well known Lennard-Jones for instance
$$U(r_{ij}) = 4 \epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}}\right)^6 \right]$$
different forcefields may have different numeric values for $\epsilon_i$, $\epsilon_j$, $\sigma_i$,$\sigma_j$ but they will all have a value. You need to use mixing rules to take the respective $i$ and $j$ parameters and create an $\epsilon_{ij}$ and $\sigma_{ij}$. We call these mixing rules a.k.a. combining rules. Wikipedia lists several of the common ones. https://en.wikipedia.org/wiki/Combining_rules
For instance for the Lorentz-Berthelot
$$\epsilon_{ij} = \sqrt{\epsilon_{i} \epsilon_{j} }$$
and
$$\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}$$
Finally, in SPC and SPC/E the hydrogens do not have Lennard-Jones parameters. It is not uncommon for hydrogens to have no dispersion parameters. This can lead to problems though, since dispersion parameters keep other atoms/particles from getting too close.
So a specific example of an SPC water molecule interacting with a Sodium ion, the electrostatics just needs to know that the Oxygen has a negative charge of $q_o =−0.82$ and the sodium has a charge of $q_{N_a} =+1$. A naive approach would be to treat this as a bare coulomb interaction which is just
$$ U(r_{ij})^{\rm electrostatic} = k \frac{q_o q_{N_a}} {r_{ij}}$$
where
$$k = \frac{e^2}{4 \pi \beta}$$
$e$ is the charge of an electron and $\beta$ is the permittivity and a global and constant value for the simulation.
A better approach is to use EWALD's summation for electrostatics. It is more complicated, but the only charge specific parameters it needs to know is still just $q_o$ and $q_{N_a}$. The other parameters it uses for damping can be considered global parameters. Note: also do this for each of the hydrogen charges interacting with the Sodium ion.
For Lennard-Jones you would need to use mixing rules for the oxygen (there is no LJ parameters for the hydrogens)
$$\sigma_{(o-N_a)} = \frac{\sigma_o + \sigma_{N_a}}{2}$$
and
$$\epsilon_{(o-N_a)} = \sqrt{\epsilon_o \epsilon_{Na}}$$
then do the standard LJ calculation for force and or potential energy.
