# Towards an effective potential for the monomer, dimer, hexamer, solid and liquid forms of hydrogen fluoride

## I Introduction

Hydrogen fluoride (HF) is one of the simplest molecules capable of forming hydrogen bonds. Despite the simplicity of the compound, the theoretical description of its behavior, especially in the liquid phase, is still far from being fully satisfactory. Much relevant theoretical work has concentrated on the problem of determining a potential model suitable for computer simulations. For this purpose a model is needed that reproduces correctly the main features of the real interaction potential and which is simple enough to be computed efficiently.

In the last two decades several empirical potentials have been developed for Molecular Dynamics (MD) or Monte Carlo (MC) simulations of liquid HF.  Of these models, only the three-site ones  are able to reproduce correctly the dipole and quadrupole moments of the monomer. These models represent the charge distribution of each monomer by fractional charges placed at three sites on the molecular axis, two on the F and H nuclei and the third one at an appropriate position X along the F-H bond. A first three-site model (called HF3) was developed by Klein and McDonald  by fitting ab initio results for the potential energy surface of the (HF)<sub>2</sub> dimer. Cournoyer and Jorgensen  proposed a second three-site model (called HFC), with a simplified non-Coulombic part consisting of a single Lennard-Jones interaction between the fluorines. The parameters of the model were fitted directly to experimental thermodynamic data for the liquid, while simultaneously providing a reasonable equilibrium geometry of the dimer. Recently, Jedlovszky and Vallauri  presented two further models, hereafter referred to as HF-JV1 and HF-JV2, respectively. HF-JV1 is a variant of HFC, with charges reproducing both the monomer dipole and quadrupole and with an accurate treatment of the long range part of the Coulombic interactions, neglected in the original work on HFC.  The HF-JV2 model includes molecular polarizability, by adding an induced point dipole moment at each F site, while keeping the charge distribution of HF-JV1 unchanged. The scalar polarizability of the molecules was set equal to its experimental value, while the two parameters of the Lennard-Jones interaction between F atoms were fitted to the experimental values of the liquid density and internal energy.

Unfortunately, none of the available models is able to reproduce, in a fully satisfactory way, both the thermodynamics and the structure of liquid and gaseous HF. In the search for a new potential, the essential ingredients can be identified by reviewing the known properties of the gas, liquid and solid phases of HF. The simplest associated form of HF is the gas phase dimer, (HF)<sub>2</sub>, whose structure and rovibrational spectrum were first characterized by the microwave spectra of Dyke, Howard and Klemperer.  The equilibrium configuration of the isolated dimer F-H$`{}_{}{}^{}\mathrm{}`$F-H<sup>′′</sup> is planar but bent, as a consequence of a competition between dipolar and quadrupolar electrostatic interactions  (whenever ambiguity is possible, H denotes the hydrogen atom involved in the hydrogen bond, while H<sup>′′</sup> is the other one). The atoms F-H$`{}_{}{}^{}\mathrm{}`$F form a nearly linear arrangement, with H placed slightly off the FF axis, and with a distance between the hydrogen-bonded fluorines $`r_{\mathrm{FF}}`$ 2.7 Å. The second hydrogen atom forms an angle $`\widehat{\mathrm{FFH}^{\prime \prime }}115^{}`$ and the H–F bonds are slightly stretched with respect to the monomer. 

The (HF)<sub>2</sub> bent dimer appears to be the basic structural motif of all the best known associated forms of HF, namely the gas phase (HF)<sub>6</sub> cyclic hexamer,  the low temperature crystal  and even the liquid.  Pairs of adjacent HF units in the hexamer exhibit a bent arrangement similar to that of the dimer. The same structure is found in crystalline deuterium fluoride (DF) which is made up of infinite zig-zag chains of DF units.  Interchain F–F distances ($`3.2`$ Å) are much larger than intrachain F–F distances ($`2.5`$ Å), an indication that interchain interactions are weaker than the interactions between adjacent HF units within the chains. The very small entropy difference between the solid and the liquid  suggests that the liquid is largely associated. Best fit analysis of dielectric constant data  and of Raman and IR spectra  support this conclusion and indicate that acyclic zig-zag chains are abundant in the liquid. Finally, the neutron scattering measurements  of the radial correlation function of liquid DF yield H–F and F–F average neighbor distances consistent with a most probable local structure similar to that of the other forms of HF.

Since all associated forms of HF share the (HF)<sub>2</sub> dimer as a common structural unit, any satisfactory potential model should reproduce the main features of the (HF)<sub>2</sub> potential energy surface. Ab initio quantum mechanical calculations have contributed significantly to clarify the interactions present in the isolated dimer (HF)$`_2,`$ by determining accurately its equilibrium properties  and by mapping out the entire potential energy surface.  Much of this work has been directed to obtain analytical models for the potential surface, usually by fitting the numerical results to a chosen functional form.  The analytical expressions of these models, which have been developed mainly to study the rovibrational states of (HF)<sub>2</sub>,  are usually very complex and therefore not suitable for simple MD or MC simulations.

Fortunately, as already mentioned, there is evidence that the main features of the real molecular interactions may be approximately reproduced by much simpler empirical models with effective potentials. In fact, it has been long known  that the bent structure of the (HF)<sub>2</sub> dimer is dominated by the classical two-body electrostatic interactions between the permanent multipole moments of the monomer (mainly the dipolar and quadrupolar terms). The incorporation of many-body polarization effects gives rather small refinements on the predicted equilibrium angles,  a finding consistent with the experimental observation that the dipole moment of the dimer is only weakly enhanced relative to that of the isolated monomer. 

The effects of the environment on the HF molecules do not seem to be large. However, a comparison of the structures in the sequence “monomer $``$ dimer $``$ hexamer $``$ liquid $``$ solid” (as shown by Tables II and IV below) evidences a tendency in which a larger degree of association implies a small increase of the H–F bond distance, together with a much larger decrease of the F–F distance (from 2.8 Å in the gas phase dimer to 2.49 Å in the solid). Unfortunately, as stressed by Röthlisberger and Parrinello,  the available three-site models  refer to rigid molecules and, consequently, they cannot reproduce the relaxation of the interatomic distances in going from the gas phase dimer to the condensed phases.

The previous observations suggest a relatively clear picture of the ingredients required for a satisfactory potential model for HF. First, a charge distribution is required which approximately reproduces the first few multipolar moments of the HF molecule. Second, the H–F bond cannot be considered as rigid and intra-molecular interactions must be included to allow for the observed variations of the H–F and F–F distances in the various aggregation forms of HF. Finally, further atom–atom interactions must be introduced to model the remaining non-Coulombic intermolecular forces. Having all this in mind, we have tried to construct a new three-site model by simultaneously using data on the gas, liquid and solid phases. For this purpose: (a) the molecules are not rigid, but the H–F bond length can vary, and, (b) the parameters are fitted to theoretical and experimental data, including the ab initio structure of the HF dimer,  the room temperature density of liquid DF  and the experimental structure of solid DF at 4 K. 

For the structure of the dimer we have used the ab initio (HF)<sub>2</sub> potential energy surface developed by Bunker, Jensen, Karpfen, Kofranek and Lishka (BJKKL).  The BJKKL calculations, in which the potential energy of the (HF)<sub>2</sub> complex has been determined for over 1000 different configurations, represent the most complete and accurate scan of the energy surface to date. These quantum mechanical results are in excellent agreement with the experiments on gas phase (HF)<sub>2</sub>.

The decision to add to the fit some data on solid and liquid HF (more precisely, DF) has been taken because of the partial failure of our preliminary models fitted only to the ab initio surface. These models reproduced well the zig-zag (HF) chains characteristic of the crystal, but gave totally wrong inter-chain distances and, furthermore, did not agree with the experimental density of the liquid.  This behavior was to be expected. In fact, the (HF)<sub>2</sub> potential surface  only accounts for the basic HF–HF interactions within a chain, and obviously exclude the weak long-range interactions responsible for the distances between different (HF) chains. Since the density of the liquid is affected by interactions between distant pairs of HF molecules in relative orientations which are not sampled in the solid, it is also understandable that only by fine-tuning the long range potential it has been possible to reproduce the experimental density of the liquid. Finally, it must be mentioned that the addition of data on solid and liquid HF to the fit gave only a small deterioration of the agreement with the ab initio data on the dimer. This fact confirms that solid and liquid data add information on regions of the potential surface that are not sampled by the dimer.

## II Methods and Calculations

### A Potential model

The potential model is represented by intra- and inter-molecular parts:

$$V_{\text{HF}}^{\mathrm{intra}}(r_{\text{HF}})=D_e\left\{1\mathrm{exp}[\alpha (r_{\text{HF}}r_e)]\right\}^2,$$
(1)   
$$V_{AB}^{\mathrm{inter}}=V_{AB}^{\mathrm{Coul}}+V_{AB}^{\mathrm{non}\mathrm{Coul}},$$
(2)   
$$V_{AB}^{\mathrm{Coul}}=\underset{iA}{}\underset{jB}{}\frac{q_iq_j}{r_{ij}},$$
(3)   
$$V_{AB}^{\mathrm{non}\mathrm{Coul}}=A_{\mathrm{FF}}\mathrm{exp}(B_{\mathrm{FF}}r_{\mathrm{FF}})C_{\mathrm{FF}}r_{\mathrm{FF}}^6.$$
(4)   
BJKKL  fitted the intra-molecular part of their own ab initio surface with a Morse potential, eq. (1). Here $`r_e`$ represents the equilibrium H-F distance, $`D_e`$ the dissociation energy and $`\alpha `$ is an effective range parameter. Because of their simplicity and accuracy, the BJKKL functional form and parameter values are adopted in this paper.

The Coulombic interactions between molecules $`A`$ and $`B`$ are modeled through three point charges for each HF monomer, two at the nuclear positions $`𝑹_\mathrm{H}`$ and $`𝑹_\mathrm{F}`$, and the third at a position $`𝑹_\mathrm{X}`$ along the H–F bond. In eq. (3) $`q_i`$ and $`q_j`$ are the fractional charges on the $`i`$th site of molecule $`A`$ and $`j`$th site of molecule $`B,`$ respectively; $`r_{ij}`$ is the distance between these sites. The motion of the site X is constrained so that it remains at the same relative position along the bond,

$$𝑹_\mathrm{X}=\beta 𝑹_\mathrm{F}+(1\beta )𝑹_\mathrm{H},$$
(5)   
where $`\beta `$ is an adjustable parameter between 0 and 1. The charges are $`+q`$ at both the H and F nuclei, and $`2q`$ at the third site to preserve neutrality. This three-site charge model is related to that of Refs. . By allowing for changes in the H–F bond length, the present model effectively accounts for a part of the polarization effects. For solid and liquid DF the Ewald’s method  has been used to ensure complete convergence of the Coulombic interactions (which have an infinite range).

The remaining non-Coulombic part of the inter-molecular potential is represented in a simplified way, using only a Buckingham “$`\mathrm{exp}6`$” atom-atom interaction between the fluorines, eq. (4). This term is meant to represent the interactions between the electronic clouds around two far away atoms. Since the hydrogens in HF are essentially bare nuclei, it makes good physical sense to avoid atom-atom interactions involving them. As a matter of fact, no improvement in the quality of the fit is found by adding similar H–H and H–F non-Coulombic interactions.

A rather well defined hierarchy of interactions may be identified in the chosen potential model. The length of the HF monomer is solely determined by the intra-molecular potential. The structure of the dimer is also influenced by the position of the charge site X and by the F–F equilibrium distance, i.e. by the position of the minimum of the $`\mathrm{exp}6`$ interaction between the fluorines. Finally, the charge and the remaining properties of the $`\mathrm{exp}6`$ model, mainly the strength of the long range attractive term $`C_{\mathrm{FF}}r_{\mathrm{FF}}^6`$, affect the structure and density of solid and liquid HF. The presence of this hierarchy implies that small changes in the charge and in the long range attraction can be compensated by the remaining parameters to maintain the correct monomer and dimer structures.

### B Potential optimization

The potential model contains five adjustable parameters, $`\beta `$, $`q`$, $`A_{\mathrm{FF}}`$, $`B_{\mathrm{FF}}`$ and $`C_{\mathrm{FF}}`$, and three parameters fixed at the ab initio values,  $`r_e`$, $`\alpha `$ and $`D_e`$.

In a first series of attempts, we fitted the present model (as well as other preliminary ones) only to the ab initio potential energy of the dimer. The $`\chi ^2`$ deviation between the model and the ab initio surface was computed, for each given combination of parameter values, with the same relative weights used by BJKKL.  The $`\chi ^2`$ was minimized by searching the parameter space with Nelder-Mead simplex method.  The resulting potential was then tested by computing some properties of the other associated forms of HF. In particular, the liquid phase was studied by isothermal-isobaric MD simulations, as described in the next Subsection. As discussed in the introduction, these preliminary models fitted only to the ab initio surface gave unsatisfactory results, so that it was decided to add more data to the fit.

For this purpose, the equilibrium geometries at $`T=0`$ K of the HF dimer, hexamer and crystal were determined as a function of the potential parameters by minimizing, with the WMIN program,  the total potential energy with respect to the structural parameters. The deviations of the calculated geometries from the ab initio dimer structure  and from the experimental DF crystal structure  at 4 K are then added to $`\chi ^2`$, with weights subjectively chosen to make the contributions from the surface, the dimer and the crystal roughly equal.

Since the new set of parameters, although more satisfactory, still did not give the correct density of the liquid, it become necessary to tune the long range interactions. In a further set of minimization runs, a range of $`q`$ and $`C_{\mathrm{FF}}`$ values was searched, again with Nelder-Mead method. The three remaining parameters $`\beta `$, $`A_{\mathrm{FF}}`$ and $`B_{\mathrm{FF}}`$ where determined as a function of $`q`$ and $`C_{\mathrm{FF}}`$ by fitting the ab initio and crystal data. Each complete set of five parameters was then used in a short MD simulation to determine the equilibrium density of the liquid at 293 K and to add to the $`\chi ^2`$ the deviation from the relevant experimental value.  This method, which involves two nested fit procedures, was found to be reasonably efficient and converged to a minimum in about fifty cycles.

The main problem encountered in the fit was the noise in the computed liquid density due to insufficient MD equilibration with our computer time constraints. To reduce this noise, the MD equilibration was done in parallel with the potential optimization, by accepting after successful $`\chi ^2`$ minimization cycles the final configuration of the MD run as the initial configuration for the next run. With this strategy the current MD configuration was always the one with the best potential parameters so far. Since the parameters change rather slowly, the simulated system tended to remain close to equilibrium. No structural data on the liquid, beside the density, has been included in the fit. This rather drastic choice avoids the repeated calculation of equilibrated radial correlation functions, which would have required even longer MD runs then those needed for equilibrated densities.

As a technical detail, it must be noticed that no attempt was done to embed potential surface calculation, MD simulation, dimer and crystal energy minimizations and the two nested Nelder-Mead procedures into a single monolithic program, which would have been unmanageably complex. Separate programs, calling each other as distinct processes at the operating system level,  were used instead. The two Nelder-Mead procedures, in particular, were actually a single program invoking a second copy of itself. No special changes were required for the surface, MD and energy minimization programs.

The optimal set of parameter values is shown in Tab. I and represents a compromise among the best results that can obtained separately for the dimer, the crystal and the liquid. As usual, no special physical significance should be attributed to the potential parameters. In fact, because of the possible compensation among different terms in the potential model, alternative slightly different sets of parameters might have been used. The model may be simply regarded as a tool to reproduce the observed data and predict new results.

### C Molecular Dynamics

The MD calculations employed 500 deuterium fluoride (DF) molecules, in a cube with periodic boundary conditions, and, using Andersen’s isothermal-isobaric (NPT) method  simulated a liquid sample in contact with a heat bath and subject to a hydrostatic pressure of 1 atm ($`10^4`$ GPa). The simulated liquid was obtained by melting and equilibrating at 293 K an initial crystalline configuration. The behavior of the system as a function of temperature was then determined by raising or lowering the bath temperature in steps of $`10`$ K. Each temperature was maintained for at least 5 ps for equilibration and further 5 ps for analysis. The equations of motion were integrated using the velocity Verlet algorithm,  with a time step of 0.25 fs. As previously described, each DF molecule carried three point charges, at D, F and X sites. Ciccotti, Ferrario and Ryckaert method for linear constraints  has been used to maintain each massless X charge at a fixed fraction $`\beta `$ of the DF bond.

## III Results

The most important properties calculated with the present potential model for the monomer, dimer, planar (HF)<sub>n</sub> rings, hexamer, crystal and liquid forms of HF (or DF) are compared in Tables II, III and IV with the available experimental and ab initio data. As described in section II C, the properties of the liquid have been determined through MD calculations, while the equilibrium geometries at 0 K of the other forms of HF are found by minimizing the potential energy.

### A Monomer, dimer and cyclic polymers of HF

The excellent results (Tab. II) for the equilibrium bond length, dissociation energy and spectroscopic parameters of the monomer, which all depend only on the parameters fixed to the BJKKL values,  indicate that the Morse model accurately reproduces the main features of the true intra-molecular potential. The spectroscopic parameters $`\nu _e`$ (harmonic frequency) and $`x_e`$ (anharmonicity constant) for the energy levels of a Morse oscillator, $`E_n=h\nu _e[(n+\frac{1}{2})x_e(n+\frac{1}{2})^2]`$, are directly obtained from $`D_e`$ and $`\alpha `$ through $`\nu _e=\alpha \sqrt{D_e/\mu }`$ and $`x_e=h\nu _e/4D_e`$, where $`\mu `$ is the reduced HF mass. 

The multipole moments computed from our three-charges model (taking the molecular center of mass as origin) are very close to the experimental and ab initio moments of the monomer.  The dipole and quadrupole values follow the trend of the ab initio results and therefore are slightly underestimated with respect to the experimental data. The octupole and hexadecupole moments are also in reasonable agreement with the ab initio calculations. This overall agreement is an indication of a good match between the electrostatic interactions in real HF and in the model.

The computed minimum energy structure of the dimer (Tab. II) compares well with the experimental and ab initio data.  The experimental F-F distance is excellently reproduced, as well as both the angles $`\widehat{\mathrm{H}^{}\mathrm{FF}}`$ and $`\widehat{\mathrm{FFH}^{\prime \prime }}`$ of the bent equilibrium configuration. Moreover, the increased length of the HF molecule in going from the monomer to the dimer, and the slight length difference between the two HF intramolecular bonds,  are well predicted. Since these length changes were the primary reason for allowing non-rigid molecules, such a behavior must be considered very satisfactory. Unfortunately, in spite of the excellent dimer geometry, the dimerization energy is clearly underestimated. This drawback was not completely unexpected. In fact, in real HF and in quantum mechanical models dimerization is accompanied by hydrogen bonding, which is not explicitly incorporated in the present classical treatment.

As a further test of the potential, we have computed the equilibrium geometry of the planar (HF)<sub>n</sub> rings, with $`n=2,3\mathrm{}8`$. As shown in Tab. III, the structural parameter computed for planar rings, with $`C_{nh}`$ symmetry, are in good agreement with the available ab initio results.  Though the model systematically underestimates the ab initio binding energies, it nevertheless reproduces correctly the relative stability of the different structures. The smallest ring, which is the cyclic dimer, has a binding energy of $`\mathrm{\Delta }E=3.18`$ kcal/mole, and is thus substantially less stable than the bent dimer (Tab. II).

For cyclic (HF)<sub>3</sub> the binding energy for each hydrogen bond, $`\mathrm{\Delta }E/3`$, is slightly less that the binding energy of the bent dimer. The bond stabilization energy, $`\mathrm{\Delta }E/n`$, increases for larger rings, up to the hexamer, and then decreases again. The particularly favorable stability of the (HF)<sub>n</sub> rings with $`n6`$ is readily understood by noticing that for these rings the $`\widehat{\mathrm{H}^{}\mathrm{FF}}`$ and $`\widehat{\mathrm{FFH}^{\prime \prime }}`$ angles (Tab. III) are close to those of the bent dimer (Tab. II). Planar (HF)<sub>n</sub> rings must satisfy the geometric constraint $`\widehat{\mathrm{FFH}^{\prime \prime }}\widehat{\mathrm{H}^{}\mathrm{FF}}=\alpha _n`$, where $`\alpha _n=180^{}360^{}/n`$ is the inner angle of the $`n`$-sides regular polygon. The hexamer, for which $`\alpha _n=120^{}`$, can be obtained by joining essentially undeformed bent dimers, and is the most stable structure.

The hexamer can be stabilized even further by allowing for non planar structures. We find that the stablest structure has a non-symmetric “chair” shape, with average bond lengths and angles (Tab. II) which compare well with the available experiments  and which are almost identical to those found in the dimer. A “boat” structure slightly above in energy, at $`4.96`$ kcal/mole, is also found, with lengths and angles close to those of the “chair” structure.

### B Crystal and liquid DF

Low temperature crystals of DF are orthorhombic, space group $`Cmc2_1`$ ($`C_v^{12}`$), with four molecules per unit cell on the $`\sigma _v`$ plane sites.  The minimum energy structure computed for the DF crystal (see Tab. IV) is close to the experimental structure at low temperature.  The discrepancies in the lengths of the cell axes partially compensate each other, yielding a density only slightly smaller than the experimental one. The increased H-F intra-molecular lengths with respect to both monomer and dimer, and the F-F distance smaller than in the dimer, are well reproduced.

In Fig. 1 the densities predicted by the present model over the whole temperature range of liquid HF (which, at atmospheric pressure, goes from the freezing point at $`83`$ C up to the boiling point at $`19.75`$ C) are compared with the experimental data. The experimental densities show a nearly linear dependence on $`T.`$ The data are from two sources, covering different temperature ranges.  The slight vertical shift between the two data sets is due to the experimental uncertainties. It is very satisfactory to note that, although our potential model has been fitted only to a single density at $`293`$ K, the straight line corresponding to its predictions is somewhat vertically shifted, but has essentially the same slope as the experimental density. This slope is not reproduced by the HFC, HF-JV1 and HF-JV2 models,  since their density straight lines (obtained by joining the corresponding results, which, unfortunately, are available only at $`203`$ and $`273`$ K) intersect the experimental curve.

Another pleasant feature of the present model, also shown in Fig. 1, is that at $`303`$ K and above the average density continued to decrease over the whole MD analysis period. At each $`T<303`$ K, the density oscillated around an equilibrium value. In our opinion, this behavior indicates that the simulated liquid boils at some point in the range $`293÷303`$ K, in good agreement with the experimental normal boiling point of HF ($`T_b=292.9`$ K).

With regard to the internal energy $`U,`$ our MD results exhibit a nearly linear dependence on $`T.`$ However, their absolute values (Tab. IV) are underestimated with respect to the available experimental data  (a similar drawback is also present in the HF3 model ).

Fig. 2 reports the partial pair radial correlation functions $`g_{ij}`$(r) computed for the liquid at 293 K, together with the real space function $`d(r)4\pi \rho _mr\left[G(r)1\right]`$. Here $`\rho _m`$ is the molecular number density, and $`G(r)`$ is a composite (or total) pair correlation function, obtained by adding the three partial pair correlation functions with the appropriate nuclear weights  and then by convoluting the sum with the same experimental resolution function used in the neutron scattering experiments.  The resulting theoretical prediction for $`d(r)`$ is compared in Fig. 2, with the corresponding neutron diffraction data for liquid DF at 293 K, which are the only available real space data. 

The first peak at $`r0.95`$ Å is due to the intra-molecular H–F bond distance and here the agreement between simulation and experiment is excellent. The second and third peak of the experimental $`d(r)`$ occur at $`r1.6`$ Å and $`r2.55`$ Å, and correspond to the hydrogen-bond $`r_{\mathrm{HF}}`$ inter-molecular distance and to the $`r_{\mathrm{FF}}`$ separation, respectively. Unfortunately, the model (with the present choice of parameters) fails to reproduce these peaks and the complex liquid structure at longer distances. The reason of such a shortcoming may be found by comparing the partial pair correlation functions $`g_{ij}(r)`$ of the present model Fig. (2) with the analogous ab initio MD results of Röthlisberger and Parrinello.  The position of the first inter-molecular peak of all our $`g_{ij}(r)`$ is shifted toward larger $`r`$ values (a similar trend occurs in the polarizable HF-JV2 model ). In addition, the height of the hydrogen-bond peak of $`g_{\mathrm{HF}}(r)`$ is nearly half the correct value, and the height of the first peak of $`g_{\mathrm{HH}}(r)`$ is also underestimated.

It is to be recalled that the HF3 potential  reproduces the three principal peaks of the experimental data  and gives the best performance for $`d(r)`$ among the available models. This quite good agreement was obtained by modeling the hydrogen-bond interaction with a Morse term.  To our knowledge, no other empirical model takes explicitly into account the hydrogen bond.

Before concluding this Section, it may be useful to summarize the performances of the best available models for liquid HF. The HF3 potential gives, as already mentioned, quite good agreement with neutron diffraction structural data at the normal boiling point (293 K),  but systematically fails to reproduce thermodynamics:  the predicted internal energies are largely underestimated (i.e., their absolute values are too small) and the pressures in MD calculations at constant volume are several kilobars too high, indicating that the model system is less strongly bound than real HF. In comparison with HF3, the HFC model  yields a slightly better thermodynamics but a slightly worse structure for the liquid, with charges corresponding to a dipole moment enhanced with respect to the monomer. Then, the predictions of HF-JV1 for the liquid phase  are not very different from the HFC ones, but HF-JV1 fails completely to reproduce the properties of the isolated dimer. Finally, whereas HF-JV1 works reasonably well only at room temperature, the polarizable model potential HF-JV2 represents a true improvement  as regards the predicted density at low temperature, i.e. at $`203`$ K. However, problems with the pair correlation functions are encountered also for HF-JV2.

## IV Conclusions

We have shown that a simple potential model suitable for MD or MC simulations can reproduce, quantitatively or semiquantitatively, many physical properties of the hydrogen fluoride over the whole set of its solid, liquid and gaseous associated forms. To cover such a wide range of environmental conditions, it has been necessary to consider a molecular model with variable bond length. The present investigation confirms the plausibility of assuming the dimer as the basic structural unit, but also stresses the need of fitting the potential to a set of experimental and theoretical data which includes information on the condensed phases.

In the liquid, the correct trend of the density is very encouraging, but the underestimated energies and the problems with the radial distribution functions indicate that some physical effect is still misrepresented by the model. Unfortunately, it is not likely that adding data on the liquid structure could improve the results with the current type of model (without an hydrogen bond term). In fact, our experience with the fit shows that those few parameter sets which give better liquid structures were incompatible with the gas and crystal data. We believe that the essential shortcoming of the model is the neglect of any explicit representation for the hydrogen bond, which cannot be reduced to purely electrostatic multipole interactions. Such an explicit modeling of the hydrogen bond is also lacking in most other three-site potentials for HF, with the exception of the HF3 model.  Although the hydrogen bonding has a quantum mechanical origin, an approximate classical treatment is however possible and may have significant consequences, as seen from the rather good structural results for the HF3 potential.  The inclusion of potential terms representing the hydrogen bond interactions appears therefore the next necessary step for more accurate HF models. More data on the liquid, including at least the radial correlation function, need to be incorporated in the fit.

In conclusion, our model cannot be considered as a definitive solution of the problem, but can be seen as a significant step towards a really satisfactory potential for MD or MC simulations. It has the merit of pointing out the importance of a variable molecular length and of the hydrogen bond. Moreover, it shows that the strategy of a simultaneous fit to data covering all the associated forms of HF can be successful and should be considered as the appropriate way to fully accomplish the difficult task of finding a potential model for such a strongly associating system.

###### Acknowledgements.

Work done with funds from MURST (Ministero dell’Università e della Ricerca Scientifica e Tecnologica) through the INFM (Istituto Nazionale di Fisica della Materia), from CNR and from the University of Bologna (“Finanziamento speciale alle strutture”). We thank Bunker, Jensen, Karpfen, Kofranek, and Lishka for providing their ab initio data.