{'fraction_non_alphanumeric': 0.035841848619436595, 'fraction_numerical': 0.009938814175233717, 'mean_word_length': 4.59867849069727, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular mean field theory for liquid water

Attractive bonding interactions between molecules typically have inherent conservation laws which influence the statistical properties of such systems in terms of corresponding sum rules. We considered lattice water as an example and enunciated the consequences of the sum rule through a general computational procedure called"Molecular mean field"theory. Fluctuations about mean field are computed and many of the liquid properties have been deduced and compared with Monte Carlo simulation, molecular dynamics and experimental results. Large correlation lengths are seen to be a consequence of the sum rule in liquid phase. Long range Coulomb interactions are shown to have minor effects on our results.


I. INTRODUCTION
Molecules in bulk interact with each other consistent with the underlying quantum mechanical dynamics, for example, hydrogen bond, variety of sulphur bonds etc. These interactions, though quite strongly attractive compared to ambient thermal energy, still fluctuate only when the concerned molecules are very close to each other (of the order of Angstrom).
These molecular interactions are countable and categorizable. For concreteness, consider water molecules. Each water molecule has two hydrogen arms and two lone pair arms. Interaction between water molecules is via these oppositely polarized bond arms and when a hydrogen arm of one molecule bonds with lone pair of a neighboring molecule, it results in "hydrogen bond". So, for any number of water molecules (ρ) the number of hydrogen bonds (HB) and number of dangling bonds (DB) i.e., lone-pair and hydrogen arms which are not hydrogen-bonded, satisfy a "sum rule". It states that the arms of water molecule are either hydrogen-bonded or not i.e., DB + 2 HB = 4ρ (1) Now, if we consider a bulk of water molecules the above equation still holds when DB, HB and ρ are appropriately defined per unit volume. In other words, the local topology of molecular interactions implies a sum rule which is also true in the bulk for any thermodynamic conditions such as temperature, pressure, volume. Furthermore, this is also independent of other interactions in the dynamical system such as van der Waals' (vdW), Coulombic etc.
These facts are not surprising since Eq.(1) is a topological constraint which is insensitive to many details of the dynamics.
In most descriptions of solid phase (ice or glasses), by design, the above sum rule is maintained explicitly. In the gas phase where the density of hydrogen bonds is negligible the sum rule is again trivially satisfied. If there are significant HB then the system can be analyzed by introducing appropriate density of dimers which still satisfy the sum rule thermodynamically.
Liquid phase is typically studied using molecular dynamics (MD) and Monte Carlo (MC) simulations of popular models [4, 13-15, 19, 22, 31, 37, 40] where in again the sum rule is maintained at every step of the dynamics and many bulk and other correlation properties [36] are well reproduced and are consistent with experiments [27,34,35,41].
Liquid phase is also described phenomenologically by a mean field which invariably is taken to be density field and fluctuations there of are then postulated to obey Ornstein-Zernike (OZ) equation [3,8,12]. Furthermore, all molecular interactions have the property that they cannot approach each other closer than certain minimum distance, usually incorporated as hard sphere repulsion in the interaction potential between the molecules. This in the liquid phase gives us the well known vdW term in the equation of state, often referred to as density saturation effect [12]. Phenomenologically, such saturation effects can be incorporated in the OZ equation and thereby understand important components of density fluctuations namely coordination shell peaks using approximations like Percus-Yevick, hypernetted-chain etc [12]. But these descriptions are not as successful as MD simulations, primarily because by design in the construction of phenemonological models one does not envisage fluctuations in HB or DB densities. Consequently, they cannot satisfy the sum rule in large subvolumes where surface effects are negligible. The sum rule argument suggests that to describe water we need atleast two fluctuating fields to be self-consistent.
Some important developments in the past on systems with hydrogen-bond' like interactions are discussed below. Wertheim's work on statistical thermodynamics of fluids with directional interactions consists of systematic graph-theoretic expansions of n-mer tree-like configurations and provides OZ-like operator equation in terms of total density and singletdensity fields of the system [42][43][44]. Another significant work on associating fluids was by Andersen who provided a Mayer-like density expansion technique to deal with directional interactions [1]. Both the descriptions work very well at high temperatures and also at densities corresponding to liquid-vapor coexistence. An extension to Wertheim's theory was proposed as SAFT [9,25] wherein the equation of state is phenomenologically extended to include contributions from complex geometries of constituent molecules which interact through directional forces. Another phenomenological description aimed at describing mainly the low temperature anamolies in fluids with characterstic orientation-dependent interactions was given by Truskett et al [39]. Here in again, the equation of state is approximated with vdW like terms which are then phenomenologically supplemented with hydrogen-bonding contribution that captures certain high-density features.
In this paper we propose to study a simple model Hamiltonian which incorporates the molecular hydrogen-bonding interactions. To accomodate the hard sphere repulsion we envisage the system on a hypercubic lattice and the model is essentially a slight generalization of the Pauling's model of water [30]. The partition function corresponding to the lattice model is analyzed by introducing appropriate discrete lattice fields. It is shown that the sum rule is automatically true in the bulk. 'Molecular mean field' (MMF) approximation extremizes the partition function in terms of the defined auxiliary fields. In addition, all the observables such as ρ, HB and DB densities are functionals of these auxiliary fields. One of the mean field equation which implies sum rule also implicates the "equation of network" i.e., a relation between equilibrium densities of HB and ρ. We study the equation of state and various mean field fluctuations in terms of correlation functions. We also considered long range Coulomb interaction and studied its consequences. Subsequently, an MC simulation study of the model is pursued and compared with the mean field theory quantitatively. We also discussed results of our analysis in the context of experiments and MD simulations. A brief discussion on future directions is presented at the end.
which imply that every water has two H arms and two L arms only. A hydrogen bond is realized when two water molecules two lattice units apart have their H and L arms meet at a site, as shown in Fig.(1).
At any site and on any link on the lattice the possible configurations are as shown in Fig. (1). We assign an energy −ν for every unpaired arm (H or L), and −λ for every hydrogen bond. We disallow for any two H (or L) arms to meet at the same site; also, more than two arms are disallowed to meet at a site as in Fig.(2).
To implement the constraints as shown in Fig.(2) in our analysis it is useful to define two where e α is a lattice unit vector in the direction α with the property e α = −e α . The discrete field b(x) counts the number of non-zero arms (approaching arms) in the neighborhood of site x, while q(x) measures the charge i.e, the difference between number of H arms and L arms meeting at site x. By construction, b(x) varies between 0 and 6 on a three dimensional hyper-cubic lattice and by imposing the condition that b(x) ≤ 2 in our analysis we ensured that no more than two arms can meet at a site. q(x) in turn varies between −b(x) to b(x).
In terms of these variables the Hamiltonian is : with constraints, The range of q(x) follows from Eq.(4). The Kronecker delta function denoted here as δ(p, q) is defined as δ(p, q) = 1 for p = q and 0 otherwise.
We now write the grand canonical partition function for our system at a finite chemical potentialμ for water and inverse temperature β: where we have introduced two auxiliary fields η(x) and φ(x) to impose the constraints Eq.(4).
The discrete fields η and φ take integer values in the range [−N, N] at every site, where N is any suitable large integer (greater than 8). The prime over summation refers to sum being restricted to constraints Eqs. (2,3,6). The introduction of auxiliary fields η(x) and φ(x) allows summation over other discrete fields W (x), H α (x), b(x), q(x) within their respective allowed range at each site x without any restrictions from the neighborhood configurations i.e., as if a single site functional Z site .
The Z site expression thus obtained is stated below. For brevity in the following expression π N η is written as η, π N φ as φ and further using fugacity definitions ν ≡ exp(βν), λ ≡ exp(βλ) and µ ≡ exp(βμ), It is useful to note the following identity in implementing this computation.
The observables in the system such as water density ρ, hydrogen bond density (HB) and dangling bond (DB) are calculated as follows : where V is the volume i.e., total number of lattice points.
η(x) and φ(x) are discrete fields varying in the range [−N, N]. By construction the partition function is independent of N for N ≥ 8. In practice it is convenient to evaluate this partition function by taking N → ∞, whereupon effective η(x) and φ(x) become continuous fields. We implement this limiting procedure and check if the sum rule is obeyed. In the N → ∞ limit, summation over η and φ is replaced by integrals. The resulting functional integral has the following trivial property : Taking derivatives explicitly in the above gives terms proportional to ν, λ, µ. Since these are being summed at all sites each of these terms can be regrouped in terms of derivatives of ν, λ, and µ as : The µ dependent term in the Eq.(12) has contributions from the four neighboring sites.
Since all sites are being summed over, the µ term in Eq.(11) gets each contribution four times. We notice that this equation is precisely the "sum rule" constraint Eq.(1). This demonstrates that the sum rule in terms of continuous auxiliary fields is automatically true.

III. MOLECULAR MEAN FIELD THEORY
Next, we evaluate the functional integral within the MMF approximation. The partition function's integrand can be envisaged as a product of field dependent phase factors at each site. When we enumerate them site-by-site the phase factors cancel exactly corresponding to physically allowed configurations. Evaluating along this procedure amounts to the standard high temperature or Mayer-like expansion. Instead we attempt an approximate method wherein we first notice that if we relax the constraints Eq.(4) the integrand still peaks for the same configurations that obey Eq.(4) strictly. Hence in the thermodynamic limit, approximating the integrand suitably around the peaking configurations we may reliably estimate the partition function. This reliability can be self consistently established by computing the variance or correlation functions.
The leading contribution to the functional is expected to come from the extremum which maximizes the integrand Z site . Furthermore, since we are choosing to describe fluid phase of the model we seek such spatial configurations which are discrete translational and rotational invariant in auxiliary fields. The Z site over a space-independent field configurationη,φ is given by : It is evident that the maximum for Z site occurs atη =φ = 0 since all fugacities are positive. Z site at the maximum is given by Z o : The extremization with respect toφ is trivially true, while that with respect toη yields, This is a consequence of the sum rule Eq.(1) within the zeroth order approximation. Various densities are calculated from Eq.(10) with the approximate partition function Z o . Furthermore, using Eq. (15) : Eliminating λ in Eqs. (16,17) we obtain : This is the "equation of network" which indeed is a different way of writing the sum rule Eq.(1) in terms of dangling-bond fugacity ν. Z o is given by : Rewriting it as the last term in Eq. (19) shows that it is also inverse of the density of void sites (every water molecule in our model necessarily occupies five sites and for every hydrogen-bond one site is double-counted and hence is compensated in the expression).
From the sum rule it follows 0 ≤ HB ≤ 2ρ. Consequently, ρ here varies between ν 5ν+2 and 1 3 . The upper bound on ρ is indeed the highest possible density in the model while the lower bound is a consequence of MMF approximation, meaning that this description is self consistent only for densities greater than ν 5ν+2 . Without loosing any generality, we choosẽ λ = 1 or λ = e β i.e, measuring all energies in the units of hydrogen-bond energy. Then we make the observation that if temperature ( 1 β ) is always positive, from Eq. (17) we can show that ρ is greater than 1 5 . Furthermore, as β → ∞, from Eqs. (16,17) we see that ρ → 1 3 , HB → 2 3 and DB → 0, showing that every water molecule's arms are all hydrogen-bonded leaving no void sites. This saturation at ρ = 1 3 is verified to be exactly true by explicit construction of such configurations. We find that the highest density is not that of a unique crystal configuration. Instead, there are infinitely many configurations corresponding to different spatial and orientational arrangements of water molecules.
The equation of state in terms of densities is given by : where P o is the pressure and β is written in terms of natural energy units of the model.
By considering the case where H and L arms are neither energetically favored or suppressed i.e.,ν = 0 and ν = 1, the entropy per site is given by, In the limit β → ∞, ρ reaches the maximum value 1 3 and the entropy at the highest density tends to a constant value ln( 3 2 ). The constant 3 2 compares exactly with that of Pauling's estimate for tetrahedral ice i.e., 2 2 6 16 ≡ 1.5 [30] and agrees well with the numerical estimate by Nagle i.e., 1.50685±0.00015 [26]. We note that these results are independent of dimension, hence in two dimensions it also compares with the exact result of Lieb on 'square ice' [21].

A. Fluctuations
Next, we evaluate the functional integral of the partition function by considering small fluctuations around the mean field and approximating them by a Gaussian i.e., one-loop correction [7]. Expanding Z site up to quadratic terms, we obtain : are the reduced fugacities; such that all of them are less than 1.
Inserting the above expression for Z site in Eq. (7) and evaluating the resulting Gaussian integral by Fourier transform in a periodic box the pressure P is given by : where and ∆ = (1−cos(k i )). Thereupon the densities ρ, HB, DB are calculated using Eq.(10).
We compare these results with that of MC simulations in a later section (see sec.IV A on MC).

B. Correlation functions
The position space correlation functions for η(x) and φ(x) fluctuations i.e., G ηη (r) and G φφ (r), to the leading order, are given by : We note that to zeroth order µ ′ ≃ ρ, λ ′ ≃ HB, ν ′ ≃ DB and hence, using Eq. (15) : The asymptotic behavior for large r of G ηη and G φφ correlators can be obtained by pursuing small-k expansion of the integrand and noting that for small k, ∆ ≃ 1 12 i k 2 i . The G ηη correlator for large r is, This shows G ηη has periodic peaks reminiscent of the coordination shell peaks whose amplitudes are exponentially falling off.
The G φφ correlator takes the following asymptotic form for large r, in addition to oscillatory behavior prominent at short distances : We make the observation from Eqs. (28,29) that the correlation functions Eqs. (26,27) of the mean field theory diverge if ρ → 0 i.e., even the local fluctuations about the mean field are very large rendering the approximation invalid. Indeed the theory fails well before that because it violates the sum-rule already at ρ = ν 5ν+2 . This is the reason why our mean field theory fails for the gas phase and thus does not describe any liquid-gas transition. To sum it, in our model MMF approximation shows density saturation and the description is more accurate at higher densities only.
Next, we compute some important physical correlations such as water-water and chargecharge.
The non-zero value of W -field at each site x, y picks only the term proportional to the fugacity µ in Z site expression Eq. Thereby, correlation of water fluctuations to leading order is given by : A similar computation for the q(x) field can be done, wherein the non-zero q value picks terms proportional to the fugacity ν in Z site at either sites (Eq.8). To the leading order, Note that the single-site average over q(x) field is zero. We remark that in this model W (x) being charge neutral gets contribution from the neutral η(x) field, while q(x) from that of the charged φ(x) field.

C. Coulomb interaction
In this section we consider the influence of long range Coulomb interaction potential between the charges H, L given in terms of charge e as : where prime over summation means x = y. In our model we envisage the charges at the tip of water's H or L arms. This can be incorporated in our analysis by using an auxiliary field technique : where the laplacian operator χ(x) = α (χ(x + e α ) − χ(x)) = α ∇ α χ(x) and m is a parameter that regulates the range of interaction. The interaction potential behaves as e −mr r for large distances r, which also reduces to Coulomb interaction when m = 0. If the lattice constant and m are both taken to be zero, then it reduces to exact Coulomb interaction for all r.
By inserting the above in our partition function Eq.(7) all the interactions of water degrees of freedom remain unchanged with the following transformation η → η, φ → φ + χ βe 2 .
The extremum of the new partition function is still atη =φ =χ = 0. The leading zeroth order term remains unchanged; the one-loop correction about the mean field gets additional contributions due to quadratic terms corresponding to φχ and χχ in the Gaussian expansion, given by : The pressure P with Coulomb interactions is, 3 [ln(P ηη (∆)) + ln(P (∆))] where which to the leading order can be written as : The two-point correlations of the fields in momentum space are given by : The • if there is no water, "insert" one provided there are four arms free around the site • if there is water, with equal probability either "delete" water or "rotate" the arms to another possible allowed configuration as dictated by Eq.(5).
The new configuration is accepted with a probability of Boltzmann weight over energy change in compliance with detailed balance condition [10,45]. All energies β,ν,μ are calculated in units of hydrogen-bond energy i.e.,λ = 1. Furthermore, we considered specificallyν = 0 case thereby making the theory essentially parameter-free other than temperature and chemical potential which are varied in steps as per needs of the simulation.
The thermodynamic observables of theμV T ensemble i.e., number of water molecules ρ, number of hydrogen bonds HB, total energy E are obtained as averages from the simulation. A significant equilibriation time is allowed before the production run begins. In The size dependence of the averages is ascertained and an optimal lattice size of 20 sites per side is found to be closely reproducing averages upto fourth decimal relative to bigger lattice sizes.
A 'successive seeding' procedure is employed to adiabatically progress in the parameter space of the model. For example, to carry out simulations at a fixed temperature and a range of chemical potentials an initial run carried out at low densities and equilibriated over long time is used as a 'seed' i.e., input configuration for the successive run. This is repeated to reach higher densities. This procedure accelerates equilibration considerably compared to any random seed configuration. Furthermore, since a range of chemical potentials are explored the step size in the parameter value is appropriately adjusted such that a quench-like situation is avoided. Thereby, the system evolves smoothly in configuration space without any unwanted domains persisting. The end averages remain unchanged when the seeding procedure is carried out in an alternative parameter space, for example, instead of chemical potential temperature can be varied in small steps. This confirms the absence of any possible bias created by our successive seeding procedure in most part of the parameter space (except near first order phase transitions where hysteresis exists).
In order to obtain equation of state using MC we employed fixed temperature simulations, followed by application of Gibbs-Duhem procedure to obtain pressure versus density curves [11,28,29]. In theμV T ensemble temperature is fixed and a range of chemical potentials are explored starting from zero density where pressure can be normalized to the value zero.
Employing the successive seeding procedure we obtain various thermodynamic averages as functions of chemical potential. The Gibbs-Duhem procedure involves integrating the ρ vs.
µ curve so that pressure can be obtained from the following relation : where theμ i corresponds to ρ = 0 and P = 0 andμ f to that of the desired ρ. A set of representative ρ vs.μ curves are shown in Fig.(4). Each such curve has two prominent shoulders where considerable slope change occurs. At very highμ values no more equilibriated configurations could be traced and the density starts shooting up to the saturation value.
We also made preliminary investigation in studying phase transitions in the model. As shown in Fig.(4) the model exhibits discontinuity in ρ(μ) for temperatures β > 2 (in units of hydrogen-bond energy). For instance, at β = 3 the system jumps from a low dense (ρ ∼ 0.025, h ∼ 1.2) to higher density (ρ ∼ 0.16, h ∼ 1.7) where h = HB ρ is average number of hydrogen-bonds per molecule. Hysteresis is also observed when the system progresses first from low-dense to high-dense state upon increasingμ and then retracing the path by decreasingμ. This indicates the presence of first-order phase transition in the region. We interpret this as liquid-gas transition in the model.

A. MMF and MC
One of the important expositions of the MMF theory is the "equation of network", a functional relationship between total water density and HB density. From Eq.(18), at ν = 1, the mean field equation of network is given by, In MC simulation we identify the liquid-gas transition for temperatures β > 2. As shown in Fig.(4), for β = 3.0 atμ ≃ −1.91 water density jumps from 0.025 to 0.16. The jump gap increases with β, so does the higher density state to which system jumps. Within our limited exploration of the phase diagram we identify that water densities ρ 0.16 correspond to liquid phase, only above which MMF theory is self-consistent (Eq.46). Furthermore, the linear relationship between ρ and HB is accurate at high pressures as confirmed by MC simulation. As noted earlier, MMF approximation complies with the sum rule of the system; the same is exactly satisfied by MC simulation at every configurational move and hence over ensemble average. Thus, the agreement between mean field (zeroth order and one loop correction) and MC over the equation of network is inescapable. However, since the mean field description holds forte in thermodynamic limit, the HB density in the above equation corresponds to lowest free energy (or high pressure) states for each water density.
Next, we compare the MMF computation of equation of state with that of MC simulation.
From theory, at ν = 1, The comparison is put forth in Fig.(7). It shows that the locus of high pressure states at each density in MC simulation has the same profile as in MMF theory (both in zeroth order and one-loop'). Further, the qualitative agreement is good only in the high density regime where the equation of network too is accurate. A quantitative comparison of equation of state between MMF and MC is unreliable because MMF pressure absurdly vanishes at ρ = 1 7 whereas, physically, the pressure is zero in this model only at ρ = 0. As discussed earlier, MMF approximation fails for small densities. Therefore, a consistent normalization between various schemes of calculation is not present. Thus, the qualitative picture obtained from MMF calculation is only indicative, nevertheless consistent with MC results.
The spatial correlation functions are also computed from the MC simulation and compared with that of analytical expressions obtained within the mean field approximation.
Firstly, we note that the underlying hyper-cubic lattice dominates all correlations especially at short distances. Then, to extract the rotational invariant part of any function f (x) we implement the following projection procedure : where | x| = i x 2 i and Θ is Heaviside step function defined as Θ(x − a) = 1 for x ≥ a and 0 for x < a. The important correlations in the model are < W (0)W (r) > and < q(0)q(r) >.
The k-space integrals in Eqs. (35,36) are computed numerically and the MMF correlation functions are compared with those of MC simulation. The positions of coordination peaks in both are in agreement. In addition, charge correlations show similar exponential fall-offs asymptotically. Fluctuations are seen to increase as we approach saturation density in both cases, however they are higher in simulation.

V. MMF, EXPERIMENTS AND MD
MMF theory predicts correlation length of charge fluctuations i.e., in qq correlation function Eqs. (33,36). These fluctuations could give rise to long distance correlations in other charged fields in the system such as in dipole, as seen in MD simulations (wherein two correlation lengths of order 5.2Å and 24Å have been observed, former being 10 times larger than the latter in strength) [17]. To relate to MMF we make the following observations.
Water molecule on average participates in 3.6 hydrogen-bonds in ambient conditions [14,23,24] i.e., DB ρ = 0.4. Using this value in m 2 of Eq.(33) the correlation length turns out to be approximately 2.02× lattice constant in MMF. We next estimate the lattice constant in physical units by comparing the position of highest coordination peak in W (0)W (r) correlation obtained from MMF to that of experimental data. This procedure gives a lattice constant equal to 2.8

1.75Å
. Hence in physical units MMF predicts a correlation length of about 3.2Å for liquid water. It should be noted that these predictions are not robust as the coefficients such as 5 24 in Eq.(33) might vary with topology of the underlying lattice. We also pursued a preliminary comparison of the equation of network from experiments and MD simulations. The density of hydrogen bonds is indirectly probed and inferred in various experiments [2,6,18,32,38] and also computed in MD and MC simulations under varying external conditions [5,16,23,24]. The latter works used different definitions for hydrogen-bond computation such as energy-based, geometry-based or hybrid. The data is put in the perspective by converting the mass densities to number densities using known radius of a water molecule. As shown in Fig.(6), in the region of high molecular density i.e., corresponding to liquid water we find that HB density and water density are linearly related to each other. A linear fit function is used for HB vs. ρ curve and compared with that of the MMF equation Eq. (18). We find that the dangling bond fugacity ν is in the range (0.06, 0.18), implying that the corresponding energyν is positive and large compared to thermal energy i.e., dangling bonds are highly disfavored in liquid water.
Another useful quantity namely fraction of water molecules with i hydrogen bonds can also be calculated. Consider a water molecule in a configuration where its i-number of arms are hydrogen-bonded to neighboring molecules and its other (4 − i) arms remain dangling type. A weight can be associated with each such configuration defined in terms of appropriate site fields and summed over all possible orientations of the molecule. We denote this weight averaged with respect to the full partition function for each i as p i . For instance, the averaged weight assigned to a molecule which is hydrogen-bonded to only two other molecules is given by : Thereupon, f i which is fraction of i-bonded molecules can be calculated from the relation, Note that the above expression is obtained to zeroth order approximation within the model.
There exist one-loop corrections to it that can be calculated from the MMF theory, but numerically they are small. These distributions agree very well with MD simulations [5].
Further, probabilities for various cluster configurations (n-mers) can also be calculated within the MMF theory along the same lines as above. Liquid water is known to form molecular clusters such as trimers, tetramers, pentamers [33].

VI. DISCUSSION
We analyzed the statistical mechanics of molecules with attractive bonding interactions.
We argued that these systems typically have sum rules built-in which have to be respected in the analysis. In addition, there are other constraints, such as those shown in Fig.(2) for water, all of which are captured by defining fields such as b(x) and q(x) as in Eq.(4). The analysis yields ultimately a functional integral in terms of auxiliary fields η(x), φ(x) which are conjugates of the discrete fields b(x) and q(x) respectively. A simple-minded mean field analysis of the auxiliary fields is shown to reproduce all major thermodynamic properties of the system, which are verified to be reasonably accurate in liquid phase and near saturation density. Principally, "equation of network" which is another manifestation of the "sum rule" plays an important role. We showed that all these consequences are known to be qualitatively true in experiments and MD simulations as well.
The MMF approximation is shown to be the leading order in a formal Feynman loop expansion of the functional integral with necessarily small parameters called reduced fugacities.
Correlation functions and their properties have been computed and shown to qualitatively agree with MC simulation. It is however observed that many of the details in correlation functions do depend upon the underlying lattice, but the analytic structure is amenable to interpretation in the continuum as well. Pressure and correlation lengths have large con-tribution from the fluctuations of DB and ρ in the form DB ρ . We note that any long range interactions such as Coulombic have very little quantitative effect on our results.
The MMF theory is seen to be inconsistent at low densities, consequently fails to describe the gas phase of the model and the corresponding liquid-gas phase transition. In the explicit model considered here we do not have a crystalline solid phase, but we may have a glassy phase. Namely, at zero temperature the lowest energy configuration is infinitely degenerate.
In this region MMF theory does not show any liquid-glass first order phase transition, however it does indicate a second order transition as seen from q(x) correlator. In Eq. (33) for β → ∞: ρ → 1 3 , DB → 0 implying m 2 → 0. In MC simulation we did observe dynamical slowing down as we reach saturation density, perhaps indicative of a phase transition.
The MMF technique can be applied to associating liquids which may constitute variety of molecules and hence with variety of sum rules. Therefore, we expect a variety of dangling bond fluctuations contributing to pressure and correlation lengths. Indeed, there may even be variety of phase transitions. Essentially, the MMF technique is amenable to analysis of all such systems.

VII. ACKNOWLEDGMENTS
We acknowledge discussions with Prof. Gautam Menon on MC simulation data analysis.