Structure of Simple Dipolar Water-Like Fluids: Primitive Model and Hard Tetrahedra

Dipolar versions of two qualitatively different types of simple short range model fluids which exhibit the phenomenon of hydrogen bonding and which could thus serve as a reference in equations of state for associating fluids have been considered: the primitive model of water descending from the TIP4P model and the fluid of hard tetrahedra. The hydrogen bonding structure exhibited by the latter model results from purely repulsive interactions whereas in the first model the “hydrogen bonding interaction” is explicitly incorporated in the model. Since the water molecules bear a strong dipole moment, the effect of the added dipole-dipole interaction on the structure of the two short-range models is therefore examined considering them both in the full and screened dipole-dipole modifications. It is found that the hydrogen bonding structure in the primitive model resulting from the site-site interactions is so strong that the additional dipole-dipole interaction has only a marginal effect on its structure and contributes thus only to the internal energy. On the contrary, even only a weak dipole-dipole interaction destroys the original hydrogen bonding structure of the hard tetrahedron fluid; to preserve it, a screened dipole-dipole interaction has to be used in the equation of state development.


INTRODUCTION
Simple models have played a fundamental role in our understanding of molecular mechanisms governing the behavior of fluids and for the development of tractable expressions of their thermodynamic properties. The use of simple models goes as far back as to van der Waals (vdW) who used, intuitively, the notion of the excluded volume (hard body model) for the development of his famous equation of state (EoS). It took nearly one century before the vdW equation was laid on a rigorous theoretical footing using a perturbation expansion which is the only theoretical tool to deal with complex realistic interaction models. It is based on results of molecular simulations carried out during the 1960s and which showed that the structure of normal (i.e., non-polar) fluids can be well estimated by that of appropriate purely repulsive hard body fluids which is a necessary condition for the perturbation expansion about a hard body reference to converge (Barker and Henderson, 1976). Further extensive simulations in 1990's and at the beginning of 2000's then extended this finding also to complex fluids: It was shown that the structure of fluids is determined, in general, by short-range interactions which however may be not only repulsive (which is the case of normal fluids) but also attractive (in the case of polar and associated fluids); for a review see (Nezbeda, 2005). These findings have extended the potential of the perturbation expansion to derive a molecular-based EoS for the entire class of fluids and write such an EoS in the form which results from a decomposition of the considered intermolecular interaction model (force field), u (1, 2), into short-range reference, u SSR (1, 2), and perturbation, u pert (1, 2), parts, u(1, 2) u SSR (1, 2) + u pert (1, 2) In the above equations symbol (1,2) stands for the complete set of generalized coordinates of molecules 1 and 2, P is the pressure, β is the inverse temperature, β 1/k B T, where k B is the Boltzmann's constant and T is the temperature, ρ is the number density, ρ N/V, z ref is the compressibility factor of the reference system, and Δz is a correction term.
The only long-range part of u (1, 2) is the Coulombic interaction. In all so far performed simulations on the effect of the range of interactions, the SRR has been obtained by a gradual switching off this interaction (for a review see (Nezbeda, 2005)), where S is a switch function, r 12 is the distance between reference sites of the molecules, and R′ and R″ are its appropriately chosen parameters (for details of the effect of various choices of R′ and R″ see (Nezbeda and Kolafa, 1999)). In addition to the constraint that the reference fluid has to reproduce the structure of the considered fluid, the choice of the reference is also subject to another constraint: the reference has to be amenable to a theoretical treatment resulting in close analytic expressions for its properties. And this is the point where simple models enter the process: The fluid defined by potential (3) remains too complex to provide analytic results and simple models are therefore used to approximate its properties.
Confining our considerations henceforth to associating fluids exemplified by water, its primary distinctive feature is its open structure resulting from hydrogen bonding (H-bonding) and this phenomenon has to be reproduced by any simple model considered for the reference in the perturbation approach. First such simple intuitive models emerged in the end of the 1980's (BoL, 1982;Dahl and Andersen, 1983;Smith and Nezbeda, 1984;Kolafa and Nezbeda, 1987). All of them employ a hard sphere with embedded attractive sites, using typically a squarewell interaction, and different models differ then in the number of the H-bonding sites and their geometry. Later on, such models (referred to as primitive models, PM) have been linked directly to realistic force fields (parent model) and constructed as their descendants using statistical mechanical tools (Nezbeda, 1997;Vlcek and Nezbeda, 2004a), see Figure 1. Essentially the same models, but without any link to realistic force fields, have also been used in SAFT equations (Vega and Llovell, 2016).
Besides the PM's where the H-bonding structure is forced explicitly by the attractive site-site interactions, a hard body model which also exhibits a water-like structure is the fluid of hard tetrahedra . This model is therefore another possible simple model to be used as a reference both in a perturbed molecular-based EoS or in a vdW-like EoS for water. This choice may be viewed as an associated fluid counterpart of hard spheres (hard body) used for non-polar fluids. Its properties have been the subject of intensive research (Haji-Akbari et al., 2011;Kolafa and Labík, 2015;Irrgang et al., 2017) and an EoS in a close analytic form is also available (Tian et al., 2019).
In order to complete the perturbed EoS development it is necessary to consider contributions of the Coulombic interactions screened at short separations, Since water molecules possess a strong dipole moment, in practical computations the leading term of the electrostatic interactions is (very accurately) approximated by the leading dipole-dipole interaction, u DD (1, 2). In the evaluation of Δz DD one should however bear in mind that, according to Eq. 5, the DD interaction should be considered only over a reduced intermolecular separations to avoid its double counting at close separations. This has been the case of molecular-based equations (Nezbeda and Pavlíček, 1996;Jirsák and Nezbeda, 2007a) but not of vdW-like equations for water (Vega and Llovell, 2016). Formally, the omission of this constraint may not cause problems. At the macroscopic level it may improve the FIGURE 1 | Schematic representation of the TIP4P geometry of water molecules and the primitive model. Sites X represent the hydrogen-like sites H moved on the surface of the oxygen sphere.
Frontiers in Chemistry | www.frontiersin.org December 2021 | Volume 9 | Article 783741 performance of the equations but, on the other hand, it also may turn such an EoS to a purely empirical one and at the molecular level such an equation may thus lose any justification. The primary motivation for this study has been the examination of the effect of adding the dipole-dipole interaction, either completely into the reference or in its screened form only, on the structure and, consequently, on the EoS development. We consider therefore two qualitatively different simple model fluids, i) the dipolar primitive model of water descending from the TIP4P (Jorgensen et al., 1983) force field with the explicit site-site interactions producing the H-bonding structure, and ii) the fluid of dipolar hard tetrahedra which itself, without the DD interaction, adopts an H-bonding-like structure solely due to purely repulsive interactions . The dipole-dipole interaction is considered both over the entire range and screened at short separations. Carrying out standard Monte Carlo simulations we focus on the structure described by the complete set of the site-site correlation functions. In the next Section we provide necessary theoretical and computational details and in the following Section results are presented and discussed. The main findings are summarized in the last Section along with a potential further development.

The Models
We are going to consider two types of short-range dipolar hardcore models, the fluid of hard tetrahedra and the primitive model of water descending from the realistic TIP4P model (Jorgensen et al., 1983).

Primitive Model of Water
Primitive models descend, in general, from realistic models by neglecting all medium (van der Waals) and long ranged (electrostatic) interactions, and account qualitatively only for the short-range, both repulsive and attractive interactions.
The parent TIP4P models have, in addition to the central uncharged oxygen site, two H-like sites bearing a positive charge and an M-site with a negative charge, see Figure 1. Denoting the sites which bear charges as P (positive charge) and N (negative charge), and the site without electrostatic interactions as O, the complete intermolecular interaction energy of the PM is given by [fore details see (Vlcek and Nezbeda, 2004a)].
where the summation in the second term of this equation runs over the pairs of like sites, in the third term over the pairs of unlike sites, and

Hard Tetrahedron Fluid
There are different conventions, and hence also different scalings, used to describe the geometry of tetrahedra, see Figure 2. Kolafa and Nezbeda (Kolafa and Nezbeda, 1995), and also Haji-Akbari et al. (Haji-Akbari et al., 2011), considered a regular tetrahedron inscribed in a cube of the edge length equal to 2a whereas Kolafa and Labik (Kolafa and Labík, 2015) then used tetrahedron's edge, h, to scale distances. For convenience, and to keep contact with the fluid of hard spheres (HS), we scale in this paper distances by diameter σ HS of the circumscribed sphere. Thus, a ( 3 √ /6)σ HS , h (2/3) σ HS , and the tetrahedron's volume is V ( 3 √ /27)σ 3 HS . Hard tetrahedron interaction potential is thus given by.

Dipolar Models
Dipolar versions of the above models are obtained by adding a point dipole of strength μ to the central site, Frontiers in Chemistry | www.frontiersin.org December 2021 | Volume 9 | Article 783741 whereR 12 is a unit vector in the direction of center-to-center vector R 12 , and e i is a unit vector in the direction of dipole i. The total interaction potential of the considered models assumes thus the form where u ref (1, 2) is either u HT (1, 2) or u PM (1, 2). However, as it is clear from Eq. 3, the formal addition of the dipole-dipole interaction to models which already include H-bonding is not fully correct; the dipole-dipole interaction has to be suppressed and the dipolar models should thus properly possess this modified form, Dipolar hard body fluids are characterized by a dimensionless energetic parameter λ which combines the strength of the dipoledipole interaction with the scaling parameter, where β 1/k B T, k B is the Boltzmann's constant and T is the temperature. The direction of the dipole moment in the PM is defined by sites O and M. Concerning hard tetrahedra, there is no a priori reason for a specific direction of tetrahedron's dipole moment. In  it was argued that the face-to-face packing of tetrahedra in dense systems resembles a hydrogen bond network in water and a geometric definition of a "hydrogen bond" was provided [for further details see ]. In this view the vertices of tetrahedron may be identified with the location of two hydrogen atoms and the lone electron pair (TIP5P-type geometry). We thus define the dipole moment as the vector from the center of the tetrahedron to the center of an tetrahedron's edge.

Simulation Details
For simulations with the site-site models it is convenient to follow the physical picture and represent the dipole by two charged sites with charges q i ±1 a distance l 0 apart; the dipole moment is then μ ql 0 . The closest approach of two PM molecules is σ HS and the value l 0 /σ HS was set to 0.01. The closest approach of two tetrahedra is σ HS /3 and a lower value of l 0 had therefore to be used, l 0 /σ HS 0.001. Test simulations were performed to ensure that the use of these values yielded, within the simulation errors, the same result as those with the point dipole.
To account for the long ranged character of the Coulombic interactions, the reaction field method (Neumann et al., 1984) was employed. The dipole-dipole interaction energy assumed thus the form.
where the summation runs over the charges on molecules a and b, the second term is the reaction field correction for the finite radius cutoff, r c , ε r is the dielectric constant of the surrounding continuum set to infinity in the simulations, and ϵ is related to the conventional parameter λ by the relation Standard simulations in an NVT ensemble (Frenkel and Smit, 2002) with 512 particles were performed at several densities and for a number of values of λ. Control quantities  were always computed and the evolution of the energy was monitored to follow the development of the system to be sure that the productive runs started from a properly equilibrated configuration. The structure of the considered fluids was the primary quantity of interest in this study. All necessary information is provided by the full pair correlation function g (1, 2) whose complete determination is however practically impossible. The usual way is to characterize the structure of fluids by a set of partially averaged, i.e., the site-site correlation functions g ij (Nezbeda and Smith, 1981), where δ is the Dirac delta distribution, and 〈. . .〉 denotes an ensemble average. Another quantity of interest obtainable from the center-to-center correlation function, g CC , is the coordination number, which differentiates associating fluids from non-polar ones; R min in this equation denotes the location of the first minimum of g CC . Low values of the coordination number (the number of molecules in the 1st coordination shell) within the range from 4 to 6 at liquid densities point to a typical water-like structural arrangement in the liquid.

RESULTS AND DISCUSSION
Using the methodology described in the preceding section we performed simulations for a series of values of parameter λ with the switching function whose parameters were taken from the previous studies (Nezbeda and Kolafa, 1999;Kolafa and Nezbeda, 2000) and for a series of densities with the main focus on the structure described by the complete set of the site-site correlation functions. When setting the parameters of the dipolar models we wanted to keep contact with real water. The dipole moment of the water molecule (gas value) is 1.85D. The size of the water molecule is not uniquely defined but its value used in various applications varies around 3Å. Confining the considerations to the room temperature we get then for parameter λ value around 3. When the liquid value of μ is used, μ ≈ 2.9D, then λ is around 7. For comparison, in their extensive study of the dipolar HS fluid Theiss and Gross (Theiss and Gross, 2019) reported thermodynamic functions up to λ 7, which was also the upper limit in our simulations.

Primitive Model
All parameters of the PM but the site-site strength, ϵ SW , are obtained from the parent TIP4P model using statistical mechanical tools and taken from (Vlcek and Nezbeda, 2004a). For ϵ SW /k B T Jirsak and Nezbeda (Jirsák and Nezbeda, 2007b) used the value 4440 K to reproduce the temperature density maximum. Another possibility is to consider the value 3300 K which yields the experimental energy of the dimer. In theoretical study (Vlcek and Nezbeda, 2004b), Vlcek and Nezbeda reported results for the PM up to β ≡ϵ SW /k B T 8. Since relating the PM to the full water model is not, in principle, correct (Nezbeda, 2020) and because of ambiguity of its value, we decided to use also the lowest temperature, β 8. As regards density, as a consequence of its relatively open structure, liquid water packing fraction is lower than that of, e.g., argon-like fluids. Majority of simulations were therefore carried out at a typical density, η 0.35, and these results are also reported here.
In previous studies on the effect of the long-range interactions on the structure of water (Nezbeda and Kolafa, 1999;Kolafa and Nezbeda, 2000) it was shown that the shortest cutoff guaranteeing the identity of the structures was the range < 4 A, 6 A > and the same values were therefore used also for R′ and R″ using for the HS diameter the same value, i.e., σ HS 2.653Å.
Changes in the structure in dependence on λ and the gradual inclusion of the dipole-dipole interaction at the given conditions are shown in Figure 3. As we see, the effect of the dipole-dipole interaction does not seems, perhaps surprisingly, significant with the results of the screened and zero DD interaction being barely distinguishable. The only noticeable difference we find in g OO for λ 7. The strong dipole-dipole attraction increases the population around the central particle but this effect has only a marginal impact on the coordination number which in all cases attains values about 3.95. Concerning the site-site correlation functions, g OH and g HH , they seem to remain intact so that also the orientational arrangement is not affected by the added dipoledipole interaction. Consequently, although the chosen temperature, β 8, is very low and corresponding to cold liquid water, the resulting H-bonding is so strong that the added dipole-dipole interaction does not seem to be able to break it.
Unlike the structure, the additional interaction has to have an impact on the energy which is shown in Table 1. The contribution of the dipole-dipole interaction is rather small,  perhaps also surprisingly small, and becomes comparable with the H-bonding energy only after the full inclusion (no screening) for the highest value of λ, λ 7.

Fluid of Hard Tetrahedra
Despite its simple and compact shape, hard tetrahedron is an extremely nonspherical body. For comparison, its convex body shape factor α is 2.2346 which corresponds, e.g., to hard prolate spherocylinders of the length-to-breadth ratio greater than 5 (Boublík and Nezbeda, 1986). It is therefore instructive to compare first the center-center correlation functions of the fluids of hard spheres and tetrahedra at different densities; see Figure 4 where the distances are scaled by the respective closest approach distance. g cc of the HS fluid remains qualitatively the same and changes only numerically with increasing density with the coordination number remaining also the same, about 12, which is a typical number for nonpolar fluids, and the same applies also to the location of the minimum, r/r closest appr ≈ 1.6. g cc of the HT fluid is rather flat without a pronounced first peak (cf. Figure 5 of ). Although density η 0.35 falls, generally, into a higher medium liquid density range, the HT fluid at this density still behaves as a rare gas with the coordination number about 20. Structural changes thus take place over a relatively narrow range of high densities. At packing fraction η 0.425, which is close to the onset of metastability (Kolafa and Labík, 2015), the minimum of g CC moves to the same location as that of the HS fluid, r/r closest appr ≈ 1.6, and the coordination number drops to about 3. With respect to the above findings, namely that the hard tetrahedron fluid at medium density behaves as rare gas, we will present and discuss its dipolar versions only for the high density. As regards values of λ, a direct comparison between λ tetrahedron and λ HS is not possible. It is necessary to realize that the latter is scaled by the closest approach which is about three times larger than that for tetrahedra. The preferably used value in this paper, λ 1, represents therefore a very strong dipole-dipole interaction.
To examine the effect of the dipole-dipole interaction and its screening, it is also necessary to set the parameters of the switching function. In the case of the PM, the dipole-dipole interaction was completely switched off directly at the end of the hard core repulsion (r/σ HS 1) but at a distance, approximately, R′/σ HS 1.5 away with a relative width of the switching range 0.75. With no other information on the switching available, the same values have been used therefore as a guidance for setting the same range, < R ′ , R ′′ > , also for the tetrahedron fluid.
The water-like arrangement of tetrahedra is exclusively only due to their shape so that it is relatively fragile and it may be thus expected that any additional (perturbing) interaction may destroy it. To confirm this intuitive presumption, we show in Figure 5 the center-center correlation functions, g CC , for a series of λ. As it is seen, by adding a dipole, regardless of its strength, the structure of the fluid changes. In fact, the g CC function strongly resembles that of the dipolar HS characterized by a very high and narrow peak close to the contact (cf. Figure 4 of (Ganzenmüller and Camp, 2007)). The minimum also moves towards contact and then g CC becomes rather flat with an indistinctive second maximum, unlike the HS case which continues to exhibit sharp peaks at multiples of the HS diameter. The coordination number drops down to about 2. In contrast to the HS case, for hard tetrahedra away from the contact the DD interaction competes with their shape which suppresses the exclusive DD interaction dominance and hence also no further sharp peaks are observed.
The low coordination number deserves a discussion. An analysis of configurations reveals that, typically, two polar tetrahedra are arranged face-to-face around the central tetrahedron which explains the very high peak of g CC close to contact, see Figure 6B. There can we also see two other tetrahedra in vicinity of the central body which would mean that there are (may be) four molecules around the central one which would correspond to the structure of low temperature water. However, the two tetrahedra are oriented with their vertices pointing to the central tetrahedron so that their center falls slightly beyond the range determined by the minimum of the correlation function g CC . It thus seems that the definition of the coordination shells based on the minima of the center-center correlation function, which suits well for systems made up of not too nonspherical molecules, may be misleading for highly nonspherical molecules as, e.g., hard tetrahedra: some tetrahedra may fall into vicinity of the central tetrahedron by their vertex but their center may be then outside the range. For comparison we show in Figure 6A the surrounding of a nonpolar tetrahedron. As mentioned above, a typical gas-like chaotic distribution is observed.
Information on the orientational arrangement is provided by the complete set of the site-site correlation functions g CV and g VV . These functions, both with the screened and full DD interaction for λ 1 are shown in Figure 7 and Figure 8. Since the model with the dipole moment is asymmetric, the vertices are not equivalent. Denoting the edge connecting the charged vertices as 1-2 and the other two vertices as 3 and 4, then there are two different center-vertex correlation functions, [C-1] and [C-3], and 3 vertex-vertex correlation functions, (Barker and Henderson, 1976), (Barker and Henderson, 1976;Nezbeda, 2005;Nezbeda and Kolafa, 1999), and (Nezbeda and Kolafa, 1999).
There are two features of g CV and g VV which catch attention. First, in the full agreement with the claim that the structure is determined primarily by the short range interactions, we see that the dipole-dipole interaction switched on beyond the repulsive range has only a marginal effect. The C-V functions are nearly identical and the differences in the V-V functions are nearly negligible. Second, it is a relatively weak effect of the full DD interaction on the orientational arrangement (V-V functions). Except for the localized (and some sharp) maxima resulting from the interplay between the DD interaction and hard core geometry [e.g., hard tetrahedra g exhibits a cusp at r/σ 2 3 ] the correlation functions of the full and screened models follow practically the same course. The full DD interaction is even suppressed at short separations for the V1-V1 and V3-V3 pairs.

CONCLUSION
Simulations of the dipolar versions of the two short-range models, the primitive model of water and hard tetrahedron fluid, have revealed differences in their response to adding the dipole-dipole interaction.
The primitive model, which has built-in hydrogen bonding, seems nearly "immune" to this additional interaction. It means, the dipole-dipole interaction does not cause practically any change in the structure of the liquid. This is a very favorable/  important finding from the point of applications. Despite its long range, the dipolar version of the primitive model can be used directly as a suitable reference for the development of a molecular-based EoS. In the vdW-like approach the inclusion of the full Δz DD term represents then an acceptable approximation to the (theoretically correct) screened interaction. Since such a dipolar reference would capture most of the interactions and their contribution, it is expected that its amendment by a contribution of dispersion forces (Aim and Nezbeda, 1983) could yield, in both approaches, a good equation of state for water. Unlike the primitive models, the hydrogen bonding structure exhibited by the hard tetrahedron fluid results from the purely excluded volume effect and is thus sensitive to any other added interaction. As discussed in the preceding section, the dipolar tetrahedron fluid may exhibit rather a complex behavior different from that of the non-polar tetrahedron model and its formal use as a reference would be thus erroneous. However, when the dipole-dipole interaction is screened, the structure of this model fluid does not change and it may be thus also used as a reference for the development of a molecular-based equation of state.
The last finding opens a possibility to follow the original vdW way of thinking which has resulted in plethora of commonly used cubic equations (Valderrama, 2003). These equations, with the hard sphere reference term, have been extended also to associating fluids by adding an associating term. The resulting equations are referred to as Cubic-Plus-Association (Yakoumis et al., 1998). This purely empirical approach may be cast into a physically acceptable/correct one if the hard sphere term is replaced by that which reflects H-bonding. The hard tetrahedron fluid may serve this purpose and such an attempt will the subject of future investigation.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the author, without undue reservation.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication. FUNDING Support for this work was provided by the Czech Science Foundation (Grant No. 20-06825S).