Fast hybrid methods for the simulation of dielectric constants of liquids

Journal of Molecular Liquids 117 (2005) 3 – 16 www.elsevier.com/locate/molliq

Fast hybrid methods for the simulation of dielectric constants of liquids J. Richardia,*,1, P.H. Friesa,*, C. Millotb a

Laboratoire de Reconnaissance Ionique, Service de Chimie Inorganique et Biologique, CEA/DSM/De´partement de Recherche Fondamentale sur la Matie´re Condense´e, CEA-Grenoble F-38054 Grenoble cedex 9, France b Universite´ Henri Poincare´—Nancy 1, Laboratoire de chimie the´orique UMR CNRS-UHP no. 7565 BP 239, 54506 Vandoeuvre-les-Nancy Cedex, France Available online 29 September 2004

Abstract The distance-dependent Kirkwood factors and dielectric constants of the SPC/E and TIP3P models of water, the H1 model of methanol, and the OPLS model of acetonitrile are investigated with the help of long molecular dynamics (MD) simulation runs and the hypernetted chain (HNC) approximation of the molecular Ornstein–Zernike (MOZ) theory. The MOZ theory predicts dielectric constants which show the same dependence on the liquid model as their simulated counterparts, are very accurate in the case of acetonitrile, but are too low by about 20% for the H-bonding liquids, because the HNC closure moderately accounts for the shortrange orientational order of their molecules. In the case of too short simulation runs, some pitfalls concerning the accuracy of the dielectric properties and their errors are emphasized. For a given liquid model, if the distance-dependent Kirkwood factor of the MOZ theory has a long-range profile similar to the converged one of an accurate simulation, rapid and reliable hybrid methods based on the combination of the MOZ and simulation results can be proposed to estimate the dielectric constant. Test calculations on the abovementioned models are carried out to assess the validity of the approach. Hybrid methods are shown to give dielectric constants of similar accuracy as those derived from long simulation runs of one million of steps or more, but with a CPU time demand at least 10 times shorter. D 2004 Elsevier B.V. All rights reserved. Keywords: Fast hybrid methods; Simulation; Dielectric constants

1. Introduction The models of intermolecular interactions, which have been developed during the last two decades, are nowadays widely used to calculate condensed phase properties. They are applied to the study of chemical and biological processes at a molecular level. A significant test of the quality of an interaction model is to compare derived theoretical liquid properties with their experimental counterparts. Such a comparison was often restricted to the site–site distribution functions, the heats of vaporization, the densities, and the

* Corresponding authors. Fax: +33 4 38 78 5090. E-mail addresses: [email protected] (J. Richardi)8 [email protected] (P.H. Fries). 1 Present address: Laboratoire des Mate´riaux Me´soscopiques et Nanome´triques, UMR-CNRS 7070, Universite´ Pierre et Marie Curie (Paris VI), B.P. 52, 4, place Jussieu, 75230 Paris Cedex 05, France. 0167-7322/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2004.08.016

diffusion coefficients. On the contrary, the theoretical macroscopic dielectric constant was rarely determined by simulations and compared to the experimental value. This is in contrast with the theoretical and practical importance of this property, which accounts for the local orientational order of the molecules and its propagation through space, and is crucial for the response of the solvent to the charge distribution of a solute. When insulators are placed between charged conductors, changes of the electric potentials or charges of these conductors are observed [1]. From a phenomenological point of view, these effects can be accounted for by simply assuming that inside the insulators, the electric field gives rise to a density of electric dipoles, named polarization [1,2]. The polarization P is proportional to the averaged electric field or macroscopic field E inside matter, i.e. P=[(e s
1)/ 4p]E, where e s is the macroscopic dielectric constant. The rigorous molecular theory of e s came up for a long time

4

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

against many conceptual and computational difficulties which were summarized at the beginning of the 1970s by Bfttcher [2], who gave a comprehensive presentation of the various approximate descriptions available at that time. Shortly after, the development of new formalisms of the statistical mechanics of molecular fluids together with substantial progresses in the proper handling of the longrange electrostatic interactions in computer simulation have led to a more and more thorough understanding of the relationships between the dielectric constant and the molecular correlation functions. The major conceptual problems were resolved at the beginning of the 1980s [3– 6]. Thus, various simple expressions of e s in terms of the liquid structure had been derived in the absence of applied electric field by calculating how the electric field due to a given molecule affects the surrounding molecules [3]. Besides, using a consistent molecular theory of the timedependent polarization induced in a liquid by a weak, spatially and temporally varying applied electric field, expressions of the frequency dependent dielectric constant had been obtained as a function of the structural and microdynamic properties of the liquid [4]. The potentiality of this approach had also been successfully tested in the case of liquid acetonitrile [7]. Unfortunately, a systematic study of the dielectric properties of fluids is still hampered by too heavy computational demands. Usually, the dielectric constant is derived from the Kirkwood factor defined PN Yby the fluctuation formula hM2i/(Nl 2), where M ¼ i¼1 l i is the total dipole moment of the simulation cell. Y l i is the dipole of the molecule i and N denotes the number of molecules of the sample. In comparison with other properties, hM2i converges very slowly during a simulation. Moreover, hM2i is very sensitive to the boundary conditions chosen in order to handle the long-range interactions. In general, the lattice sum techniques (e.g., Ewald sum or Ladd method) or a reaction field method must be used to obtain a consistent dielectric constant. Therefore, even on a modern workstation, the accurate determination of the dielectric constant of simple water model, such as the SPC/E model, takes several weeks. This should be compared with the few hours required to obtain well-stabilized site–site distribution functions, energies, and densities for the same model. Finally, in a ferrofluid, i.e. a colloidal suspension of ferromagnetic particles, the computation of the magnetic susceptibility faces similar problems when the concentration of the particles and their magnetic dipole–dipole interactions become strong enough to give rise to orientational correlations between the magnetic dipoles [8]. Besides the fluctuation formula, other methods were envisaged to facilitate the simulation of dielectric properties. They include the use of the ion–ion potentials of mean force [9] and of functions which characterize the dipole–dipole orientational correlations [10]. The dielectric constant can be extracted from the long-range dependence of these functions, but the corresponding simulations require rather large systems leading to a strong reduction of the efficiency

of this method. Several workers studied the polarization of the system as a function of an applied electric field [11–13]. The applied field method was claimed to be two to four times more efficient than the usual fluctuation method [12], though this can be called in question by another study [13]. Variations of these methods are the transient response methods, in which the time-dependent polarization response is monitored [14]. Another technique is based on the calculation of the wave vector-dependent permittivity tensor e˜˜ s ðkÞ [5,7,15,16], which is computed for several of the smallest wave vector lengths k accessible by simulation. Then, the dielectric constant is derived from the kY0 limit of the transverse components of e˜˜ s ðkÞ. Other methods are based on the calculation of the probability distribution function of the total dipole moment ( P(M)) or of its square ( P(M 2)). Kusalik et al. [17] showed that this approach becomes five times more efficient than the fluctuation formula when it is combined with a functional form for P(M 2). It was also found that sampling P(M) by the umbrella method can improve the efficiency of the simulation of dielectric constants [18]. In addition to the simulation methods, integral equations of the statistical mechanics of fluids [19] can be used to compute the pair distribution function from a given intermolecular pair potential. In the early 1970s, for rigid anisotropic particles, Blum and Torruella [20] pioneered the convenient expansions of the two-body spatial correlation functions in rotational invariants. Their formalism is particularly suitable to include the simplifications of molecular symmetry and leads to a tractable solution of the molecular Ornstein–Zernike (MOZ) relation [21], which is central in the integral equation methods. It applies to polar molecules and one of its first successful applications to solvents of practical interest was the calculation of the dielectric constant of a model of water by Carnie and Patey [22] who showed the importance of the electric quadrupole and molecular polarizability. The MOZ theory was then solved in conjunction with closures of the hypernetted chain (HNC) type for dipolar hard spheres [23–25]. It was shown that, when compared with bexactQ simulation data, these HNC type approximations give good or even excellent results for configurational energies, dielectric constants, and other studied structural properties. This preliminary success was the starting point of the considerable efforts which in the nineties have proved the accuracy and usefulness of the MOZ–HNC type approximations for fluids of particles showing a large variety of anisotropic shapes and charge distributions. Indeed, for models of molecular fluids including homonuclear diatomic systems [26,27], hydrogen chloride [28], SO2 and H2S [29,30], acetonitrile [31,32], acetone and chloroform [33], and solutions of ions in acetonitrile and acetone [34,35], the integral equations produce structural and thermodynamic properties, which are in good agreement with the simulation data. In particular, the simple HNC approximation is accurate for those properties, such as the dielectric constants, which are

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

rather insensitive to the details of the repulsive part of the intermolecular potential. Furthermore, if the molecular polarizability is taken into account by a self-consistent mean field (SCMF) procedure, the HNC approximation gives dielectric constants which for simple molecular models show the right experimental trends for such different polar solvents as acetonitrile, acetone, and chloroform [36,37], formamide, NMF and DMF [38], THF and methylene chloride [39]. When the intermolecular electrostatic interactions are accurately described, the theoretical HNC predictions often become remarkable as for acetonitrile and acetone [40]. Finally, the HNC theory can serve to investigate the important and difficult question of the impact of molecular shape and polarity on chiral discrimination in fluids [41]. The HNC potential of mean force between a probe solute and a complex of paramagnetic gadolinium Gd3+ in water is also useful for an accurate description of the relative motion of the two solutes [42]. This leads to an unexpected application to the relaxation contrast agents in medical magnetic resonance imaging (MRI): The magnetic dipole–dipole Hamiltonian, which couples nuclear spins on the probe solute with the electronic spin of the Gd3+ complex, is modulated by the intermolecular motion of these two solutes, and this induces a paramagnetic relaxation enhancement of the nuclear spins. The measurement of this enhancement allows an indirect determination of the Gd3+ electronic relaxation, which is one of the factors affecting the efficiency of the contrast agents. Unfortunately, recent studies of models of water and methanol [43,44] showed that at ambient conditions, the HNC dielectric constants quantitatively differ from their simulation counterparts. This probably occurs because the strong local orientational ordering of the molecules caused by the H-bonds at room temperature cannot be accounted for by the spherically symmetric bridge functions used in the chosen HNC or HNC type approximations. However, this drawback of the HNC theory disappears together with the H-bond network at higher temperatures, and in all cases, the HNC approximation still predicts the main structural features and trends. This discussion clearly infers that the MOZ–HNC approximation gives at least a correct semiquantitative picture of the liquid structure of polar solvents. The above discussion reminds us that computing the macroscopic dielectric constant e s, which is of central interest for polar solvents, is still one of the most difficult tasks of the statistical mechanics of liquids. To sum up, numerical simulations demand very long computer times in order to avoid discouraging statistical scatter and MOZ theories give results which are sometimes only semiquantitative, for example for liquids of associating molecules. In Section 2, the employed molecular dynamics (MD) techniques, the MOZ theory, and the respective ways to calculate the Kirkwood factor g K yielding e s are presented. Section 3 deals with interaction models of water (SPC/E and TIP3P), methanol (H1), and acetonitrile (OPLS). The simulated dielectric constants are compared to the accurate

5

values obtained previously for the same water models [45]. The statistical convergence of the simulated short and longrange contributions of the dipole–dipole correlations to g K is carefully analyzed. For simulation runs of 100 000 steps, it will be shown that the short-range contribution is properly converged, but that the considerable scatter, which can be displayed by the long-range contribution, prevents the accurate computation of the Kirkwood factor and of e s. By performing long simulation runs of 300 000 to 2 000 000 steps, the bexactQ long-range dipole–dipole contributions to g K will be computed for various liquids and shown to be small quantities, of similar sizes to those of their MOZ counterparts. Section 4 stems from the previous observations and introduces a new class of efficient methods for calculating the Kirkwood factor g K. These methods are hybrid since they combine the simulation which rapidly gives the converged short-range dipole–dipole contribution to g K, together with the MOZ–HNC theory which provides a reasonable approximation of the small long-range contribution. The reliability and accuracy of the hybrid methods are checked by comparison with the precise simulation data of the long runs performed for this work. The parameters of the intermolecular potentials were proposed in the literature and correspond to the SPC/E [46] and TIP3P [47] model of water, the H1 model of methanol [48], and the OPLS model of acetonitrile [49]. The dielectric constants are calculated at 258 C using the experimental densities (0.997 g/cm3 for water [47], 0.787 g/cm3 for methanol [16], and 0.777 g/cm3 for acetonitrile [49].

2. Methods 2.1. Molecular dynamics simulations In order to study the statistical accuracy of the Kirkwood factor, very long MD simulations were carried out. The equations of motion of the NVE ensemble were integrated numerically for short time periods. The time step was 1.5 fs for water and 2.0 fs for methanol and acetonitrile. We employed a crude method for performing constant temperature MD [50] by simply recalibrating the velocities every 100 steps, so as to fix the kinetic temperature of the system. The rotation matrix defining the orientation of each rigid molecule [50] was represented by a quaternion [51,52], the components of which satisfy equations of motion which were solved using a leap-frog evolution algorithm introduced by Fincham [50]. To obtain reliable dielectric constants, the long-range electrostatic potential should be correctly handled. For Stockmayer fluids, it was shown that the dielectric constant obtained with the reaction field method can exhibit a strong dependence on the size of the simulation cell even with more than 512 molecules [17]. The dependence is much smaller when lattice sum techniques are used. Therefore, the long-range electrostatic interaction was taken into account

6

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

by truncating the Ladd expansion after the terms of the 10th order [53,54], which gave equivalent results to the Ewald method for Stockmayer liquids [55]. The Ladd approach is a systematic expansion of the lattice sums leading to a modified expression of the electrostatic interaction tensors. Details are given elsewhere [56]. Altogether, nine simulation runs were carried out using the Ladd method. They differed in the liquid model, in the length of the simulation, in the number of molecules and in the dielectric constant of the reaction field. In most cases, a dielectric constant of the reaction field close to that found for the studied model was chosen. Additional runs were carried out with the reaction field method in order to get only the short-range orientational correlation function needed by the hybrid methods. Therefore, as explained later, these runs can be limited to rather small numbers of steps. In Table 1, the parameters of the simulation runs are given. For acetonitrile, we did not use numbers of simulation steps as large as for the other liquid models, because the statistical errors of the dielectric constants observed for these simulation runs of moderate lengths are already small (see Table 2). 2.2. The molecular Ornstein–Zernike theory At the beginning of the 1970s, Blum and Torruella [20,21] pioneered an efficient rotational invariant formalism to solve the MOZ theory. Within this context, Carlevaro et al. [57] developed the analytical approach for fluids of molecules coupled by more and more sophisticated potentials, which they designed in order to capture the essential physical features while keeping sufficiently simple for an analytical solution. However, only thanks to the computer technology progress of the nineties, it became possible to solve the MOZ theory with the accurate HNC approximation for realistic molecular models similar to those employed in the simulations. Details of the method are given elsewhere [20,21,23,29,33,36,37]. Table 1 Parameters of the MD simulation runs: the number of simulation steps, the number of molecules in the simulation cell, the method applied to treat the long-range electrostatic interaction and the dielectric constant e ext of the external surrounding of the simulation box and its replicas Run

Liquid model

Steps

Molecules

Method

e ext

1a 1b 1c 1d 2a 2b 2c 3a 3b 3c 4a 4b 4c

water (SPC/E) water (SPC/E) water (SPC/E) water (SPC/E) water (TIP3P) water (TIP3P) water (TIP3P) methanol (H1) methanol (H1) methanol (H1) acetonitrile (OPLS) acetonitrile (OPLS) acetonitrile (OPLS)

2 000 000 2 000 000 300 000 300 000 1 000 000 1 000 000 100 000 1 000 000 300 000 100 000 300 000 300 000 100 000

512 256 256 256 512 256 256 512 256 256 512 256 256

Ladd Ladd Ladd RF Ladd Ladd RF Ladd Ladd RF Ladd Ladd RF

70 70 l l 92 92 l 24 24 l 18 24 l

In comparison with computer simulations, the MOZ theory is a very fast method. Typically, the calculation of the properties of the water models studied here takes less than an hour on a usual workstation. However, because of the HNC approximate closure, which is used to solve the MOZ theory, this method can only give approximate liquid properties. Recent comparisons between the MOZ and simulation results indicate the accuracy of the dielectric constants, which we want to calculate here. For aprotic solvents such as acetonitrile [32], acetone [33] and DMF [38], the dielectric constants obtained by the MOZ(HNC) theory and by simulations agree within the error of the simulations. In contrast, for H-bonding liquids such as water and methanol [43], the MOZ theory underestimates dielectric constants by about 20% with respect to the simulation values. This discrepancy is too large to hold the MOZ(HNC) approximation for an accurate alternative to simulations for computing dielectric constants of H-bonding liquids. 2.3. The evaluation of the dielectric constants The macroscopic dielectric constant e s can be calculated from the Kirkwood formula: ðes
1Þð2es þ 1Þ ¼ y D gK ; 9es

ð1Þ

where y D=4pql 2/(9k BT) is the dipole parameter and g K the Kirkwood factor. l is the absolute value of the molecular dipole, q the number density, and T the temperature. In a simulation, an intermediate Kirkwood factor g KV can be computed from the formula g VK ¼

X hM2 i Y ; with M ¼ l i: N l2 i

ð2Þ

In Eq. (2), h. . .i indicates the average over the finite simulation box. g KV depends on the value of the dielectric constant e ext of the external dielectric continuum surrounding the simulation box and its replicas. Then, it must be corrected to obtain the Kirkwood factor: gK ¼

ð2eext þ es Þð2es þ 1Þ ¼ g VK ð2eext þ 1Þ3es

ð3Þ

When e ext=e s, no correction is required. We have g K=g VK and e s can be calculated from Eq. (1). In the general case, since e s is unknown, we have e ext pe s. Then, substituting the expression (3) for g K in the Kirkwood formula (1), we obtain the modified Kirkwood equation adapted to simulations ðes
1Þð2eext þ 1Þ ¼ 3yD g VK : ð2eext þ es Þ

ð4Þ

Within the framework of the MOZ theory, g K is obtained as follows. We introduce the R-dependent Kirkwood factor

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

7

Table 2 Dielectric constants of models of liquid water (SPC/E and TIP3P), methanol (model H1), and acetonitrile (OPLS) derived from the simulations, MOZ theory, and MD/MOZ hybrid methods for different lengths of the simulation runs Liquid model

Run

Water (SPC/E)a

Water (TIP3P)c

Methanol (H1)e

Acetonitrile (OPLS)f

1a 1b 1c 1d 1db 2a 2b 2cb,d 3a 3b 3cb,d 4a 4b 4cb

MD/MOZ hybrid

Simulations

Truncation

Tail-correction

72.6F1.0 72.0F1.4 74.0F4.5 59.1F2.6 62.5F0.4 91.4F1.3 91.4F1.3 67.1F0.9 24.8F0.7 24.2F1.0 23.0F1.8 20.3F1.0 20.3F1.0 26.2F0.8

71.6F1.0 71.0F1.4 73.0F4.5 58.1F2.6 65.1F0.4 96.3F1.3 96.3F1.3 85.0F0.9 24.5F0.7 23.8F1.0 23.1F1.8 18.9F1.0 19.0F1.0 20.4F0.8

MOZ theory

65.7F6.2 74.4F7.5 65.7F12 88.1F25

58.2

96.4F14 104.6F8.2 156F125 25.4F3.9 20.1F2.5 55F143 17.8F2.1 17.8F2.1

79.1

19.2

18.0

The errors were calculated from the variances according to Eqs. (11) and (12), using small runs of l sr = 20 000 steps and f = 2. a Dielectric constants for the SPC/E model in the literature: 70.7F0.8 [58], 69.6F1.5 [59], 68.2F5.8 [45]. The first two errors given in the literature are very small compared with the limited lengths of the performed simulations. For instance, in the case of the second value reported by Svishchev and Kusalik [59], three independent runs of 50 000 steps were performed, giving values of 71.0F6.1, 68.5F6.6, and 69.2F7.2. The three resultspwere combined and ffiffiffi predicted a dielectric constant of 69.6F1.5. Due to the statistical laws, the statistical error of a three times longer run is expected to be 1= 3 times the errors of the single runs. This corresponds to an error of about 4.0. b Use of the first maximum of G K as the truncation distance. c Dielectric constants for the TIP3P model in the literature: 102F3 [60], 96.9F6.7 [45]. d For TIP3P water and OPLS acetonitrile, the short runs (2c) and (3c) of 100 000 steps clearly lead to unphysical dielectric constants and unreasonable statistical errors. e Dielectric constants for the H1 model of methanol in the literature: 24F2 [16]. f Dielectric constants for the OPLS model of acetonitrile in the literature: 18F1 [61].

G K(R), which can be written in terms of a weighted integral of the rotational invariant coefficient g 110 00 (R) as Z R 4p 110 g00 ð R VÞR V2 dR V; ð5Þ GK ð RÞ ¼ 1
pffiffiffi q 3 0 where R (or RV) is the distance between the centers of two molecules. The Kirkwood factor is simply the long-distance limit of G K(R), i.e. g K=G K(R=l). Note that the coefficient 110 g 00 (R) is a direct result of the procedure employed to solve the MOZ theory, using the rotational invariant expansions of the pair correlation functions. The rotational invariant coefficient g110 00 (R) can also be computed from the simulation trajectories as: pffiffiffi 3 110 g00 ð RÞ ¼
qð N
1ÞdV + * N N X X lˆ i lˆ j ; ð6Þ i¼1 j¼1; j p i R
dRVRij VR þ dR

al. [45], who used similar statistics. For similar statistics, note that our error estimate of e s is larger by a factor of 2. Indeed, it was derived from an error of g VK , proportional to twice the square root of the variance of g VK , and not just the square root of this variance as usual [50] (see Section 3.2). The MOZ dielectric constant, also reported in Table 2, shows the same dependence on the liquid model as its simulated counterpart. The value for acetonitrile is excellent, but too low by about 20% for each model of H-bonding liquid, because the HNC closure moderately accounts for the short-range orientational order of its molecules [43,44]. Thus, the long simulation runs of this work confirm previous observation.

3. Statistical analysis of the R-dependent Kirkwood factor G K(R) obtained by simulations 3.1. The slow statistical convergence of g KV

3

where N is the number of molecules, dV=(4p/3)[(R+dR)
(R
dR)3], and lˆ i the unit vector of the total dipole of molecule i. Our simulated dielectric constants are given in Table 2. For the long runs (1a), (1b) and (2a), (2b) concerning the SPC/E and TIP3P water models, they are in excellent agreement with the accurate values obtained by Hfchtl et

Let us try to understand why the dielectric constant cannot be reliably computed from short simulation runs of 20 000 to 100 000 steps. For that purpose, the long simulation runs for the TIP3P model of water (run (2a)) and for the acetonitrile model (run (4a)) were splitted into small runs of 100 000 simulation steps. For each one of

8

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

these runs, we computed the R-dependent Kirkwood factor as + * N X N X 1 Y Y li d lj ; GK ð RÞ ¼ 2 ð7Þ l N i
1 j¼1 Rij VR

where the summation index i runs over all the N molecules in the simulation box whereas j designates those molecules which satisfy the minimum image convention with i and have an intercentre distance R ij VR with i. In particular, j can be equal to i. Let l box be the side length pffiffiffi of the simulation box. For a simulation, if we take R ¼ 3lbox =2 the summation index j runs over all the molecules in the box. Therefore, the value pffiffiffi of G K(R) at R ¼ 3lbox =2 simply corresponds to the factor g VK , which yields the Kirkwood factor g K according to Eq. (3). In each simulation using a finite value of the dielectric constant e ext of the surroundings, e ext was chosen to be close to e s of the model of liquid (e extce s), so that g VK is very similar to g K. However, an accurate calculation of e s requires that g K be defined from g VK using Eq. (3). Besides, how much the molecules located within a distance R from a reference molecule contribute to the Kirkwood factor can be estimated by the R-dependent Kirkwood factor G K(R). Typical profiles of G K(R) corresponding to blocks each of 100 000 steps were selected from the long simulation runs. The curves reported in Fig. 1 were obtained from the run (2a) for TIP3P water and run (4a) for OPLS acetonitrile. For the sake of comparison, we give the G K(R) functions of the whole simulation runs and of the MOZ theory. Obviously, the parts of G K(R) at smaller distances obtained by the runs of 100 000 steps do not strongly differ. Thus, for the TIP3P model, the largest deviation of the second maximum of G K(R) at 6 2 is 0.21 (see left arrow in Fig. 1(a)). This should be compared with the large fluctuations of G K(R) at larger distances. Thus, g VK , which is read at the right side of Fig. 1(a) (see right arrow), varies between 2.64 and 4.21. These values of g VK correspond to the dielectric constants between 68 and 135. The strong fluctuation of the part of G K at larger distances is due to a combined effect. The weak orientational correlation between two molecules at larger distances leads to a small YY value of the average hl i lj i. In contrast to the small average Y Y Y Y correlation, the fluctuations of l i d lj about hli d lj i are not strongly decreased at larger distances. However, these fluctuations have to be averaged out in order to obtain Y Y Y Y hl i d lj i. According to Eq. (7), the fluctuations of hli d lj i are amplified by the great number of particles j around a molecule i at larger distances. We conclude that the difficulty of extracting a small averaged property amidst large fluctuations explains why long simulation runs are necessary to obtain an accurate dielectric constant. In Fig. 1, let us focus on the long-range part of G K, which is so difficult to compute. As the dipole–dipole correlation becomes smaller with increasing distance of the molecules, the function G K should show decreasing

Fig. 1. Several R-dependent Kirkwood factors G K(R) obtained by averaging 100 000 time steps of a simulation for the TIP3P model of water (run 2a, 512 molecules) and for the OPLS model of acetonitrile (run 4a, 512 molecules). In comparison, the G K(R) functions of the whole simulation run and of the MOZ theory are given.

oscillations at larger R. This behavior is observed for the functions G K after the whole run. The broad peak after 15 2 reported in the simulation of water is an artifact due to the limited size of the box. It is worth noting that for the short runs the behavior of the G K at larger R can be notably different from the expected one. For instance, the lowest curve of the runs of 100 000 steps shows a monotonic decrease between R=8 and 15 2. 3.2. Assessing the long-range of G K(R) by the MOZ theory A comparison of the G K functions obtained by the MOZ theory and the whole simulation run has shown that the deviations between both methods for models of water and methanol are mainly due to the drastically underestimated first peak of the MOZ function G K(R) [43]. This can be seen

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

9

in Fig. 1(a) for TIP3P water. The relative heights of the second and third maxima with respect to the first minimum after the first peak are in good agreement with the simulations. This shows that for TIP3P water, the MOZ theory correctly predicts the dipole–dipole correlation between molecules separated by a larger distance. The absence of the broad peak after 15 2 observed in the simulation of the water models confirms that it is an artifact of the simulation method. For TIP3P water, we conclude that after the first peak, the oscillations of the G K(R) function, which have notable amplitudes, are predicted with a rather good accuracy by the MOZ theory. The same property is shown in Fig. 1(b) for OPLS acetonitrile and holds for the other liquid models. It should be noted that Stell et al. [3] employed the function G K(R) of a linearized MOZ–HNC approximation to estimate the range of the dipole–dipole correlations in fluids of dipolar hard spheres. They also stressed that knowing the range of these correlations making a significant contribution to e s is necessary for a proper usage of the simulated dielectric results. 3.3. Defining the long-distance region According to the previous property of the oscillating profile of the MOZ function G K(R), we propose to use this function to define a long-distance region beyond a truncation distance R t, after which the changes of G K(R) have a minor influence on the value of the dielectric constant. We employ the MOZ function instead of the simulation results for two reasons. On the one hand, G K(R) is affected by artifacts due to the limited size of the box. On the other hand, as discussed above, its behavior is very noisy and can be very misleading after short simulation runs. To study the influence of the truncation of G K(R) on the dielectric constant of the solvent e s, an R-dependent dielectric constant E s ð RÞ is defined from Eq. (3) using the values of G K(R) instead of g K: ðE s ð RÞ
1Þð2E s ð RÞ þ 1Þ ¼ yD GK ð RÞ: 9E s ð RÞ

ð8Þ

For the four studied models, the results of the MOZ theory are shown in Fig. 2. In comparison, the E s ð RÞ obtained by the long simulation runs (a) are given. In general,pthe ffiffiffi correction using Eq. (3) is only valid for G K (R) at R ¼ 3lbox =2. As the dielectric constants of the external surrounding and studied model are very close and as we are only interested in the overall evolution of E s ð RÞ, we renounce a correction of G K(R). Using the MOZ functions E s ð RÞ, we define truncation distances R t which satisfy the two truncation conditions (TC): !

(TC1) The R-dependent dielectric constant E s ð RÞ calculated for R values larger than R t should agree with e s to

Fig. 2. The R-dependent dielectric constant E s ðRÞ defined in Eq. (8) obtained by the MOZ theory and by the simulations 1a to 4a (a) for the TIP3P and SPC/E models of water, (b) for the model of methanol and (c) for acetonitrile. The arrows mark the truncation distances defined in Section 3.

10

!

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

within 10%. Then, we can conclude that the part of G K(R) after R t is of minor importance for the calculation of the dielectric constant with respect to the statistical precision of the simulation results. An accuracy of 10% was chosen, because it is similar to the relative statistical error observed for the long simulation runs (see Table 2). (TC2) The cut-off R t should be the smallest distance, which corresponds to a local maximum of E s ð RÞ and fulfills condition (TC1). The local maximum must also exist at a very near value Rsim in the function E s ð RÞ of t the simulation. The Rsim value gives the truncation t distance of the simulation. Of course, we could use any other characteristic point of E s ð RÞ (e.g., local minima or points of inflection) to define R t. The only requirement is that the characteristic point does exist in the functions E s ð RÞ of both the simulation and MOZ theory, because the results of the two methods should be combined at this point as shown later.

runs i and i+1 do not show any appreciable correlation provided that l srz5,000. In this case indeed, the structures of the last molecular configurations of run i are only correlated to the structures of a small fraction of the following l sr configurations of run i+1. This observation is coherent with the value of the Debye relaxation time s D in these liquids [45,63], which is similar to the time 7.5 to 10 ps of a 5000 steps run, during which a notable loss of dipole–dipole orientational correlations is then expected. Thus, the statistical tools apply to the (nearly) independent random t variables g K,i V and G K,i . A proper statistical analysis can be carried out for sufficiently large samples with n srz20–40, which sets a lower bound of 100 000 steps to the simulation runs for computingPdielectric constants. Then, P the mean P P values g K V¼ 1=nsr ig VK;i and GtK ¼ 1=nsr i GtK;i should t have Gaussian distributions. How g K,i V and G K,i spread can be measured by the usual unbiased estimators of their variances s g2KVand s G2 KV defined as

For the three liquids studied here, we observed that using local maxima leads to well defined R t values. For both models of water (Fig. 2(a)), the second maximum of E s ð RÞ fulfills conditions (TC1) and (TC2) sim (SPC/E: Rsim t =5.84 2, TIP3P: R t =6.05 2). In the case of methanol (see Fig. 2(b)), already after the first maximum, the R-dependent dielectric constant oscillates around the dielectric constant with a maximal deviation of 3 %. We chose the second maximum at R t=8.84 2 (Rsim t =8.72 2) to define the truncation distance, because we observed that the first maximum of E s ð RÞ of the MOZ theory at 6 2 can become a shoulder in the simulation. For acetonitrile, the behavior of E s ð RÞ is very different from those of the other solvents studied here. Even at large distances, E s ð RÞ can strongly deviate from the dielectric constant e s. The R value of the third maximum satisfies the conditions for the truncation distance (R sim t =14.84 2). The value of Rsim for acetonitrile is about twice as large t as those found for the other liquids made of molecules of similar sizes.

s2gVK ¼

3.4. Statistical errors

and

In the case of our simulations, we want to estimate the statistical accuracies of the factors g VK and G Kt=G K (R tsim), the G K (R) value at the truncation distance. To carry out the statistical analysis, the complete runs were split into n sr short runs of l sr steps. For each short run i, the values t of the random variables g VK,i and G K,i were calculated. These random variables are the sum of random quantities obtained at the l sr steps which are correlated to the extent that they stem from successive configurations, the structures of which are correlated during some characteristic time. The central limit theorem does not hold and t there is no reason why g K,i V and G K,i should have Gaussian distributions. However, we checked that these random variables associated with two successive short

nsr   X 1 P 2 gVK;i
gKV nsr
1 i¼1

ð9Þ

and s2G t ¼ K

nsr  X P 2 1 GtK;i
GtK : nsr
1 i¼1

ð10Þ

We are now in a position to discuss the accuracy of g VK, and G Kt . This will be done for the simulation run (1a) of the SPC/E model, which clearly shows the major problems that may arise. Choose l sr=5000. For run (1a), t n sr=400. The values of g K,i V and G K,i range in the intervals [0.4,8.8] and [1.7,3.4], respectively. We have the mean P P values gKV ¼ 2:37 and GtK¼ 2:56, and the variances s g KV2= 1.7 and sG2 tK = 0.09.The related standard errors are s g KV=1.3 and s GtK=0.3. Now, for a simulation of a given length (n sr l sr) steps, P estimators of the errors of P gKV and GtK are given by sgVK DgVK ¼ f pffiffiffiffiffi ð11Þ ffi nsr sGtK DGtK ¼ f pffiffiffiffiffi ffi; nsr

ð12Þ

respectively. For run (1a), where n sr=400, Dg KV=f 0.065 and DG Kt = f 0.015. The factor f determines the probability P P to find the true g VK in the interval g VK
DgVK ; g VK þ DgVK . P In the case of a Gaussian distribution of g VK, i.e. for reasonably large n sr, this probability is 68%, 95%, and 99.7% for f=1, 2, and 3, respectively. Similar results hold for G Kt. To avoid a systematic overoptimistic estimation of the accuracy of the computed averages, we chose f=2 instead of the usual value [50] f=1. Then, the true values g VK and G Kt are in the intervals 2.37F0.13 and 2.56F0.03, respectively, with a probability of 95%.

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

11

The variation de s of the dielectric constant in terms of the variation dg VK of the intermediate Kirkwood factor is simply obtained by differentiating Eq. (4). We obtain: des ¼

ðes
1Þ2 dgVK 3yD gV2K

ð13Þ

When e extce s and e sH1, we have e sc(9/2)y Dg K, g Kcg VK , so that de s simplifies to des ¼

27 yD dg KV : 4

ð14Þ

This expression gives an approximation of the error of the dielectric constant in terms of the error of dg VK . The error of the dielectric

P constant was also estimated

P by the difference Des ¼ es g VK þ DgVK Þ
es ; where es g VK þ DgVK Þ denotes P the dielectric constant calculated from the factor g VKþDgVK. Both methods give very similar results. For instance, taking l sr=5000 and applying Eq. (14) to the run (1a) of SPC/E water, for which we have y D=6.2424 and Dg VK = 0.13, we obtain De s=5.5, i.e. e s=66F5.5. This shows that the P statistical accuracy of g VK is sufficient, but it required an intensive computation of 2 000 000 steps. Would it have been possible to notably reduce the computation length in a safe manner? The following discussion shows that such a reduction could lead to misleading results. As previously, consider the long simulation run (1a) of the SPC/E model, which is now split into 20 shorter runs of 100 000 steps. Again, choose l sr=5000, so that n sr=20. For each run of 100 000 steps, P the values g VK and Dg VK were calculated and the results are shown in Fig. 3. Obviously, Dg VK strongly fluctuates and P does not always indicate the reliability of the computed g VK values. For example, the 8th run of 100 000 steps provides a Kirkwood factor of 1.55F0.27, which corresponds to the dielectric constant 37F8 according to Eqs. (1) and (3). The comparison with the result 66F6 obtained after the whole simulation shows that the error gives a completely wrong impression of the accuracy of the computed dielectric constant. The statistical interpretation of the misleading prediction of the 8th run is immediate. According to the long run (1a), the btrueQ variance of g VK,i is s g KV2=1.7. From Eq. (11), a reasonable estimate of the error of any short run with n sr=20 is Dg VK = 0.6, as given by the yaxis value of the dashed line of Fig. 3(b). The error corresponds to a factor f=2. For n sr=20, the probability of P finding g VK outside the confidence interval [2.37F0.6] is as high as 5%. Among 20 other values, the presence of the P value g VK ¼ 1:55 of the 8th run, which is just 0.2 unit smaller than the lower bound 1.77 of the confidence interval is quite possible. Furthermore, the values of Dg VK are only estimators of the errors, and for small samples of size n sr=20, the probability that the btrueQ Dg VK is twice as large as the estimate calculated from Eqs. (9) and (11) is of the order of a few percents. This explains the misleading Dg VK =0.27 of the 8th run, which corresponds to an

Fig. 3. (a) The Kirkwood factors and (b) corresponding errors calculated from Eq. (11) for different simulation runs of 100 000 steps, using small runs of l sr=5000 steps and f=2. The runs of 100 000 steps were derived from the run (1a) of the SPC/E model of water. In Fig. 3(a), the Kirkwood factor obtained after the whole simulation run is marked by a dashed line. In (b), we show the error for a run of 100 000 steps obtained by Eq. (11) using the standard error of the whole run (marked by whole run).

accumulation of g VK,i values in a narrow region markedly below g VK =2.37. To sum up, for random variables such as g K,i V , which spread over a large interval [0.4,8.8] and have a mean value P g VK ¼ 2:37 and a standard error s g KV=1.3 of nearly the same sizes, large samples should be used in order to avoid not only a poor estimate of g VK , but also the wrong and dangerous impression that this estimate is reliable since its apparent error Dg VK is small. As a practical rule, Dg VK P should be made smaller than g VK=10 in order to guarantee a satisfactory relative precision. Here, runs longer than 1 000 000 steps are required. On the other hand, for random t variables such as G K,i , the variance is 10 times smaller than

12

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

that of g VK,i and runs of only 100 000 steps can be used with good confidence to get reasonable estimates of G Kt of satisfactory relative precision. In Table 3, we give the statistical errors sg KV and s G Kt for the small runs which correspond to l sr=20 000 and are obtained by splitting the long simulation runs (1a), (2a), (3a), and (4a) of the SPC/E water, TIP3P water, H1 methanol, and OPLS acetonitrile, respectively. A splitting into shorter (l sr=5000) and longer blocks (l sr=100 000) gave similar errors. The estimates of the errors Dg VK and DG Kt for runs of 100 000 and 1 000 000 steps calculated from these variances according to Eq. (11) and (12) are also shown in Table 3. In the case of the runs (1a), (2a), and (3a) for the models of water and methanol, we have Dg VK z3.5DG Kt, so that for each of these models simulation runs at least 10 times longer would be necessary to reach a mean intermediate Kirkwood factor g VK with a statistical error similar to that of G Kt. For acetonitrile, Table 3 shows that Dg VK is still significantly larger than the corresponding DG Kt, even though their difference is smaller. This quantitative statistical analysis confirms the qualitative picture derived from Fig. 1. For each of the studied models, the noisy contribution of g VK –G Kt is the main reason why an accurate dielectric constant cannot be computed from short simulation runs. Because of the property (TC1) of the truncation distance R t, the residual g VK –G Kt should have a minor effect on the dielectric constant. In addition, this long-distance term is the contribution to g VK which is the most affected by the effects of the limited size of the box as shown in Fig. 4. In this figure, the R-dependent Kirkwood factors G K(R) of the SPC/E model are displayed for the very long runs (1a) and (1b) with 512 and 256 molecules, respectively. The two functions are in good agreement up to the third maximum. Then, they start to deviate and yield intermediate Kirkwood factors (256 molecules: 2.58F0.17, 512 molecules: 2.37F0.15), which are different with a high probability. We conclude that the contribution g VK –G Kt is not well known due to statistical and systematic errors. In what hybrid methods which employ P follows, we propose P G Kt instead GKPof g VK to obtain dielectric constants. These methods are justified as long as the systematic error Table 3 Square roots s g KV and s G tK of the variances and errors of g VK and G Kt for the simulation runs (1a), (2a), (3a) and (4a) of models of liquid water (SPC/E and TIP3P), methanol (H1), and acetonitrile (OPLS) Liquid model

s g KV

s G tK

100 000 steps DG Kt Dg VK

1 000 000 steps Dg VK DG Kt

Water (SPC/E) Water (TIP3P) Methanol (H1) Acetonitrile (OPLS)

0.75 1.10 0.71 0.14

0.18 0.16 0.20 0.10

0.67 0.98 0.63 0.13

0.21 0.31 0.20 0.04

0.16 0.14 0.17 0.09

0.05 0.04 0.05 0.03

The errors were calculated from the variances according to Eqs. (11) and (12), using small runs of l sr=20 000 steps and f=2. The variations between the variances of the liquids are partly due to the different lengths of the simulation runs.

Fig. 4. The R-dependent Kirkwood factor G K(R) of the SPC/E model for the run (1a) with 512 molecules and for the run (1b) with 256 molecules. P

P

introduced by approximating g VK by GtK (about 0.2 for run P (1a)) is smaller than the statistical error on g VK (when 212 ffi for run (1a), N MD l sr=20 000, f=2, we have DgVK ¼ pffiffiffiffiffiffi NMD being the number of simulation steps).

4. The MD/MOZ hybrid methods 4.1. Approximate expressions of the dielectric constant We propose two methods which avoid a direct computation of the part g KV –G Kt by simulations. Both methods employ information derived from the MOZ theory. For this reason, they belong to the hybrid methods combining the results of simulations and of Ornstein–Zernike theories. The first method, which we call the btruncation methodQ’, rests on the fact that the part of gVK
G Kt is usually small (see Fig. 1), in particular, for water and methanol. Therefore, this part is neglected and the dielectric constant is computed from G Kt derived from simulations: P ðets
1Þð2ets þ 1Þ ¼ yD GtK t 9es

truncation method

ð15Þ

In the case of the truncation method, the function G K (R) of the MOZ theory is used to find a reasonable  MOZregion for the truncation distance. The ratio E s RMOZ =es observed t in the MOZ theory indicates the likely systematic error made by the truncation. For example, for the TIP3P model, we predict from Fig. 2 that the dielectric constant e ts should be too low with an error V5%, whereas for the methanol and the SPC/E model the truncation method should yield excellent estimates of e s. G Kt is not corrected according to Eq. (3), because the first part of G K(R) up to the second peak is not strongly affected by a different dielectric

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

constant of the reaction field. In the case of lattice sum techniques applied to SPC/E water, this can be seen in Fig. 5 where the functions G K(R) obtained for 256 molecules by run (1a) with e ext=70 and (1c) with e ext=l are compared. The alternative technique, which we call the btailcorrectionQ method, employs the residual g KMOZ
G Kt,MOZ obtained by the MOZ theory for the same liquid as an approximation of this part of g K. This yields the formula: ðetcs
1Þð2etcs þ 1Þ 9etcs nP o ¼ yD GtK þ gKMOZ
Gt;MOZ K

tail
correction method ð16Þ

The dielectric constants computed by both methods are compared to the simulation results and MOZ values in Table 2. The errors of the dielectric constant obtained by simulations were estimated by the difference P Des ¼ es ðgVK þ DgVK Þ
es as explained in Section 3. The errors of both hybrid methods were derived from the errors of DG Kt. For instance, in the case of the tail-correction t t method the dielectric constants e tc s (G K+DG K) were comPt P t t puted using GK þ DGK instead of GK in Eq. (16). The P difference between etcs GtK þ DGtK and etcs yields the statistical error De tc s . Both hybrid methods yield dielectric constants in good agreement with the simulation results. Only the dielectric constant for the TIP3P model obtained by the truncation method is underestimated in comparison with the simulation result. This agrees with our predictions. 4.2. Influence of the simulation parameters The comparison of the dielectric constants derived from the hybrid methods for the runs with 512 and 256 molecules

Fig. 5. The R-dependent Kirkwood factor G K(R) of the SPC/E model for the run (1b) and (1c) using Ladd sum with e ext=70 and e ext=l and for the run (1d) using the reaction field method (e ext=l).

13

(runs (a) and (b)) indicates that the dependence of the hybrid results on the size of the box is very small. Usually, the dielectric constant e ext of the external dielectric surrounding employed during the simulation was taken to be close to the dielectric constant of the studied model. Fig. 5 shows that the choice of e ext has an influence on the G K functions and, therefore, also on the result of the MD/MOZ hybrid method. When the hybrid method is used to study a liquid model of unknown dielectric constant, the MOZ theory usually gives a valuable estimate for e ext. Nevertheless, the results should not depend too much on the chosen e ext. In order to test this requirement, we calculated the dielectric constant by the hybrid methods for the run (1c) with e ext=l and compared the obtained values to their counterparts of run (1a) with e ext=70 (see Table 2). For instance, in the case of the tail-correction method, the dielectric constant increases from 71.6 to 73.0. This indicates that the dielectric constant of the surrounding e ext has little influence on the results of the hybrid methods. However, under the conducting boundary conditions expressed by e ext=l, there is a long-range increase of G K(R) beyond the second maximum as shown previously [45] or in Fig. 5. This is consistent pffiffiffi with the

 fact that according to Eq. (3), g VK ¼ GK R ¼ 3lbox =2 ¼ 32 gK is notably larger than g K. The long-range increase of G K(R) precludes the use of a too large truncation distance. This is particularly true in the case of the reaction field method. 4.3. Statistical convergence of the hybrid methods According to Eq. (11), the statistical error of the dielectric constant obtained after a n sr V l sr run of steps can be estimated from its counterpart De of a run of n sr l sr pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi steps by the formula nsr =nsrV De. From Table 2, after 100 000 steps, the statistical errors of the hybrid methods for the runs (a) of the SPC/E and TIP3P models of water, the methanol model, and the acetonitrile model, are 4.5, 4.1, 2.2, and 1.7, respectively. These values are smaller than the estimates 8.8, 14, 3.9, 1.2 found by the usual simulation method after 1 000 000 steps. To illustrate the efficiency of the new hybrid methods, consider the dielectric constants of TIP3P water and methanol for the runs of 100 000 steps which were generated by splitting the long runs (2a) and (3a). The results obtained by simulation and the tailcorrection method are compared in Fig. 6. For instance, in the case of water, the largest and smallest values of the dielectric constant obtained by the hybrid method are 93 and 99, while the simulation results vary between 68 and 138! As discussed in Section 3, the statistical analysis of short simulation runs can give a false impression of the accuracy of the dielectric constant (see Figs. 1(a) and 2(a)). In Fig. 7, for SPC/E water, we compare the errors of the dielectric constants of the simulation and tail-correction method for the runs of 100 000 steps deduced from the long run (1a). The errors Dg KV and DG Kt were computed according to Eqs. (11) and (12), respectively, using the standard errors

14

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

SPC/E water by comparing the results of runs (1d) and (1c) in Table 2. Recall that the truncation distance is at the second maximum of E s ð RÞ displayed in Fig. 2(a). We observe that the results of the reaction field methods c59 are notably different from those of the lattice sum techniques c72, because of the completely changed behavior of G K(R) already at intermediate R, as shown in Fig. 5. Now, in the context of the reaction field method, we study whether the use of a smaller truncation distance can improve the stability of the hybrid methods. For this purpose, consider the SPC/E model of water. As observed in Section 3, the function G K(R) calculated by the MOZ theory has a correct oscillatory behavior already after the first maximum. Therefore, in the case of the tail-correction method, the first maximum of E s ð RÞ displayed in Fig. 2(a) can be chosen as the truncation distance R t to estimate the dielectric constant. For the simulation run (1a) using the lattice sum techniques, this gives a dielectric constant of 67.7F0.2 instead of 71.6F1.0, which indicates that there is a systematic difference of 5% due to the shorter truncation distance. Attempts with the other models of water and methanol yield similar deviations. Moreover, the shorter is the truncation distance, the less affected are the hybrid results for the dielectric constants when the reaction field method replaces the lattice sum techniques. For instance, when R t is taken at the first maximum, the run (1d) with the reaction field method gives a dielectric constant of 65.1F0.4 (result 1db in Table 2) instead of 67.7F0.2 for the run (1a) with the lattice sum techniques. The small difference between these values should be compared to the differences between their counterparts 58.1F2.6 (run (1d) in Fig. 6. Dielectric constants of the TIP3P model of water and the H1 model of methanol for different simulation runs of 100 000 steps. The results by simulations and by the tail-correction hybrid methods are shown. The runs of 100 000 steps were obtained by splitting the runs (2a) and (3a).

computed for the runs of 100 000 steps with l sr=5000. These values are compared to the errors obtained using the standard errors of the long simulation runs. Obviously, the errors of e tc s , calculated from the runs of 100 000 steps and P deduced from those of GtK are much more reliable than the P errors of e s derived from those of gVK . Thus, another advantage of the hybrid methods is that the accuracy of their results can be derived from short simulation runs. 4.4. Is the faster reaction field method applicable? Finally, let us examine whether the function G K(R) obtained by a simulation using the reaction field method is appropriate for the hybrid method. If it were true, the hybrid methods would become even more efficient, because the reaction field approach is about three times faster than the Ladd sum techniques. We test this possibility in the case of

Fig. 7. Errors of the dielectric constants of the SPC/E model of water. The values were calculated by Eqs. (11) and (12) for different simulation runs of 100 000 steps. The results by simulations and by the tail-correction hybrid methods are shown. The runs of 100 000 steps were obtained by splitting the run (1a). We also show the error for a run of 100 000 steps obtained by Eqs. (11) and (12) using the standard error of the whole run (marked by whole run).

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

Table 2) and 71.6F1.0 (run (1a) in Table 2) when R t is taken at the second maximum. During a simulation, the E s ð RÞ profile up to the first maximum converges even more rapidly than its extension up to the second maximum, which corresponds to more distant molecules. Therefore, the statistical error of the dielectric constant after 100 000 steps has a value of 0.7 using the first maximum compared to 2.6 with the second maximum. This shows that combining the reaction field method with a small truncation distance can give an estimate of the dielectric constant with an accuracy of 10% even after very small runs of 40 000 steps. The dielectric constants of the other liquid models obtained by the same procedure support this conclusion. For the models of liquids studied here, a run of 40 000 steps using the reaction field method takes a CPU time of less than 1 h on a DEC a 4100. On the same computer, 36 days are required by the run (1a) in order to obtain a dielectric constant with a statistical error of about 10%. To be complete, note that using an atomic cut-off instead of the molecular cut-off of this work may result in a favourable reduction of the systematic error caused by the reaction field method. Indeed, recent simulations have shown that the G K(R) function of the reaction field method and its counterpart of the lattice sum techniques are nearer when an atomic cut-off is used [62].

5. Conclusion The R-dependent Kirkwood factor G K(R) obtained by simulation and the MOZ theory have been investigated in detail for several liquid models. For each model, using both the functions G K(R) of the MOZ theory and simulation, we derived a truncation distance Rsim t , beyond which G K(R) is of minor importance for the dielectric constant. A detailed statistical analysis of the intermediate Kirkwood factor g KV and G Kt=G K(Rsim t ) using the results of long simulation runs has shown that during a simulation, the residual g KV–G Kt has a slow statistical convergence responsible for that of the dielectric constant. The residual g KV–G Kt is also affected by systematic errors such as those caused by the limited size of the simulation box. To reduce the CPU time needed to obtain well-converged dielectric constants, we propose the following hybrid method, which avoids the direct computation of the longrange part of G K(R). First, G K(R) is computed by a short simulation run (50 000 to 100 000 steps) and the MOZ theory. Second, a characteristic truncation distance (local extremum, inflection point) is derived from the R-dependent dielectric constant E s ð RÞ obtained from the G K(R) function of the MOZ theory. The simulated G K(R) profile should have a truncation point at a similar distance Rsim and the t variation of the G K(R) function of the MOZ theory after Rsim should have a minor effect on the dielectric constant. t Furthermore, when e ext=l, the simulated G K(R) at large R

15

markedly exceeds the value, it would have in an infinite liquid, and this precludes the use of a too large truncation distance. Third, the dielectric constant is calculated either by Eq. (15) (truncation method) or Eq. (16) (tail-correction and E tcs gives an procedure). Comparing E ts ¼ as Rsim t idea of the precision of the truncation method. The statistical error of the dielectric constant can be derived from that of G Kt given by Eq. (12). Following Allen and Tildesley [50], bComputer simulation is an experimental science in so far as the results may be subject to systematic and statistical errorsQ. The MOZ formalism rests on a closure approximation. For a new class of fluids, it provides properties, the accuracy of which is questionable until they have been tested against their counterparts obtained from well-converged simulation data and/or experiment on a similar actual system. The situation is analogous to that encountered when applying the popular density functional theory [64] (DFT) of quantum chemistry, for which density functionals such as the exchangecorrelation energy are approximate. Thus, for a new class of fluids, a conservative use of the hybrid methods should begin by comparing the MOZ results to those of a reference long simulation run with the lattice sum techniques. In particular, within the statistical errors, the MOZ theory and simulation should give similar G K(R) functions at large distances. Then, the hybrid methods should be applied by using the indicated recipes and its predictions checked again the simulated values. If the accuracy of a hybrid method is satisfactory, i.e. its systematic error is smaller than that expected from the statistical error on the intermediate Kirkwood factor g KV of the reference simulation run, this hybrid method can be reasonably extended to similar fluids in rather near thermodynamical states, using short simulation runs. For all the liquid models studied here and previously the HF3 model of hydrogen fluoride [65], we have shown that the dielectric constants obtained by the MD/MOZ hybrid methods are in good agreement with the simulated values. The advantage of the truncation method is that we can obtain the dielectric constants without solving the MOZ theory for the new interaction model when the region of Rsim and the t possible errors were established for a similar model. With the help of the hybrid methods, the dielectric constant can be calculated with a precision of at least 10% after 100 000 steps. Moreover, the dependence on the size of the simulation cell is small. Lattice sum techniques give more accurate estimates of e s than the reaction field method. Nevertheless, the latter can be employed with a very short truncation distance and still gives a quite reasonable approximation of e s provided that the tail-correction method is employed. For a very short truncation distance sim Rsim t ,G K(R) (R V R t ) has a very rapid convergence, so that the statistical error nearly vanishes for a run of just 100 000 steps. Unfortunately, combining the reaction field method and a short truncation distance causes a systematic error of about 10%.

16

J. Richardi et al. / Journal of Molecular Liquids 117 (2005) 3–16

To sum up, the MD/MOZ hybrid techniques enable us to calculate the dielectric constant of a liquid model with the same accuracy as a very long simulation, but at least 10 times more rapidly. In other words, rather accurate dielectric constants of molecular liquid models can be estimated by a relatively short standard MD (or Monte Carlo) simulation routinely performed to obtain thermodynamic and structural properties. Furthermore, it should be emphasized that the long simulation runs of this work and the hybrid methods give dielectric constants with similar precisions of c10%. A decrease of the statistical error of the simulated values to 1% would result in unmanageable computations. The efficiency of the hybrid methods will also be particularly valuable for more complicated models, in particular those including the polarizability explicitly and requiring several weeks of CPU on modern workstations to simulate the dielectric properties [56,66]. Moreover, for liquids showing very long-range orientational order such as formamide and N-methylformamide (NMF) [38], the MOZ theory and hybrid methods might even be the only approaches with reasonable computational time demands.

Acknowledgment Johannes Richardi gratefully acknowledges the attribution of a Marie Curie postdoctoral fellowship by the European Commission.

References [1] A. Frqhling, Cours d’Electricite´, Vol. I, Dunod, Paris, 1966, pp. 92 – 115. [2] C.J.F. Bfttcher, 2nd editionTheory of Electric Polarization, Vol. I, Elsevier, Amsterdam, 1973. [3] G. Stell, G.N. Patey, J.S. Hoye, Adv. Chem. Phys. 48 (1981) 183. [4] P. Madden, D. Kivelson, Adv. Chem. Phys. 56 (1984) 467. [5] E.L. Pollock, B.J. Alder, Physica 102A (1980) 1. [6] M. Neumann, O. Steinhauser, G.S. Pawley, Mol. Phys. 52 (1984) 97. [7] D.M.F. Edwards, P.A. Madden, I.R. McDonald, Mol. Phys. 51 (1984) 1141. [8] Z. Wang, C. Holm, H.W. Mu¨ller, Phys. Rev., E 66 (2002) 021405. [9] L.X. Dang, B.M. Pettitt, J. Phys. Chem. 94 (1990) 4303. [10] P.G. Kusalik, Mol. Phys. 73 (1991) 1349. [11] D.J. Adams, E.M. Adams, Mol. Phys. 42 (1981) 907. [12] H.E. Alper, R.M. Levy, J. Chem. Phys. 91 (1989) 1242. [13] X. Ni, R.M. Fine, J. Phys. Chem. 196 (1992) 2718. [14] P.G. Kusalik, Mol. Phys. 76 (1992) 337. [15] B.M. Ladanyi, M.S. Skaf, Annu. Rev. Phys. Chem. 44 (1993) 335. [16] T. Fonseca, B.M. Ladanyi, J. Chem. Phys. 93 (1990) 8148. [17] P.G. Kusalik, M.E. Mandy, I.M. Svishchev, J. Chem. Phys. 100 (1994) 7654. [18] Z. Kurtovic´, M. Marchi, D. Chandler, Mol. Phys. 78 (1993) 1155. [19] J.-P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. [20] L. Blum, A.J. Torruella, J. Chem. Phys. 56 (1972) 303. [21] L. Blum, J. Chem. Phys. 57 (1972) 1862. [22] S.L. Carnie, G.N. Patey, Mol. Phys. 47 (1982) 1129.

[23] P.H. Fries, G.N. Patey, J. Chem. Phys. 82 (1985) 429. [24] F. Lado, M. Lombardero, E. Enciso, S. Lago, J.L.F. Abascal, J. Chem. Phys. 85 (1985) 2916. [25] E. Lomba, C. Martı´n, M. Lombardero, Mol. Phys. 77 (1992) 1005. [26] C. Martı´n, M. Lombardero, E. Lomba, J. Chem. Phys. 98 (1993) 6465. [27] C. Martı´n, M. Lombardero, M. Alvarez, E. Lomba, J. Chem. Phys. 102 (1995) 2092. [28] C. Martı´n, E. Lomba, M. Lombardero, F. Lado, J.S. Hbye, J. Chem. Phys. 100 (1994) 1599. [29] F. Lado, E. Lomba, M. Lombardero, J. Chem. Phys. 103 (1995) 481. [30] M. Alvarez, F. Lado, E. Lomba, M. Lombardero, C. Martı´n, J. Chem. Phys. 107 (1997) 4642. [31] J.A. Anta, E. Lomba, M. Alvarez, M. Lombardero, C. Martı´n, J. Phys. Chem., B 101 (1997) 1451. [32] J. Richardi, P.H. Fries, R. Fischer, S. Rast, H. Krienke, J. Mol. Liq. 73,74 (1997) 465. [33] J. Richardi, P.H. Fries, R. Fischer, S. Rast, H. Krienke, Mol. Phys. 93 (1998) 925. [34] J. Richardi, P.H. Fries, H. Krienke, J. Chem. Phys. 108 (1998) 4079. [35] R. Fischer, J. Richardi, P.H. Fries, H. Krienke, J. Chem. Phys. 117 (2002) 8467. [36] P.H. Fries, W. Kunz, P. Calmettes, P. Turq, J. Chem. Phys. 101 (1994) 554. [37] P.H. Fries, J. Richardi, H. Krienke, Mol. Phys. 90 (1997) 841. [38] J. Richardi, H. Krienke, P.H. Fries, Chem. Phys. Lett. 273 (1997) 115. [39] J. Richardi, P.H. Fries, H. Krienke, J. Phys. Chem., B 102 (1998) 5196. [40] J. Richardi, P.H. Fries, H. Krienke, Mol. Phys. 96 (1999) 1411. [41] I. Paci, N.M. Cann, J. Chem. Phys. 115 (2001) 8489. [42] P.H. Fries, G. Ferrante, E. Belorizky, S. Rast, J. Chem. Phys. 119 (2003) 8636. [43] J. Richardi, C. Millot, P.H. Fries, J. Chem. Phys. 110 (1999) 1138. [44] M. Lombardero, C. Martı´n, S. Jorge, F. Lado, E. Lomba, J. Chem. Phys. 110 (1999) 1148. [45] P. Hfchtl, S. Boresch, W. Bitomsky, O. Steinhauser, J. Chem. Phys. 109 (1998) 4927. [46] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [47] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [48] M. Haughney, M. Ferrario, I.R. McDonald, J. Phys. Chem. 91 (1987) 4934. [49] W.L. Jorgensen, J.M. Briggs, Mol. Phys. 63 (1988) 547. [50] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [51] D.J. Evans, Mol. Phys. 34 (1977) 317. [52] D.J. Evans, S. Murad, Mol. Phys. 34 (1977) 327. [53] A.J.C. Ladd, Mol. Phys. 33 (1977) 1039. [54] A.J.C. Ladd, Mol. Phys. 36 (1978) 463. [55] M. Neumann, Mol. Phys. 60 (1987) 225. [56] C. Millot, J.-C. Soetens, M.T.C. Martins Costa, Mol. Simul. 18 (1997) 367. [57] C.M. Carlevaro, L. Blum, F. Vericat, J. Chem. Phys. 119 (2003) 5198. [58] M. Rami Reddy, M. Berkowitz, Chem. Phys. Lett. 155 (1989) 173. [59] I.M. Svishchev, P.G. Kusalik, J. Chem. Phys. 98 (1994) 728. [60] J.W. Essex, Mol. Simul. 20 (1998) 159. [61] R.D. Mountain, J. Chem. Phys. 107 (1997) 3921. [62] P.H. Hqnenberger, W.F. van Gunsteren, J. Chem. Phys. 108 (1998) 6117. [63] M.D. Zeidler, Ber. Bunsenges. Phys. Chem. 69 (1965) 659. [64] D. Joubert, J. Chem. Phys. 110 (1999) 1873. [65] P.H. Fries, J. Richardi, J. Chem. Phys. 113 (2000) 9169. [66] J.-C. Soetens, M.T.C. Martins Costa, C. Millot, Mol. Phys. 94 (1998) 577.