Do we have water models so far which are flexible? How good can they model specific properties of water and are they currently used with force fields and in actual simulations?
Allow me to begin with a disclaimer: what you have asked is a very contentious question among people who work on water models. I think much more contentious than it should be. Not necessarily because there's an obvious answer, but because there are many people who work on water models who really believe in their own model. The problem is that water is involved in almost everything and no model has been used in all those scenarios.
People use the MB-Pol potential, developed by the Paesani group, as if it gives you the correct answer for both static and dynamical quantities for water. Indeed, the potential comes very close to doing that for many observable properties and for reference calculations carried out at the coupled-cluster level of theory. The reason I phrased that first sentence as I did is because I have very little doubt that in ten years, there will be some new model which is somehow shown to be slightly more accurate than MB-Pol which becomes the new standard. This is simply the way of things in the water-modelling world.
It is almost certain that more models have been developed for water than for any other substance. This is because water itself has many unique properties which need to be understood and because the effect of water as a solvent is often necessary when trying to study some other material.
Because one is sometimes only interested in the dielectric properties of water and how they affect some other molecule, there are many rigid, fixed-charge models of water which rely on the Lennard-Jones potential (or something similar) to describe intermolecular interactions. The TIP and SPC families of potentials fall in this class. These are often used in biomolecular simulations. These are definitely not the most accurate water models made, so I will not discuss them further here.
I do not know the historical reasons for first creating rigid water models rather than flexible ones, but I suspect it is because introducing flexible bonds and angles causes the potential to be much more expensive to calculate during molecular dynamics (MD) since one essentially introduced many more degrees of freedom into the system. Also, having flexible bonds which vibrate at the high-frequency of an OH stretch greatly decreases the time step one can use in an MD simulation. (There are ways around this using multiple time-stepping algorithms). The most accurate models we have are most definitely flexible models, as these can describe the infrared spectrum of water, more accurately describe hydrogen-bonding, nuclear quantum effects, etc.
A Subtlety in Fitting Potentials: A detail which I think is important to mention is to make a distinction between water models which actually try to model the potential energy surface (PES) of water, and those which are trying to approximate some kind of mean free energy surface via a potential energy function. The TIP4P and SPC models fall into the latter category, as they are parameterized against experimental data. Since the actual potential energy values are not experimentally accessible in any meaningful way, virtually any model which uses experimental data in its parameterization is likely not modelling the PES, but actually has many complex effects "baked in". All the potentials I discuss from here on will be models of the actual PES of water. This means that, in principle, the dynamics of these potentials should be calculated using path integral simulations to capture nuclear quantum effects. Some papers make this mistake of comparing classical simulations of one of these models to experiment without addressing how quantum nuclei would change the answer.
In my eyes, there are two categories of water models which are both flexible and are true models of the PES. These are what I will call classical models and machine-learning (ML) models. The classical potentials are called as such because they are trying to model the quantum interactions of electron densities via some parameterizable classical model. these typically have only a handful of parameters. What I mean by classical is that the model is most often based on the mutlipole expansion in some way.
The best example I am aware of for this category are the Thole-type water models. The most popular of these are TTM2.1-F  and TTM3-F . The former of these has more accurate energetics and the latter of these has more accurate dynamics. I know it seems like that shouldn't be possible, but I promise it's true. I think there is an argument to be made (and which is somewhat ongoing in the literature) that TTM2.1-F is about as accurate energetically as any other existing water model, but it has too stiff of bonds, and thus does not reproduce the red-shifts of hydrogen bonded OH stretches properly. These models also both reproduce the opening of the HOH angle going from gas to liquid to ice. No classical model before TTM2-F had done this because the dipole should increase going from gas to liquid to ice, and this is not possible unless the charges on each atom are allowed to vary as a function of the bond length and angle.
Another popular family of classical potentials used for water is the AMOEBA model . This has been updated many times and, as far as I know, is quite popular in biomolecular simulations, as I think it is fast and can be mixed with other substances. I will describe in a minute a massive shortcoming of many potentials which is that they can only really describe water interacting with itself.
Another excellent feature of classical water models is that even though they attempt to be quantitatively accurate representations of the PES of water, they are still very fast to evaluate. On the other hand, some of the ML potentials I will discuss are much slower than these classical models (but much faster than DFT, for instance).
I am sure there are other classical water models worth mentioning, and I will likely update this answer as they get brought up.
Machine Learning Potentials:
Machine Learning (ML) models of water are the type of model being developed most actively. These models almost all try to do some kind of direct fit of the PES of water and hence tend to have literally tens of thousands of parameters. They tend to be slower to evaluate than classical potentials, but also generally more accurate (though it depends on the actual model).
As far as I'm aware the first of what I'm calling ML water models is the HBB2 potential from Joel Bowman's group . This model includes a spectroscopically accurate 1-body potential, an explicit fit of the 2-body interactions and the 3-body interactions among water molecules. The fit is done via permutationally-invariant polynomials, which at the time was not really considered an ML method, but many would now consider it to be one.
The model which built on HBB2 by generating more and more accurate training data is known as MB-Pol and is developed by Francesco Paesani's group. This model includes the same 1-body potential as HBB2, an explicit 2-body potential fit to 40,000 CCSD(T)/CBS water dimer geometries, an explicit 3-body potential fit to 12,000 CCSD(T)/aug-cc-pVTZ water trimer geometries, and a classical treatment of the 4-body and higher-order terms via the multipole expansion. This model has been shown to accurately reproduce many experimentally determined properties of water in all phases.
There is also a new crop of potentials being produced which are fitted via neural networks (NNs). While these are useful potentials in that they are the only ones discussed which can dissociate (whether or no they do this accurately is not something I've seen addressed in a convincing way). I won't go into much more detail as these are not the most accurate models specifically for water because NNs require massive amounts of data to give reliable results, so any ab initio NN potentials tend to be trained on DFT calculations. DFT functionals are often parameterized against experiments and/or CCSD(T) calculations, and DFT tends to not do well over the entire phase diagram of water, so I have no problem saying these potentials are not the most accurate ones thus far. They are useful if you need a dissociable model though.
I found quite an interesting paper recently which tries to survey the history of modelling water computationally which is called: "Modelling Water: A Lifetime Enigma". As someone who studies water, I believe this is a very fitting title.
There are two important features which I haven't addressed in detail. One is that water can dissociate. For most simulations of gas and condensed phases of water, dissociation doesn't really matter. It happens so rarely as to be completely negligible. However, there are some problems where the entire effect is predicated on water dissociating. These would include proton transport, water splitting, acid-base reactions, etc. So, the lack of highly-accurate dissociable water potentials is a problem. One recent model is known as Rex-Pon. I am not sure how accurate it is.
The second problem is the general utility of these models. The question specifically asked
are they currently used with force fields in actual simulations?
The answer to this is that it depends. By force fields, I assume you mean the simple molecular mechanics models which usually depend only on a couple of parameters and are common in biomolecular simulations. The problem here is that these potentials are not necessarily compatible with the water models described here. Most potentials need to live in their own world in order to give reliable results. Virtually all of the classical and ML potentials I described can only describe water-water interactions. I am pretty sure AMOEBA is a general potential (i.e. it can describe arbitrary atoms) as is Rex-Pon. The only problem is that general potentials tend to be less accurate, so one has to be very careful to know the limitations of your model and decide if it's useful for the problem at hand.
Another way to answer this question is to say that the whole question of modelling water is an absolute mess. In modelling water, everything is possible, but nothing is certain.
: Reddy, S. K., Straight, S. C., Bajaj, P., Huy Pham, C., Riera, M., Moberg, D. R., ... & Paesani, F. (2016). On the accuracy of the MB-pol many-body potential for water: Interaction energies, vibrational frequencies, and classical thermodynamic and dynamical properties from clusters to liquid water and ice. The Journal of chemical physics, 145(19), 194504.
: Paschek, D. (2004). Temperature dependence of the hydrophobic hydration and interaction of simple solutes: An examination of five popular water models. The Journal of chemical physics, 120(14), 6674-6690.
: Fanourgakis, G. S., & Xantheas, S. S. (2006). The flexible, polarizable, Thole-type interaction potential for water (TTM2-F) revisited. The Journal of Physical Chemistry A, 110(11), 4100-4106.
: Fanourgakis, G. S., & Xantheas, S. S. (2008). Development of transferable interaction potentials for water. V. Extension of the flexible, polarizable, Thole-type model potential (TTM3-F, v. 3.0) to describe the vibrational spectra of water clusters and liquid water. The Journal of chemical physics, 128(7), 074506.
: Ren, P., & Ponder, J. W. (2004). Temperature and pressure dependence of the AMOEBA water model. The Journal of Physical Chemistry B, 108(35), 13427-13437.
: Shank, A., Wang, Y., Kaledin, A., Braams, B. J., & Bowman, J. M. (2009). Accurate ab initio and “hybrid” potential energy surfaces, intramolecular vibrational energies, and classical ir spectrum of the water dimer. The Journal of chemical physics, 130(14), 144314.
: Partridge, H., & Schwenke, D. W. (1997). The determination of an accurate isotope dependent potential energy surface for water from extensive ab initio calculations and experimental data. The Journal of Chemical Physics, 106(11), 4618-4639.
: Ouyang, J. F., & Bettens, R. (2015). Modelling water: A lifetime enigma. CHIMIA International Journal for Chemistry, 69(3), 104-111.
: Naserifar, S., & Goddard, W. A. (2019). Liquid water is a dynamic polydisperse branched polymer. Proceedings of the National Academy of Sciences, 116(6), 1998-2003.