Effective Force Field for Liquid Hydrogen Fluoride from Ab Initio

J. Phys. Chem. B 2005, 109, 6573-6586

6573

Effective Force Field for Liquid Hydrogen Fluoride from Ab Initio Molecular Dynamics Simulation Using the Force-Matching Method† Sergei Izvekov and Gregory A. Voth* Department of Chemistry and Center for Biophysical Modeling and Simulation, UniVersity of Utah, 315 South 1400 East Room 2020, Salt Lake City, Utah 84112-0850 ReceiVed: September 23, 2004; In Final Form: NoVember 23, 2004

A recently developed force-matching method for obtaining effective force fields for condensed matter systems from ab initio molecular dynamics (MD) simulations has been applied to fit a simple nonpolarizable two-site pairwise force field for liquid hydrogen fluoride. The ab initio MD in this case was a Car-Parrinello (CP) MD simulation of 64 HF molecules at nearly ambient conditions within the Becke-Lee-Yang-Parr approximation to the electronic density functional theory. The force-matching procedure included a fit of short-ranged nonbonded forces, bonded forces, and atomic partial charges. The performance of the forcematch potential was examined for the gas-phase dimer and for the liquid phase at various temperatures. The model was able to reproduce correctly the bent structure and energetics of the gas-phase dimer, while the results for the structural properties, self-diffusion, vibrational spectra, density, and thermodynamic properties of liquid HF were compared to both experiment and the CP MD simulation. The force-matching model performs well in reproducing nearly all of the liquid properties as well as the aggregation behavior at different temperatures. The model is computationally cheap and compares favorably to many more computationally expensive potential energy functions for liquid HF.

I. Introduction Hydrogen fluoride is the simplest liquid exhibiting the ability to form strong hydrogen bonds, thus making it a challenging target for molecular modeling. Moreover, the development of reliable empirical potentials suitable to simulate HF in condensed phases has been hampered by a scarcity of available experimental data on the liquid properties of HF because of its extremely high reactivity. The difficulties in designing empirical potentials to describe the properties of HF in condensed phases are no surprise given the long history of potential development for water, which is a much less associative hydrogen-bonded liquid. Nevertheless, in the past decades several empirical potentials have been developed to simulate HF in different phases. The earliest models contained two interaction sites.1-3 However, the inability of these models to describe even qualitatively4 the structural properties of liquid HF led to a widespread belief that only models with additional interaction sites can be successful in a realistic representation of properties of HF in different phases. This situation has been commonly rationalized by the importance of the HF quadrupole moment, which might be responsible for the formation of zigzag chains in the solid and liquid5-9 or cyclic structures in the gas-phase clusters.10 In essence, a strategy to develop a better empirical potential for HF in condensed phases has followed closely the advances in the modeling of water. Water models with additional interaction sites such as TIP3P, TIP4P, and TIP5P,11-13 for example, have resulted in an improvement of the simulated water properties. The first three-site model for HF was proposed †

Part of the special issue “David Chandler Festschrift”. * Author to whom correspondence should be addressed. Phone: (801) 581-7272. Fax: (801) 581-4353. E-mail: [email protected].

by Klein and McDonald.1,14 The model was parametrized by fitting to the potential energy surface (PES) of the HF dimer. Subsequently, several alternative three-site models have been proposed15-20 that differ by the chosen analytical form of the interactions and the way the parametrization was accomplished. The models that appear to be the most successful in describing properties of liquid HF were fitted to reproduce the experimentally derived thermodynamic data over a wider range of temperatures and densities.16,17,19 However, these models cannot account for the properties of the isolated dimer. Due to the unusually large cooperativity of the HF molecules in the condensed phase, the models fitted to equilibrium properties of the HF dimer1,2,14,19 fail to adequately reproduce properties of liquid HF. The importance of an explicit accounting for polarizability in order to yield better thermodynamical properties of water has been known for almost two decades.21,22 In general, polarizability may be accounted for either by adding flexibility to the molecular geometry21 or by assigning fluctuating atomic charges or dipoles to the interaction sites22 or both. As one might expect, adding such features to three-site models of HF have improved the resulting simulated properties of the liquid.17,19 The polarizable three-site model by Jedlovszky and Vallauri17 (referred to as the JV-P model) seems at present to be the best one in terms of its ability to reproduce the thermodynamics and structure of liquid HF over a wide range of temperature and pressure regions, including supercritical states. Attempts to further account for many-body polarization effects in the condensed phase by incorporating peculiarities of the PES of HF clusters in the vacuum, as obtained from ab initio calculations and experiment, or by introducing smear charges to simulate more explicitly the molecular electronic structure have also been made with varying degrees of success.23-27 However,

10.1021/jp0456685 CCC: $30.25 © 2005 American Chemical Society Published on Web 01/19/2005

6574 J. Phys. Chem. B, Vol. 109, No. 14, 2005 all of these models again assume the presence of additional interaction sites and are therefore computationally relatively costly. One can argue that the majority of these models are “hypothetical” in a sense that they might not represent correctly the interactions at the atomic level in condensed phases of HF. This is even true for models that were fitted to equilibrium properties of clusters, because the interaction in the condensed phase is different from that in small clusters due to polarization effects. The ability of hypothetical models to mimic the properties of HF over a wider region of thermodynamic states and phases is often a matter of how large was the set of HF properties selected to carry out the model parametrization. The task of fitting a model to a possibly larger set of HF properties is obviously more promising if the number of force field parameters is increased (e.g., by introducing additional interaction sites). One indication that a model is hypothetical is that slight changes in model potential parameters or geometry trigger a sharp deterioration in the simulated properties. Such a “poor definition” of the model should not be the case if the model correctly reflects the interatomic forces in the condensed phase. This phenomenon is well-known from water modeling. Many hypothetical water models (e.g., SPC, SPC/E, and TIPnP) show significant (even unphysical) sensitivity toward parametrization.28 Ab initio molecular dynamics (MD) simulations29 show the promise to circumvent the limitations stemming from the uncertainties in the empirical potentials that are inherent in the classical MD simulation technique. For example, there have been a number of ab initio MD simulations of liquid HF.9,30,31 Despite their quite limited statistics, these simulations have provided important insight into the structure and dynamics of liquid HF on a molecular scale that was not achievable by other means. The first Car-Parrinello (CP) MD simulation by Ro¨thlisberger and Parrinello9 was actually for liquid DF, because it is computationally cheaper compared to HF. (This is because the use of deuterium instead of hydrogen permits a CP MD simulation with a larger fictitious electronic mass that, in turn, allows one to integrate the CP MD of the system with a larger time step.) Subsequent ab initio simulations have also been carried out for HF. The most recent of such simulations30,31 demonstrated fairly good agreement with experiments on the structure of HF at different thermodynamic (liquid and supercritical) states. However, as has been recently reported by Grossman et al. for water,32 there is a noticeable dependence of the properties of liquid water on the fictitious electronic mass. They found that water becomes more structured and diffuses more slowly as the fictitious electronic mass decreases. It is interesting that the reported changes of the water properties in the CP MD simulations with a lower fictitious electronic mass look similar to the changes that happen if the ionic temperature in the simulation decreases (e.g., the oxygen-oxygen radial distribution function becomes more structured). One possible explanation of this finding could lie in the fact that in a CP MD simulation the actual inertia of ionic masses is larger as a result of contribution from the inertia carried by the electronic orbitals.33 The fictitious electronic mass in all previous CP MD simulations of liquid HF was above the limit suggested for hydrogen-bonded systems. A promising way to derive an effective empirical force field, which achieves the goal of compromise between an accuracy comparable to ab initio MD simulations and low computational cost, is to fit the force field parameters to a trajectory and force database obtained from ab initio MD simulations. The first such

Izvekov and Voth method was developed by Ercolessi and Adams,34 using leastsquares fitting of the potentials to force data obtained from ab initio calculations. Their method is based on trying to match as closely as possible the first principles forces with those predicted from the classical potentials, and for this reason it has been named the force-matching method. The original force-matching method was successfully applied mainly to elemental systems.35 The application to many-component systems, which include a variety of interactions, appears to be difficult. Recently, we have proposed a new approach for force matching,36 which has proven to be efficient in defining computationally cheap and reliable effective force fields for condensed-phase molecular systems using force and configuration data obtained from ab initio (e.g., Car-Parrinello) MD simulations. The first application36 of this method resulted in a three-site nonpolarizable force field for liquid water at ambient conditions from a CP MD simulation, which showed excellent agreement both with the ab initio MD simulations used to perform the fit and with experiment. In the present work, we have applied our force-matching methodology to derive an effective nonpolarizable two-site pairwise force field for HF using a CP MD simulation of liquid HF at nearly ambient conditions. The CP MD simulation was carried out using an appropriately small value of the fictitious electronic mass. This article is organized as follows. In section II, we outline the details of our approach, fitting procedure, and the CP MD simulation that was used to perform the force matching. Then, section III reports an application of the method to HF and begins by outlining the details of the fitting procedure. The method is applied to derive the effective force field of HF using data from a CP MD simulation of liquid HF at T ) 300 K and F ) 962 kg/m3. The performance of the new model is first tested for the geometry and energetics of the HF dimer, then results on the structural, molecular, transport, vibrational, and thermodynamic properties of liquid HF from the MD simulations employing the new model are presented. In section IV, conclusions from these studies are given. II. Method and Computational Detail The force-matching method used in the present paper is described in detail in ref 36. The method is a variant of the force-matching approach originally proposed by Ercolessi and Adams34 and can be applied if a force field depends linearly on the fitting parameters (a property which can be often achieved by using proper interpolation). In this case, the problem of finding the best least-squares fit to configuration and force data can be accomplished through the solution of an overdetermined system of linear equations. Below, we outline our force-matching method to fit a pairwise central (i.e., dependent only on the atom-atom separation) force field to trajectory and force data from an ab initio MD simulation. The atom-atom force fpi (rij) (acting from particle j on particle i) is broken into a short-ranged part and a longranged Coulomb part as

(

fpi (rij) ) - f(rij) +

)

qiqj r2ij

nij

(1)

where rij is the modulus of the vector rij ) rj - ri connecting two atoms, qi is the partial atomic charge (subject to fit), and nij ) rij/rij. The short-ranged term f (r) is represented by thirdorder polynomials (cubic splines)38 connecting a set of points {rk} (which meshes the interatomic separation up to cutoff radius rkmax), thus preserving continuity of its functions and their first two derivatives across the junction

Force Field for Liquid HF from MD Simulation

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6575

f (r,{rk},{fk},{f′′k}) ) A(r,{rk})fi + B(r,{rk})fi+1 + C(r,{rk})f′′i + D(r,{rk})f′′i+1 (2)

TABLE 1: State Point, Temperature T, Density G, and Pressure p of Hydrogen Fluoride Taken from Refs 18 and 48a

r ∈ [ri,ri+1] where A, B, C, and D are known functions of r and {rk} and {fk} and {f′′k} are tabulations of f (r) and its second derivatives on a radial mesh grid {rk}. A key property for the success of the method is that a spline representation depends linearly on its parameters, which are tabulations of f (r) and its second derivative, {fk, f′′k}, on a radial mesh {rk}. The parameters {fk, f′′k} and qij ) qiqj are to thus be obtained from the fit. Equalization of known ab initio forces Fref il to forces predicted using the representation in eq 1, which act on the same atoms in the lth atomic configuration sampled along ab initio trajectories, results in the following set of linear equations K Nβ

(

∑∑ f(rRil,βjl,{rRβ,k},{fRβ,k},{f′′Rβ,k}) + β)1 j)1

-

qRβ 2 rRil,βjl

)

ref nRil,βjl ) FRil

(3)

R ) 1, ..., K i ) 1, ..., NR with respect to {fRβ,k,f′′Rβ,k,qRβ}, which are force parameters subject to a fit. In eq 3, {Ril} labels the ith atom of kind R in the lth atomic configuration; rRil,βjl is the distance between atoms {Ri} and {βj} in the lth atomic configuration; qRβ is a product of partial charges qR and qβ of atoms of kind R and β; Nβ and K are, respectively, the number of atoms of kind β and the total number of atomic kinds in the system. We will assume that in eq 3 the index l runs over a sufficiently large number L of the atomic configurations such that the equations overdetermine the force parameters. Standard equations that are linear with respect to {fRβ,k,f′′Rβ,k} must be included in eq 3 to ensure that the f (r) first derivative, f ′(r), is continuous across the boundary between two intervals.38 The partial charges qR in eq 3 can be recovered by solving the system of (nonlinear) equations

qRβ ) qRqβ

(4)

where qRβ are parameters obtained from eq 3. Some constraints on the force parameters can also be imposed. For example, it was found to be necessary to use additional constraints aimed at fixing charges on the species to desired values. These constraints can be written as linear equations with respect to force parameters36 and added separately into eq 3. It is also useful to constrain the short-ranged term to zero and optionally to have the first two derivatives to be zero at the cutoff radius (that is the most remote point rRβ,max[k] in a respective spline mesh). However, normally if the grid extends to sufficiently large distances, then the short-ranged term goes to zero smoothly at the cutoff radius resorting to constraints. In the future, we assume that these constraints are added to eq 3. In the original approach by Ercolessi and Adams, the potential parameters are fitted to the whole configuration and force data set using the least-squares algorithm. In our approach, the linear dependence of the force field on the parameters allows one to carry out a fit in a manner that is more consistent with the determination of the mean force field. In our case, eq 3 has to be solved for different smaller sets of atomic configurations, and then the solutions so obtained are averaged over a large number of such sets. Typically, the number of atomic configurations used to build each set of eq 3 is sought to be the smallest possible to ensure that the equations overdetermine the force parameters.

a

state point

T (K)

F (kg/m3)

p (bar)

I II III IV V VI bp

300 373 473 473 473 473 292.7

962 796 796 647 398 326 959

2 12 319 166 84 78 1.013

Values for boiling point (bp) are also given.

It is instructive to give the following interpretation of this method. A reference (total) force acting on each particle in each configuration along the ab initio MD trajectories can be decomposed into a sum of pairwise forces acting on that particle. Obviously, such a decomposition is not unique for a single configuration because the number of atoms in a single set is less than a number of distinct atomic pairs. By addition of more different configurations which are picked up along trajectories, it is possible to determine such a force decomposition, Fij(rij), uniquely. This task is accomplished by solving eq 3. The subsequent averaging, 〈‚‚‚〉, of the force decompositions over trajectories results in a mean (or effective) force field, 〈Fij(rij)〉. III. Results and Discussion A. Details of the ab Initio MD Simulation and ForceMatching Procedure. 1. Details and Results of the CarParrinello MD Simulation. We have carried out a forcematching fit of force field parameters (eq 3) to trajectories and forces obtained from a CP MD simulation of a system of 64 HF molecules. A standard implementation of the CP method was used within the plane-wave-basis density functional theory (DFT),39 using the Becke-Lee-Yang-Parr (BLYP) gradientcorrected exchange-correlation (XC) energy functional40 and the Troullier-Martins (TM) pseudopotential.41 In the past, the BLYP functional has been used extensively for simulations of properties of liquid water32,36,42-45 and liquid HF.9,30,31 The plane-wave cutoff in our simulations was 70 Ry, and the CPMD code, version 3.5, was used.46,47 As was briefly mentioned in the Introduction, the recent simulations by Grossman et al.32 have explored the dependence of the simulated properties of liquid water on the value of the fictitious electronic mass. They found that the fictitious electronic mass used in all previous CP MD simulations of liquid water might cause the radical distribution functions (RDFs) to be artificially less structured. For light water (H2O), the best choice of fictitious electronic mass is less than 350 au. One may expect a similar dependence in a CP MD simulation of liquid HF because the main driving force of structure formation in both systems is hydrogen bonding. It should be noted that due to limited statistics this finding may still be inconclusive and will require further verification. However, for the CP MD simulation of liquid HF in the present work, we used a fictitious electronic mass of 340 au, and the MD was integrated with a time step of 3 au. All previous CP MD simulations of liquid HF used an electronic mass of 800-900 au, which is well above limit suggested by Grossman et al. The system was simulated at phase point I at which the neutron diffraction experiments by Pfleiderer et al.48 have been carried out (Table 1). State point I represents nearly ambient conditions. The corresponding box size for 64 HF molecules in the supercell is 1.3041 nm. The system was equilibrated for 4 ps starting with an initial guess obtained from a classical MD

6576 J. Phys. Chem. B, Vol. 109, No. 14, 2005

Izvekov and Voth

Figure 1. Radial distribution functions of liquid HF from the CP MD simulation of a system of 64 HF molecules at state point I (Table 1).

simulation using the nonpolarizable three-site model of Jedlovszky and Vallauri.16 During the equilibration run, the desired temperature was maintained using a simple velocity rescaling method. The thermostat was then removed, and the system was allowed to evolve in the microcanonical ensemble (i.e., without any temperature control imposed). The resulting production run was for 12 ps. The average ionic temperature in the production portion of the simulation was T ) 299 K. The partial radial distribution functions (RDFs) gFF, gFH, and gHH, together with the total radial distribution function gtot, are shown in Figure 1. The gtot is measured in neutron diffraction experiments and can be expressed through partial distribution functions48

gtot ) 0.4966gFH + 0.2104gFF + 0.2930gHH

(5)

Only the intermolecular parts of the RDFs are shown. In general, gFF and gFH show reasonable agreement with those recently reported by Raugei and Klein31 from a CP MD simulation of 54 HF molecules at T ) 290 K (the gHH is not shown in their paper) in which a fictitious electronic mass of 800 au was used. In particular, there is a well-pronounced first minimum in the gFF of a characteristic shape. It should be noted that such a minimum is absent in the earliest CP MD simulation by Ro¨thlisberger and Parrinello.9 However, the total RDF from our CP MD simulation is overstructured compared to the experimental one.48 As will be discussed later, when the results of classical MD simulations using the effective force field from the force-matching method are presented, the system in our CP MD is still rather far from equilibrium due to the limited duration of the simulation (which was just around 10 ps). Due to the strong hydrogen bonding, the liquid HF is hard to equilibrate in the simulation, which becomes even more difficult for a system with a small number of molecules as is a case for the present CP MD simulation. The MD simulation using the forcematching model, which can be run for much longer times, suggests that a well-equilibrated system has less structure. In particular, the second minimum in gFF is less pronounced, and

the first peak is less intense. The self-diffusion coefficient derived from the CP MD simulation is 3.0 × 10-9 m2/s, which is significantly smaller than the experimental value of about 10.0 × 10-9 m2/s (and the MD simulation using the forcematching model), a feature that can also be explained by the limited equilibration. The fact that the system is not completely in equilibrium should have little effect on the quality of the force-matching procedure, however. 2. Details of the Force-Matching Procedure. The short-ranged force component of the effective HF force field (eq 1) was represented by a spline over a mesh with a grid spacing of approximately 0.0025 nm in the region of the intramolecular degrees of freedom and 0.005 nm in the region of the intermolecular motions. The interatomic separation to which the meshes extend represents the cutoff radius for the short-ranged interaction. To achieve smoothness, boundary conditions on the short-ranged force and its first two derivatives were applied to make them zero at the cutoff radius. We have used the same cutoff of 0.78 nm for all FF, HF, and HH pairs. An attempt to use different values of the cutoff for different interactions (e.g., a smaller cutoff for the HF and HH interactions) resulted in a pronounced deterioration of the simulated properties. This effect was linked to the fact that the FF short-ranged interaction was not counterbalanced by the HF interaction beyond the HF cutoff radius. Similar behavior was observed if the cutoff for the HF interaction exceeded the FF cutoff. Similar to the water case,36 this behavior seems to indicate that the short-ranged interaction at larger atom-atom separations as obtained from the fit is due to the polarization effects and an associated redistribution of the effective charges over the fluorine and hydrogen. The cutoff for the short-ranged interaction must therefore be of an adequate value to accommodate the polarization effects properly. A cutoff of 0.78 nm chosen here does satisfy such a requirement. To fit the force field parameters to data from the CP MD simulation, the overdetermined system of linear equations (eq 3) was solved repeatedly for 1360 different sets of atomic configurations, with each set consisting of two configurations (l ) 1, 2 in eq 3). The latter was equidistant in time and

Force Field for Liquid HF from MD Simulation

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6577

Figure 2. Effective pairwise forces in liquid HF as functions of interatomic separation calculated by the force-matching method. The force parameters were fitted to the trajectory and force data obtained from a Car-Parrinello MD simulation of a system of 64 HF molecules as described in the text and with a cutoff for the short-ranged interaction of 0.78 nm (FF(64,0.78) model). Indices “b” and “nb” label, respectively, the bonded and nonbonded components of the force field.

separated by approximately 4.36 fs. This corresponds to a total length of trajectories of about 12 ps of the CP MD simulation, and then the resulting solutions were averaged over all sets. To match the electrostatics performed in the CP MD simulation, the Coulomb forces in the force-matching procedure were calculated using the metallic boundary conditions to the Ewald sum.49 The Ewald sum is necessary to correctly calculate the coefficients before parameters qRβ in eq 3. The force field obtained from solving eq 3 often substantially fluctuates among different sets of atomic configurations. However, although different radial meshes for short-ranged parts of the forces or different samplings of the atomic configurations were used, upon averaging over an adequately large number of atomic configuration sets, the effective force field was seen to be quite reproducible. B. HF Force Field. In Figure 2, the force field fitted to data from the CP MD simulation of a system of 64 HF molecules with an applied cutoff of 0.78 nm is shown. This force field is referred to as the FF(64,0.78) model. The intermolecular shortranged part of the FF(64,0.78) force field, f(r), is displayed in Figure 3a. The spline parameters ri, fi, and f′′i from the fit (eq 2) are summarized in Table 1S in the Supporting Information and can be used to implement the model. Alternative implementations using the fit of the spline data by polynomial series are presented later. Several conclusions can be drawn upon inspection of the plot. First, the f(r) is relatively complex compared to water.36 As discussed later, in contrast to water, an analytical fit of a spline representation of f(r) by a truncated polynomial series appeared to be difficult. A second observation is that the f(r) is rapidly vanishing around 0.4 nm. Above this separation the interaction is dominated by a Coulomb term. However, the short-ranged interaction at larger interatomic separations may contribute substantially to the thermodynamic properties of liquid HF. The MD simulation in the NVT ensemble with a cutoff of 0.4 nm did not show any meaningful differences in the simulated RDFs and self-diffusion properties compared to a MD simulation using a cutoff of 0.78 nm. However, the thermodynamic properties were poor in the simulation with a smaller cutoff apparently as a result of missing short-ranged attractive interactions. It was therefore decided to keep a cutoff of 0.78 nm to ensure correct thermodynamic, structural, and dynamical properties. The short-ranged interaction at larger separations, when it is already almost vanishing, might originate in part from polarization effects. Another possibility is an artificial mixing between the short-ranged and Coulomb terms, which is an

Figure 3. (a) Short-ranged nonbonded component of the force field shown in Figure 2. The spline parameters are given in Table 1S in the Supporting Information. (b) Least-squares fit to the short-ranged forces shown in part a using a truncated Chebyshev series (eq 8). Insert shows the enlarged region of FF force in the area of the minimum and its fit by the Chebyshev expansion (thick solid line in part b). The coefficients of the expansion in the Chebyshev polynomials are given in Table 3.

TABLE 2: Partial Charges, Intramolecular Potential Parameters, and Properties of the Gas-Phase Monomer from the FF(64,0.78) Modele FF(64,0.78) qF, qH (au)

ab initio (BLYP)

expt

-0.4073062, 0.4073062 0.09472 0.0933a 0.09168b

r0FH (nm) kFH (kJ/(mol Å2)) 4045.19 µ0 (D) 1.858 ω (cm-1) 3341.01

1.80a

6206.29b (3988.21c) 1.826d 4138.32b (3317.39c)

a Reference 9. b Reference 52. c Reference 8 (liquid HF, 293 K). Reference 53. e Shown are qF and qH, the partial charges on fluoride and hydrogen, r0FH, the HF equilibrium distance, KFH, the intramolecular force constant, µ0, the molecular dipole moment, and ω, the harmonic frequency. For the intramolecular force constant and frequency entries from experiment, two values are shown. The first is for the HF monomer in the gas phase, and the second (in brackets) is as estimated from Raman and IR spectra of liquid HF.

d

inherent in the force-matching procedure and manifests itself more strongly for smaller systems. The small magnitude of the short-ranged forces beyond 0.4 nm may also indicate less significance of polarization effects in liquid HF compared to water.36 Other evidence of the relative unimportance of polarization effects in liquid HF is as seen in Table 2, where the condensed-phase dipole moment estimated from the forcematching partial charges is very close to the dipole moment of the gas-phase monomer. For water, the difference is more pronounced.36 The short-ranged component of the water forcematching force field is also much more evident in the background of the Coulomb forces at larger separations.36 Inserts to Figure 3 show an enlarged region of the FF shortranged force around the distance at which the force is zero (i.e., the corresponding potential has a minimum). The force profile has a characteristic shape with well-pronounced oscillations. These oscillations apparently are not due to, for example, the limited statistics as they have a tendency to become well-defined when even more configurations are used in the fit. However, it

6578 J. Phys. Chem. B, Vol. 109, No. 14, 2005

Izvekov and Voth

is difficult to say if the same is true for the ripples in the plateau around 0.2 nm of the HF short-ranged force profile (Figure 3a). This behavior becomes less prominent as the force averaging proceeds. However, whether it vanishes completely after an averaging over much longer trajectories remains unclear. At small separations, r < rcore, which are not explored in the CP MD simulation, the short-ranged forces should be extrapolated. Fortunately, the functional form to represent the forces at small separations appears to be an unimportant factor in the simulations of liquid HF over a wide range of thermodynamic conditions. In our simulations, the forces were assumed to be constant for r < rcore and continuous, i.e.,

fij(r < rijcore) ) fij(rijcore)

(6)

The intramolecular HF force is approximated by a linear function (i.e., the HF potential is assumed to be harmonic) as 0 FFH intra ) -kFH(r - rFH)

(7)

where the r0FH is the zero of the FHb spline fit shown in Figure 2 (i.e., the equilibrium HF distance in a vacuum). However, because the FHb force deviates in reality from being harmonic (which can be inferred already from Figure 2), its parametrization can be accomplished in two ways. The first one is a fit in a least-squares sense, while the second one is by a function having the same slope at zero as does the spline approximation. The latter parametrization produces a vibrational frequency equal to that derived from the spline fit. This frequency appears to be larger compared to the first method of parametrization and very close to the experimental position of approximately 3400 cm-1 of the band associated with the intramolecular and in-chain vibrations in the IR spectrum of liquid HF.8 The second method for parametrization was therefore chosen, and the intramolecular force parameters are summarized in Table 2. As seen from the data, the derived intramolecular force constant and corresponding frequency of stretching vibrations are in virtually exact agreement to experimental data at T ) 293 K for liquid HF. Partial charges derived from the fit are shown in Table 2. The resulting dipole moment of the gas-phase monomer is very close to the experimental value, suggesting that the force field that was parametrized to the liquid phase at state point I may have reasonable transferability to other states/phases. The fact that the effects of the environment on the HF molecules do not seem to be large has been noted by other authors.19 Similar to the water case,36 one might also expect that the intermolecular part of the FF(64,0.78) force field will work well if the intramolecular interaction is described by a more elaborate potential, determined independently from the force-matching. The HF bond length, r0FH, is close to the experimental value and ab initio BLYP calculations for the monomer in a vacuum. The discrepancy is in the range of bin size of the mesh used in the region of the intramolecular movements (approximately 0.0025 nm). More compact analytical representations for short-ranged nonbonded forces, f(r), which are easier to implement can be obtained by fitting the spline representation of the f(r) (shown in Table 1S in the Supporting Information) by a truncated series of polynomials. Such a polynomial fit smooths out roughness in the original force profiles, which is a remnant of the limited CP MD simulation data used to carry out the fit. The most rapid convergence is provided by the Chebyshev series.50 The shortranged forces were represented as

16

fij(r) ) A0ij/2 +

∑AnijTn(rjij) n)1

(8)

where Tn(rj) is the Chebyshev polynomial of degree n and

jrij )

(2/r - 1/rijcore - 1/rcut) (1/rijcore - 1/rcut)

The Chebyshev expansion is performed with respect to the variable 1/r that dictates the choice of jrij. The coefficients were then determined from a least-squares fit of the tabulation of f(r), which was obtained using a spline representation on a fine mesh. The resulting coefficients are given in Table 3 with the values of rijcore and rcut summarized in the footnote. The order of the expansion is the smallest one to yield a satisfactory reproduction of the properties observed in the MD simulation using the original spline representation. The respective FF and FH profiles are plotted in Figure 3b. The peculiarities of the FF(64,0.78) model make it difficult to have an approximation by polynomials as good as for water.36 For example, some structural properties, such as the characteristic split of the second hump in gFF, cannot be fully reproduced in the MD simulation using the representations by polynomials, which will be discussed later. A further increase in the order of the Chebyshev expansion remedies the problem slowly. The FF(64,0.78) model can be implemented more straightforwardly using an expansion in powers of 1/r. Table 4 summarizes the coefficients of the least-squares fit of the shortranged force terms using powers of 1/r up to 16 as basis functions. A higher order expansion does not show any meaningful improvement in the simulated properties. This representation shows less quality in the simulated properties of liquid HF compared to the representation by the Chebyshev expansion (eq 8). The performance of the polynomial representations of the FF(64,0.78) model will be discussed in more detail later. TABLE 3: Coefficients Anij of the Least-Squares Fit of Short-Ranged Nonbonded Forces fij(r) of the FF(64,0.78) Force Field Using the Expansion in the Truncated Chebyshev Series (Eq 8)a n

AnFF

AnFH

AnHH

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

0.193342488978e-1 0.170695211522e-1 0.117913019161e-1 0.571709492757e-2 0.112514819828e-2 -0.772549113885e-3 -0.655577763399e-3 0.123526028262e-3 0.520439383903e-3 0.388859135800e-3 -0.333842573042e-4 -0.498698104328e-3 -0.716838763770e-3 -0.567003900883e-3 -0.269039513926e-3 -0.484950667050e-4 0.433827972583e-4

0.621985153829e-3 0.109026063736e-3 -0.606850882024e-4 0.103623019267e-2 0.105964488573e-2 -0.134667022858e-3 -0.329883071012e-3 0.146849787350e-3 0.296795728745e-3 0.191634768985e-3 0.215389228541e-4 -0.188075944017e-3 -0.183335626847e-4 0.241548229705e-3 0.128823377093e-3 -0.298262317853e-4 0.626369657173e-4

-0.268472919071e-2 -0.867358654931e-3 0.172482739597e-2 0.121445463118e-2 -0.746712122832e-3 -0.969526094470e-3 0.102579469178e-3 0.602754074649e-3 0.283195537118e-3 0.171839984403e-4 0.132095715940e-4 -0.105013965380e-3 -0.292424029834e-3 -0.214605799048e-3 -0.718183451234e-5 0.612342076714e-5 -0.159159766495e-4

a Fitted data are spline representations of the f (r), which are shown ij in Table 1S in the Supporting Information and Figure 3a. The following FH HH radii, rijcore and rcut, were used: rFF core ) 4.2 au, rcore ) 2.25 au, rcore ) 2.8 au, and rcut ) 14.74365 au. Atomic units were used. The leastsquares fit is plotted in Figure 3b. The forces were extrapolated to small separations r < rcore as in eq 6. The partial charges are shown in Table 2. A cutoff of 0.78 nm must be applied to this expansion.

Force Field for Liquid HF from MD Simulation

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6579

TABLE 4: Coefficients Anij of the Least-Squares Fit of Short-Ranged Nonbonded Forces fij(r) of the FF(64,0.78) Force Field 16 Anij/rn a Using the Expansion fij(r) ) ∑n)2 n

AnFF

AnFH

AnHH

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

4063.54855662332 -336360.890958584 12343366.2070739 -264698607.262515 3677721011.97818 -34636485585.0036 224154592385.060 -983927047839.102 2786305842800.11 -4436942782374.86 2131577820637.39 3030645295947.51 -866029041759.284 -3324110599402.83 -3308086287978.46

-261.892583603018 20047.8950986942 -684444.695203338 13819820.8866951 -184549722.206600 1726585010.85318 -11689875203.4294 58294017890.4409 -215572710357.268 589466124012.250 -1175482409379.57 1660743758449.04 -1574161540435.74 897501162822.199 -232568485731.251

-206.174243434062 15313.5408951016 -502027.369816623 9591242.63437633 -118687327.548650 998919440.279957 -5829409306.93262 23497961691.8941 -63327395771.1718 104060624192.090 -74662502439.6124 -38932796993.1336 75060967782.0464 44550205038.8726 -83752945658.0223

a Fitted data were spline representations of the fij(r) which are shown in Table 1S and Figure 3a. Atomic units were used. At small separations FH HH r < rcore the fij(r) was extrapolated in accordance with eq 6. The following core radii were used: rFF core ) 4.2 au, rcore ) 2.25 au, and rcore ) 2.8 au. The partial charges are shown in Table 2. The cutoff of 0.78 nm must be applied to this expansion.

TABLE 5: Properties of the Gas-Phase HF Dimer, H1-F1‚‚‚H2-F2 from the FF(64, 0.78) Modele FF(64,0.78) |F1‚‚‚F2| (nm) |F1‚‚‚H2| (nm) |H1-F1| (nm) |H2-F2| (nm) R[∠(F1‚‚‚H2-F2)] (deg) R1[∠(H2-F2‚‚‚F1)] (deg) β[∠(H1-F1‚‚‚H2)] (deg) β1[∠(H1-F1‚‚‚F2)] (deg) µdimer (D) µmean (D) Upot (kJ/mol)

0.2560 0.1597 0.0948 0.0967 173.2 4.2 111.1 108.5 2.97 1.88 -20.05

JV-Pa

ab initio (BLYP)

expt

171.4 5.5 110.2 107.1

0.2797 0.272c 0.1846 0.0952 0.0958 171.9 5.3 10 ( 6c 111.9 109.2 117 ( 6c

-20.35

-17.99b

0.2730 0.1764

-19.08d

a References 17. b Reference 9. c Reference 54. d Reference 55. e The bond lengths and angles, total dipole moment µdimer, average molecular dipole moment µmean, and the interaction energy Upot are shown. The dimer properties from the model by Jedlovszky and Vallauri (JV-P)17 and ab initio calculations shown for a comparison.

Figure 4. (a) Effective nonbonded atom-atom potential derived from the FF(64,0.78) parametrization (solid lines). (b) Short-ranged part of the potential shown in part a.

In the simulations reported here, the short-ranged potentials that are necessary to determine the simulated energetics and thermodynamics were calculated by integrating the respective terms in the force fields and then shifting them to zero at the cutoff radius. These manipulations might be viewed as a source of error in the simulated HF energetics. However, the way in which the potentials are calculated should not influence the behavior of HF in the MD simulations because in our framework the short-ranged forces are not derived from the potentials. Moreover, as shown below, the calculated thermodynamic properties and energetics using potentials obtained in this way show reasonable agreement with experiment. Potentials calculated using the short-ranged parts of the FF(64,0.78) forces integrated out in the manner just described are shown in Figure 4. The classical MD simulations with the new force-matching model were performed using the DL_POLY2.9 computer code.51 The tabulated force fields and potentials for the various

representations of the model are available in a DL_POLY code format upon request. C. HF Dimer. The simplest associated form of HF is the gas-phase dimer, (HF)2, the structure of which has been studied extensively by experimental techniques54,55 and ab initio methods.56,57 In the equilibrium configuration of the isolated dimer H1-F1‚‚‚H2-F2, where “‚‚‚” stands for a hydrogen bond between two HF molecules, the two molecules are coplanar but bent. The atoms F1‚‚‚H2-F2 are arranged nearly linearly with the hydrogen placed slightly off the FF axis. The distance between the fluorine atoms is rFF ≈ 0.27 nm. Such a bent structure is believed due to a competition between the dipolar and quadrupolar interactions.18,57 The quadrupole moment of the monomer has been regarded as an essential ingredient of the model to reproduce the bent structure of the HF dimer and consequently to ensure the formation of zigzag chains in the liquid and solid phases. The absence of a molecular quadrupole moment in previous two-site models was suggested to cause their failure to predict the correct structure of liquid HF as well as the reason for the better performance of three-site models, which possess a molecular quadrupole moment. Surprisingly, as seen from Table 5, the FF(64,0.78) force field with just two sites predicts a bent structure of the HF dimer, which is in a good agreement with both experiment and ab initio calculations. The dimer properties from the JV-P model (the

6580 J. Phys. Chem. B, Vol. 109, No. 14, 2005 model assumes a rigid molecular structure), which is presently one of the best for simulations of the condensed phases of HF, are shown for comparison. The dimer ab initio calculations were carried out using the same cutoff, XC functional, and pseudopotentials as in the CP MD simulation. More accurate studies of the PES of the gas-phase dimer are available.23 As will be shown later, the FF separation in (HF)2 from the FF(64,0.78) model, which is shorter than the experimental result, is identical to a position of the first peak in the liquid-state gFF. This suggests that the discrepancy of the properties of (HF)2 between those from experiment and ab initio calculations can be largely explained by the fact that the FF(64,0.78) model represents the effective interactions in the liquid phase. It is reasonable to expect then that the performance of the FF(64,0.78) force field will deteriorate once the HF molecules are absent from their liquid environment. The concept that the FF(64,0.78) parametrization bears a “memory” of the liquid phase can be also invoked to explain the shorter F1‚‚‚H2 distance as well as the slightly more pronounced difference in intramolecular bond lengths compared to those of both experiment and ab initio calculations. It is interesting that the observed larger difference of intramolecular bond lengths is not because of a softer intramolecular force constant k (eq 7 and Table 2), which results in a lower frequency of stretching vibrations of the HF monomer in a vacuum. A harder intramolecular potential, such that the frequency of the stretching mode in a vacuum coincides with the experimental results, resulted in a similar difference of intramolecular bond lengths within (HF)2. The dimer properties in Table 5 are identical to those if the FF(64,0.78) short-ranged forces are described using the analytical fits by polynomials (Tables 3 and 4). D. Liquid HF. 1. Details of the MD Simulations. MD simulations were carried out for 256 HF molecules in a cubic supercell. The use of the least-squares fit of the FF(64,0.78) short-ranged forces by polynomials (Tables 3 and 4) instead of the spline fit resulted in slightly different properties of liquid HF. The difference is especially noticeable in the gFF. The representation using an expansion in Chebyshev polynomials works generally better than the expansion as powers of 1/r. The results presented here are from simulations using the spline representation. It will be explicitly stated when the properties of the liquid HF using the analytical fits are different. The Coulomb forces were calculated using the Ewald method, the temperature was weakly coupled to Nose´-Hoover thermostat with a relaxation time of 0.1 ps, and the pressure in the simulations in the NPT ensemble was weakly coupled to a barostat with a relaxation time of 0.5 ps. Because the force-matching model is computationally efficient, it is possible to study properties of the liquid HF on a multi-nanosecond time scale. The MD system was typically equilibrated for 1 ns, and a simple velocity scaling method was used to control the system temperature during the equilibration.49 Production runs were then 2-3 ns long. The MD time step was 0.5 fs. 2. Structural Properties. In Figure 5, the partial radial distribution functions from a MD simulation using the FF(64,0.78) model are compared to those from the CP MD simulation used to perform the force-matching. It is seen that there is apparent disagreement between the distribution functions from these two simulations. The reason for disagreement, which was briefly discussed earlier, is an insufficient equilibration of the system prior to the CP MD production run. Due to the strong aggregation in the HF liquid phase, the system is hard to equilibrate, and a much longer equilibration run in the CP MD

Izvekov and Voth

Figure 5. Comparison of the radial distribution functions gFF, gFH, and gHH for the liquid HF at state point I from a CP MD simulation (thin line) and from a MD simulation in the NVT ensemble using the FF(0.78,64) model (thick line).

simulation may be needed. This conjecture is supported by the observed dynamics of structure formation in the force-matching MD simulation of a system of 64 HF molecules initiated from the atomic configuration, which is the final one in the CP MD production run. To perform the MD simulation on such a small supercell, the force-matching force field was reparametrized using a smaller cutoff of 0.63 nm to the short-ranged term (since the cutoff should be less than half of the supercell). The comparison of the RDFs from the CP MD simulation and from first 10 ps of the force-matching MD simulation of a system of 64 HF molecules is shown in Figure 6. There is a strong resemblance in the RDFs, indicating that indeed the system has still not reached equilibrium at the end of the CP MD simulation. In the force-matching MD simulation, the RDFs adopt their equilibrium shapes, which are close to those shown in Figure 5, after approximately 40-50 ps of the simulation. However, we would like to stress again that an incomplete equilibration of the CP MD system does not harm the quality of the forcematching procedure. This is because the interatomic forces are essentially the same in the different atomic configurations sampled from the CP MD simulation. It is interesting to note that qualitatively the RDFs from the MD simulation using the FF(64,0.78) model look similar to the RDFs predicted using the three-site model developed by Klein, McDonald, and O’Shea14 (HF3 model) by fitting to the PES of (HF)2. For example, there is a characteristic split in the second hump of the gFF from both models. Such a split is absent in the gFF from almost all other empirical models developed so far. The second hump in the gFF originates from correlations between nearest-neighbor F atoms, which belong primarily to different chains. However, the agreement in the structure of the liquid HF between the two models is only qualitative. For example, the amplitudes of the RDF maxima are different. The second hump in the gFF of the FF(64,0.78) model is more symmetrical and does not end with a steep extended slope as predicted by

Force Field for Liquid HF from MD Simulation

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6581

Figure 6. Comparison of partial radial distribution functions gFF, gFH, and gHH and the total radial distribution function gtot (eq 5) for the liquid HF from the same CP MD simulation as in Figure 5 (thin line) and after 10 ps of the force-matching MD simulation of the same system as in the CP MD simulation (thick line). The force-matching MD simulation was initiated from the configuration that is the final one in the CP MD production run.

the HF3 model. Many models (e.g., the HFC model by Cournoyer and Jorgensen,15 the models by the RHNC approach,20 and the QCT potential27) do not show meaningful structure in gFF beyond the first peak, especially at close to ambient conditions. Here, we further discuss the RDFs from the FF(64,0.78) model at state point I, which are shown in Figure 5. The pair distribution functions that give the information about the structure of hydrogen bonding are gFF and gFH. The gFF has a first peak at 0.256 nm. It is interesting to note that this distance coincides with the FF separation in the gas-phase dimer predicted by the FF(64,0.78) model (Table 5). The same property was observed for the water force-matching force field (i.e., the OO distance in the gas-phase dimer is very close to the position of a first peak in the gOO).36 The coordination number up to first minimum in the gFF is 2. The gFH RDF has a first intermolecular peak at 0.16 nm. Again, this distance is in almost exact agreement to the length of the hydrogen bond in the gas-phase dimer from the FF(64,0.78) model (Table 5). The coordination number at the first minimum beyond the intramolecular peak is 1, indicating that each HF molecule is involved in a hydrogen bond. Importantly, the first peaks in the gFF and gFH are virtually of the same height, pointing to highly correlated F-F and F-H movements of the hydrogen-bonded neighbors. The gFH has a pronounced second hump, which is followed by a small dip. This maximum corresponds to F-H nearest-neighbor correlations between different chains and is observed experimentally in the gtot. As was mentioned earlier, the analytical fit of the spline representation of the FF(64,0.78) short-ranged component by polynomials is not able to completely capture all of its features. The expansion in the truncated Chebyshev series (eq 8) produces a better fit, however, compared to the expansion into powers of 1/r. Consequently, the analytical representations of the FF(64,0.78) model perform somewhat differently in the MD

Figure 7. Comparison of gFF from MD simulations using different representations of the short-ranged forces of the FF(64,0.78) model. Shown are the original spline fit (thick solid line), the fit of the spline data using a truncated Chebyshev series (thin solid line), and the fit of the spline data using expansion into powers of 1/r (dashed line). Insert shows the enlarged region of the shoulder after a first peak. The coefficients of analytical fits are given in Tables 3 and 4.

simulation. In Figure 7, gFF from an approximation of the FF(64,0.78) model using the polynomial expansions is compared to gFF functions using the original spline representation. The expansion in the Chebyshev polynomials produces gFF that agrees well with gFF simulated using the spline fit, with the only exception being that the second hump is unsplit. It is difficult to say which part of the original spline representation is responsible for such a split, which obviously is not properly approximated by the Chebyshev expansion. The expansion as powers of 1/r produces an even more different structure. However, the positions of the major maxima and minima in the RDFs by the different representations are almost the same. The differences in the liquid structure of HF from the different fits are much less pronounced in gtot. There is considerable disagreement among the different models and approximations for the partial radial distribution functions of liquid HF. The

6582 J. Phys. Chem. B, Vol. 109, No. 14, 2005

Izvekov and Voth

Figure 8. Comparison of the total pair distribution function gtot of HF as obtained from neutron diffraction experiments48 (thick line) and from MD simulations under NVT conditions using the FF(64,0.78) model (thin line) at six different thermodynamic state points (Table 1).

total RDF was separated experimentally into partial RDFs only recently,58,59 so many authors have published just gtot from their simulations. The total radial distribution functions of DF (eq 5) at different state points (Table 1) obtained from the most recent neutron diffraction experiments of Pfleiderer et al.48 are compared in Figure 8 with the corresponding functions of the FF(64,0.78) model. The simulated gtot clearly shows four peaks that come from the first peaks of gFH, gHH, and gFF and a second peak of gFH as the atom-atom separation increases (Figure 5). The agreement with experiment is seen to vary with the state conditions. Reasonable agreement is achieved for the I-IV states. The simulated gtot routinely shows a small second peak that nearly coincides with the position of the first dip in the experimental gtot. However, the simulated gtot at T ) 300 K is in slightly better agreement with the experimental gtot at T ) 293 K reported by Deraman et al.4 and more recently by McLain et al.,58,59 which shows such a small second peak.20 The partial RDFs exhibit reasonable agreement with those first reported by McLain et al. using neutrons. However, the first intramolecular peaks in the simulated gFF and gHF are more intense (by about 15%). The height of the first peak in gHH agrees well with the experimental results. The positions of the peaks are also in good agreement with experiment. To conclude this section, we will present an analysis of the effect of rigidity in the HF molecular structure on the simulated liquid structure of HF. Almost all empirical models developed so far assume a rigid HF molecular structure. In Figure 9, gFF of a rigid version of the FF(64,0.78) model is compared to gFF of a flexible one. The most noticeable consequences of the rigidity of molecules are seen to be a less intense first peak and somewhat different second hump. The gFH and gHH also have their first peaks proportionally less intense as in gFF. 3. Molecular Properties and Aggregation. The characteristic feature of liquid HF responsible for many of its properties is the apparent existence of hydrogen-bonded chains of monomers. The bent structure of the dimer in a vacuum should also be

Figure 9. Comparison of gFF of rigid (thick line) and flexible (thin line) versions of the FF(64,0.78) model from the MD simulation in NVT ensemble of the liquid HF at state point I (Table 1).

somehow manifested in the molecular complexes formed in the condensed phase. Indeed, a best-fit analysis of dielectric constant data and Raman and IR spectra suggests that aggregation of monomers into acyclic zigzag chains is common in the liquid. The same structure is found in solid DF.5 The ability of the present force-matching model to reproduce the bent structure of the gas-phase dimer forms the basis for an adequate treatment of complexation in the liquid phase. To identify chains, we have used a simple algorithm identical to that used in several previous works.9,30 The algorithm scans all monomers that do not belong to any chain previously identified. For each such monomer, it walks to both ends of the prospective chain. The next monomer in the chain is taken to be that which forms the shortest F‚‚‚H bond with the latest monomer added to the chain and with a hydrogen bond length within a preselected cutoff radius rcut. This algorithm allows one also to identify rings as well as chains having loops (i.e., a ring connected to the end of chain). However, the latter was never observed in our MD simulations using the force-matching model. The results yielded by the search algorithm are also dependent on the choice of rcut. A value of rcut of 0.225 nm was used, which is the position of the first minimum in gFH. Figure 10 shows histograms of distributions of chain lengths

Force Field for Liquid HF from MD Simulation

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6583

Figure 11. Probability densities of distributions for the hydrogenbonding angles between monomers in chains at T ) 203 K (thin line) and T ) 300 K (thick line) (Table 5).

Figure 10. Histograms of (a) chain and (b) ring length distributions, P(n), at T ) 203 K (thin line) and T ) 300 K (thick line).

TABLE 6: Properties of Chains and Rings from MD Simulations in the NVT Ensemble at Different Temperatures and Experimental Densitiesa T ) 203 K property

T ) 300 K

chain

ring

chain

ring

F-H dintra (nm) F‚‚‚H dbond (nm)

0.09691

0.09688

0.09654

0.09664

0.16278

0.16499

0.16997

0.17477

R (deg) β (deg) FF dinter (nm)

162 119 0.3518

156 118 0.3554

FH dinter (nm)

0.3058

0.3091

HH dinter

0.3425

norm dend-end (nm)

7.45 137 34.27 0.4871 0.07696

(nm) L Lmax N dend-end (nm)

0.3441 6.44 17 1.15

4.87 73 52.40 0.4351 0.08190

3.94 20 1.48

F-H The notation is as follows: dintra , average intramolecular FH bond F‚‚‚H length within chain/ring; dbond , average hydrogen bond length within chain/ring; R and β, average intrachain R and β angles (Figure 11 and FF FH HH Table 5); dinter , dinter , and dinter , average distances between nearest atoms of corresponding types in the adjacent chains; L, average number of monomers per chain/ring; Lmax, maximal number of monomers in chain/ring observed in the MD simulation; N, average number of chains/ rings in the MD simulation; dend-end, average chain end-to-end distance; norm , average chain end-to-end distance per monomer. dend-end a

and ring sizes at two temperatures. The first temperature is close to the freezing point and the second one corresponds to nearly the boiling point (Table 1). Table 6 summarizes some more properties of the chains and rings. The properties of the rings were obtained by averaging over only those configurations along the MD trajectories in which the rings were observed. A ring formation is a relatively rare event. In the MD simulation at T ) 203 K, rings were observed just in 8% of the total simulation time, while a higher temperature favors ring formation. For example, in the MD simulation at T ) 300 K, the rings were observed in 15% of the total number of configurations sampled along the trajectories. The average chain length (L entry in Table 6) is 7.45 at T ) 203 K, but it decreased

to 4.87 at T ) 300 K. A similar tendency was observed for rings. A best-fit analysis of Raman and IR spectra8 suggests that the average length of chains at T ) 293 K is L ) 6-7, which increases to L ) 8 as the freezing point (T ≈ 190 K) is approached and then becomes infinite in the solid phase. The average intramolecular and hydrogen-bond distances (dF-H intra and dF‚‚‚H bond ) within the rings are slightly larger compared to those of the chains. The average number of chains and rings (N) grows with temperature. The increase in chain number N at higher temperature is likely due to a fracturing of long chains. The larger number N of rings accompanied by their smaller average size (L) at a higher temperature suggests an enhancement in the conformational dynamics of the chains. The end-to-end distances, dend-end, at T ) 203 K and T ) 300 K, are seen to be close despite the fact that the average chain length at T ) 203 K is more than 50% larger. This indicates that longer chains at lower temperatures tend to be more curled. An increase in the norm , at T ) 300 K end-to-end distance per chain monomer, dend-end is likely due to changes in intrachain angles between the monomers. As seen from Table 6, an average distance dFH inter of about 0.3 nm between nearest H and F atoms in the adjacent chains almost coincides with the position of a second peak in gFH (Figure 5). This peak comes from non-hydrogen-bonded FH nearest neighbors. The dFF inter is close to the position of the end of the shoulder after the first peak in gFF, suggesting a broader dispersion in the FF nearest-neighbor interchain distances than just the width of the shoulder. This distance is around 0.35 nm and slightly raises with temperature. The experimental result shows that in solid DF the interchain distance is longer at 0.312 nm.5,48 The orientation correlations of monomers within the chains is an important property to validate the model. Figure 11 shows the probability density of the cosine distributions ∠F‚‚‚H-F (R) and ∠H-F‚‚‚H (β) at different temperatures (the abscissa is shown in degrees for convenience). The average values of the angles are summarized in Table 6. The distributions qualitatively agree with those determined by some other advanced models16,27 (e.g., JV-P). The probability density of R has a maximum at 180° with an average value of about 160°. The average β is closer to the position of the maximum in the probability distribution (about 120°). This picture correlates well with the experimental data on the geometry of chains in the solid. Chains in solid DF are parallel with an R angle of 180° and β of about 120°.5 The temperature increase naturally results in a broader angle distributions. However, the difference in average R and β angles between T ) 203 K and T ) 300 K is relatively small (6-10°).

6584 J. Phys. Chem. B, Vol. 109, No. 14, 2005

Izvekov and Voth

TABLE 7: Self-Diffusion Coefficients (10-9 m2/s) of Liquid HF at Six Different Temperatures from MD Simulations in the NVT Ensemble at the Experimental Densities and from Experimenta T (K)

FF(64,0.78)

expt

203 213 233 253 273 300

1.95 2.50 3.90 5.95 7.10 10.1

1.33, 1.27 1.80, 1.67 3.00, 2.70 4.65, 4.05 6.74, 5.72 10.3, 8.47

a For the experiment, two values are shown; the best fits are from refs 63 and 64, respectively.

Figure 12. Power spectra of liquid HF from the CP MD simulation (thick solid line) and the same force-matching MD simulation as in Figure 5 (thin line) and Figure 6 (dashed line).

4. Dynamical Properties. In Table 7, the self-diffusion coefficients of liquid HF using the FF(64,0.78) model (spline representation) are compared with the experimental results. The MD simulations were carried out in the NVT ensemble at experimental densities60,61 (the thermostat relaxation time was increased to 0.5 ps). The polynomial expansions of the FF(64,0.78) force field yield slightly different results (typically the self-diffusion is slightly lower). The length of the simulation to determine the self-diffusion coefficients was 2 ns for all temperatures. The self-diffusion coefficients were deduced from the slope of the linear part of the mean-square displacement (MSD). The MSD was calculated by averaging squared displacements of F atoms over all trajectory samples with a time length of integer multiples of 2 ps and a maximal length of 1 ns. The time step in the sampling was 1 ps. The accuracy in estimates of the self-diffusion will be better at higher temperatures. Due to the fact that self-diffusion in liquid HF is often underlaid by complex and uneven dynamics involving parallel and rotating zigzag chains,62 an accurate determination of the self-diffusion coefficients, especially at low temperatures, may require much longer than a few nanoseconds MD simulation.

However, the agreement with the experiment is reasonable at T ) 300 K. There is also discrepancy among different experimental results on the self-diffusion in the liquid HF. The first published results were by O’Reilly,63 which were regarded as reliable up to (15%. More recent data64 corrected the selfdiffusion to lower values with error bars of (5%. The best exponential fits to both data are summarized in Table 7. 5. Vibrational Spectrum. Figure 12 shows a comparison of the power spectra of liquid HF from both the CP MD simulation and the MD simulations using the FF(64,0.78) model. The power spectra were obtained via Fourier transformation of the respective velocity-velocity autocorrelation functions. The F-H stretching band in the spectrum from the CP MD simulation is seen to be broader and centered at approximately 3000 cm-1. Narrowness of the F-H stretching band in the power spectrum of the FF(64,0.78) model is a result of the simple harmonic approximation to the intramolecular FH potential. However, as shown in ref 36 for water, a combination of a more sophisticated intramolecular potential with the intermolecular part of the force-matching force field can significantly improve the agreement of the O-H stretching band shapes between the CP and force-matching MD simulations, but we have not adopted such a procedure in the work. For HF, the anharmonicity of the intramolecular potential is even stronger.23 The band around 500-600 cm-1 originates from the intrachain vibrational motions of hydrogen in the F-H‚‚‚F stretching mode.8 This band as predicted by the FF(64,0.78) model at T ) 300 K is broader compared to the CP MD simulation. On the same figure is plotted the power spectrum from the same force-matching MD simulation, which was used to obtain the RDFs shown in Figure 6. In this simulation, the system has a similar degree of equilibration as in the CP MD simulation. The agreement with the CP MD is seen to be much improved, especially in the region of the first slope of the band. The upper slope of the band in the CP MD spectrum has a submaximum that is clearly absent in the band from the FF(64,0.78) simulation. However, this discrepancy could again be due to the harmonic approximation of the F-H intramolecular potential. 6. Thermodynamic Properties. In Table 8, some properties of simulated liquid HF are compared at selected temperatures. Shown are the properties for different representations of the FF(64,0.78) short-ranged interaction. Because the MD simulations in the NPT ensemble in the supercritical region (states II-VI) are characterized by slow and large fluctuations in the system volume and pressure, it is difficult to determine the thermodynamic properties of these states accurately.

TABLE 8: Thermodynamic Properties for a Different Representation of the FF(64,0.78) Modelc T ) 203 K p ) 1 bar property

spl

Chb

rvp

F (kg/m3) Upot (kJ/mol) Uintra (kJ/mol) UCoul (kJ/mol) Usr (kJ/mol) CP (J/(mol K)) RP (10-5 1/K)

1174 -24.6 1.8 -27.4 1.0 51 210

1186 -24.7 1.9 -27.8 1.0

1169 -24.5 1.9 -27.7 1.3

T ) 273 K p ) 1 bar expt 1176a -31.9

56.48b 230b

T ) 300 K p ) 2 bar (state I)

spl

Chb

rvp

expt

spl

Chb

rvp

expt

1022 -21.1 2.0 -24.4 1.3 65 231

1038 -21.4 2.3 -25.1 1.4

1033 -21.1 2.0 -24.5 1.4

1015a -29.01

949 -19.6 2.1 -23.0 1.3

948 -19.7 2.1 -23.1 1.3

949 -19.6 2.1 -23.1 1.4

962a -27.7

71.13b 225b

a References 60, 61, and 48. b Reference 15. c The energy properties were obtained from MD simulations in the NVT ensemble and experimental densities. The densities were obtained from the MD simulations in the NPT ensemble. The “spl”, “Chb”, and “rvp” denote representations of short-ranged forces by, respectively, spline, expansion as Chebyshev series, and expansion as powers of 1/r. The notation is as follows: T, temperature of the simulation; p, pressure; F, density; Uintra, intramolecular energy; Upot, configurational energy; UCoul, Coulomb energy; Usr, short-ranged interaction energy; CP, isobaric heat capacity; RP: thermal expansion coefficient.

Force Field for Liquid HF from MD Simulation A first observation is that the densities are reproduced well in the whole temperature range, with an error to within 2%. The average potential energy, Upot, is higher compared to experiment. It is interesting that at T ) 300 K (state I), which is nearly the boiling point (Table 1), the Upot is very close to the dimerization energy in a vacuum (Table 5). The disagreement in Upot with experiment is systematic. For example, if one corrects the simulated Upot by its difference with the experimental value at T ) 203 K, for example, 7.3 kJ/mol, then one will have Upot values at T ) 273 K and T ) 300 K of -28.4 and -26.9 kJ/mol, respectively. The systematic deviation in the simulated energetic properties from the experiment, together with the fact that the densities are reproduced accurately, suggests that the properties that rely on the derivatives of the system internal energy (e.g., isobaric and isochoric heat capacities CP and CV) may be in better agreement with experiment. As seen from Table 8, this is indeed the case. The constant pressure heat capacity CP ) (∂H/∂T)P and thermal expansion coefficient RP ) 1/V (∂V/∂T)P, where H and V are the enthalpy and volume of the system, were calculated using numerical differentiation. For that purpose, the MD simulations in the NPT ensemble were carried out around the temperature of interest with a temperature step of 10 K. These quantities from the MD simulation were not calculated using polynomial representations of the FF(64,0.78) model. However, judging from the density and average configurational energy using the different representations, these properties should be not much different. The expansion of the short-ranged force as powers of 1/r yields almost the same thermodynamic properties compared to the spline representation, despite the fact that the structural properties are noticeably different (especially in gFF). IV. Conclusions In this paper, we have applied a newly developed method for force matching to generate an effective two-site force field for HF using trajectory and force data obtained from a CarParrinello MD simulation. This simulation was for a system of 64 HF molecules at nearly ambient conditions using the BLYP energy functional and a fictitious electronic mass of 340 au. The fictitious electronic mass in the CP MD simulation was chosen as suggested recently for hydrogen-bonded systems.32 The force-matching procedure naturally includes a fit of partial charges and intramolecular forces. The short-ranged interaction noticeably drops to zero around 0.4 nm. However, we have used a cutoff of 0.78 nm to the short-ranged interaction in order to yield correct thermodynamic properties. The resulting model has been referred to as the FF(64,0.78) model. Several representations by polynomial expansions of the short-ranged part of the FF(64,0.78) model were defined, which can be used for a more efficient implementation of the model. The performance of the FF(64,0.78) force field was first tested on the gas-phase HF dimer. The model is able to reproduce well the energetics and the (bent) structure of the dimer, despite the fact that the force field does not account for the molecular quadrupole moment. However, the hydrogen-bond length within the dimer is very close to the position of the first maximum in the condensed-phase gFH and, thus, deviates from the experiment. Generally, the bent structure of the dimer more closely resembles the arrangement of the monomers inside the zigzag chains in the liquid that is apparently a consequence of fitting the force field to the HF liquid phase. The dipole moment calculated using the force-matching partial charges is very close to the dipole moment of the monomer in a vacuum. The force

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6585 constant evaluated from the intramolecular part of the FF(64,0.78) force field agrees very closely with a best fit of Raman and IR spectra of liquid HF at ambient conditions.8 The structural properties of liquid HF at different state conditions obtained from the force-matching model are in satisfactory agreement with the available experimental data. The overstructuring of the RDFs in the CP MD simulation compared to the MD simulation using the force-matching model and to experiment can be explained by the fact that a much longer equilibration run is needed in the ab initio MD. This fact was demonstrated by presenting the RDF evolution over 10 ps time segments in the MD simulation with force-matching model. An analysis of the aggregates shows the formation of zigzag chains with an average length of four to five monomers at conditions close to the boiling point, while about seven to eight monomers at temperatures closer to the freezing point. This behavior also agrees well with the available interpretations of the experimental data. The intrachain orientations of the monomers resemble an average arrangement of monomers inside chains, as observed experimentally in solid DF. The FF(64,0.78) model shows an ability to reproduce the density of liquid HF over a broad range of temperatures (at least 203-300 K) with an error to within 2% of the experimental data. The calculated internal potential energy, Upot, of liquid HF at different temperatures is systematically higher than the experimental values by about 7 kJ/mol. Interestingly, the Upot at T ) 300 K, which corresponds to conditions close to the boiling point, nearly coincides with the dimerization energy in a vacuum. The systematic deviation of Upot from the experimental result, together with the fact that densities are reproduced correctly, results in thermodynamic properties relying on the derivatives of Upot and the system volume (e.g., CP) agreeing well with experiment. The simple two-site force-matching model presented in this paper seems able to reliably predict the behavior of liquid HF. This result is another demonstration of the utility of the newly developed force-matching method for obtaining effective force fields for condensed-matter systems from ab initio MD simulations. Acknowledgment. This research was supported by the United States Office of Naval Research. The computations were supported in part by the United States National Science Foundation (NSF) Cooperative Agreement No. ACI-9619020 for computing resources provided by the National Partnership for an Advanced Computational Infrastructure at the San Diego Supercomputer Center. An allocation of computer time from the Center for High Performance Computing at the University of Utah is gratefully acknowledged. Supporting Information Available: The spline parameters of the short-ranged part of the FF(64,0.78) force field. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) (2) 4942. (3) (4) P. Mol. (5) 1998. (6) (7) (8)

Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1979, 71, 298. Jorgensen, W. L.; Cournoyer, M. E. J. Am. Chem. Soc. 1978, 100, Jorgensen, W. L. J. Am. Chem. Soc. 1978, 100, 7824. Deraman, M.; Dore, J. C.; Powles, J. G.; Holloway, J. H.; Chieux, Phys. 1985, 55, 1351. Johnson, M. W.; Sa´ndor, E.; Arzi, E. Acta Crystallogr. 1975, B31, Cole, R. H. J. Chem. Phys. 1973, 59, 1545. Sheft, I.; Perkins, A. J. J. Inorg. Nucl. Chem. 1976, 38, 665. Desbat, B.; Huong, P. V. J. Chem. Phys. 1983, 78, 6377.

6586 J. Phys. Chem. B, Vol. 109, No. 14, 2005 (9) Ro¨thlisberger, U.; Parrinello, M. J. Chem. Phys. 1997, 106, 4658. (10) Janzen, J.; Bartell, L. S. J. Chem. Phys. 1969, 50, 3611. (11) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (12) Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381. (13) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (14) Klein, M. L.; McDonald, I. R.; O’Shea, S. F. J. Chem. Phys. 1978, 69, 63. (15) Cournoyer, M. E.; Jorgensen, W. L. Mol. Phys. 1984, 51, 119. (16) Jedlovszky, P.; Vallauri, R. Mol. Phys. 1997, 92, 331. (17) Jedlovszky, P.; Vallauri, R. J. Chem. Phys. 1997, 107, 10166. (18) Jedlovszky, P.; Mezei, M.; Vallauri, R. J. Chem. Phys. 2001, 115, 9883. (19) Della Valle, R. G.; Gazzillo, D. Phys. ReV. B 1999, 59, 13699. (20) Martin, C.; Lombardero, M.; Anta, J. A.; Lomba, E. J. Chem. Phys. 2001, 114, 355. (21) Lie, G. C.; Clementi, E. Phys. ReV. A 1986, 33, 2679. (22) Rick, S. W.; Stuart S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (23) Klopper, W.; Quack, M.; Suhm, M. A. J. Chem. Phys. 1998, 108, 10096. (24) Sum, A. K.; Sandler, S. I.; Naicker, P. K. Fluid Phase Equilib. 2002, 199, 5. (25) Wierzchowski, S. J.; Kofke, D. A. J. Chem. Phys. 2003, 119, 6092. (26) Wierzchowski, S. J.; Kofke, D. A.; Gao, J. J. Chem. Phys. 2003, 119, 7365. (27) Liem, S. Y.; Popelier, P. L. A. J. Chem. Phys. 2003, 119, 4560. (28) Hernandes, M. Z.; da Silva, J. B. P.; Longo, R. L. J. Comput. Chem. 2003, 24, 973. (29) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (30) Kreitmeir, M.; Bertagnolli, H.; Mortensen, J. J.; Parrinello, M. J. Chem. Phys. 2003, 118, 3639. (31) Raugei, S.; Klein, M. J. Am. Chem. Soc. 2003, 125, 8992. (32) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 120, 300. (33) Tangney, P.; Scandolo, S. J. Chem. Phys. 2002, 116, 14. (34) Ercolessi, F.; Adams, J. Europhys. Lett. 1994, 26, 583. (35) Tangney, P.; Scandolo, S. J. Chem. Phys. 2002, 117, 8898. (36) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. J. Chem. Phys. 2004, 120, 10896. (37) Lawson, C. L.; Hanson, R. J. SolVing Least Squares Problems; Prentice Hall: Englewood Cliffs, NJ, 1974. (38) De Boor, C. A Practical Guide to Splines; Springer-Verlag: New York, 1978. (39) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133.

Izvekov and Voth (40) Becke, A. D. Phys. ReV. A 1988, 38, 3098. Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (41) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (42) Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1993, 99, 9080. (43) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 105, 1142. (44) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 3572. (45) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2002, 116, 10372. (46) Hutter, J.; Ballone, P.; Bernasconi, M.; Focher, P.; Fois, E.; Goedecker, S.; Marx, D.; Parrinello, M.; Tuckerman, M. CPMD, version 3.5; MPI fu¨r Festko¨rperforschung: Stuttgart, 1997-2001. (47) Marx, D.; Hutter, J. In Modern Methods and Algorithms of Quantum Chemistry; John von Neumann Institute for Computing: FZ Ju¨lich, 2000; pp 301-449. (48) Pfleiderer, T.; Waldner, I.; Bertagnolli, H.; To¨dheide, K.; Fischer, H. E. J. Chem. Phys. 2000, 113, 3690. (49) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 1996. (50) Rivlin, T. J. ChebysheV Polynomials; Wiley: New York, 1990. (51) Forester, T. R.; Smith, W. DL_POLY User Manual; CCLRC, Daresbury Laboratory: Daresbury, Warrington, U.K., 1995. (52) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure, IV, Constants of Diatomics Molecules; Van Nostrand Reinhold: New York, 1979. (53) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids; Clarendon Press: Oxford, 1984. (54) Howard, B. J.; Dyke, T. R.; Klemperer, W. J. Chem. Phys. 1984, 81, 5417. (55) Pine, A. S.; Howard, B. J. J. Chem. Phys. 1986, 84, 590. (56) Kofranek, M.; Lischka, H.; Karpfen, A. Chem. Phys. 1988, 121, 137. (57) Klopper, W.; Quack, M.; Suhm, M. A. J. Chem. Phys. 1998, 108, 10096. (58) McLain, S. E.; Benmore, C. J.; Siewenie, J. E.; Urquidi, J.; Turner, J. F. C. Angew. Chem., Int. Ed. 2004, 43, 1952. (59) McLain, S. E.; Benmore, C. J.; Siewenie, J. E.; Molaison, J. J.; Turner, J. F. C. J. Chem. Phys. 2004, 121, 6448. (60) Simons, J. H.; Bouknight, J. W. J. Am. Chem. Soc. 1932, 54, 129. (61) Sheft, I.; Perkins, A. J.; Hyman, H. H. J. Inorg. Nucl. Chem. 1973, 35, 3677. (62) Ring, J. W.; Egelstaff, P. A. J. Chem. Phys. 1970, 52, 5973. (63) O’Reilly, D. E. J. Chem. Phys. 1970, 52, 5974. (64) Karger, N.; Vardag, T.; Lu¨demann, H.-D. J. Chem. Phys. 1990, 93, 3437.