J. Phys. Chem. B 2001, 105, 7783-7791
7783
Fluctuating Charge Study of Polarization Effects in Chlorinated Organic Liquids Estanis Llanta* Departament de Fı´sica i Enginyeria Nuclear, UniVersitat Polite` cnica de Catalunya, Campus Nord B4-B5, Barcelona 08034, Spain
Koji Ando† School of Chemistry, UniVersity of Birmingham, Edgbaston, Birmingham B15 2TT, U.K.
Rossend Rey‡ Departament de Fı´sica i Enginyeria Nuclear, UniVersitat Polite` cnica de Catalunya, Campus Nord B4-B5, Barcelona 08034, Spain ReceiVed: January 31, 2001; In Final Form: May 22, 2001
The role of polarization has been systematically studied for the methylchloromethane family ((CH3)4-nCCln) of organic liquids. Special attention has been paid to dynamical properties (translational diffusion and rotational relaxation) for which good agreement with experiment is obtained, although structural and electrostatic effects have been addressed as well. Molecular dynamics simulations have been performed, with and without the inclusion of polarization, over a substantial range of temperatures. Polarizability has been handled with the chemical potential equalization (fluctuating charge) method, with a transferable set of parameters having been developed for the methyl group. As a general rule, the inclusion of polarization does not affect structure and slightly slows down the dynamics at all temperatures. The overall muted influence of polarization contrasts with the substantial induced dipole moments obtained, exceeding what had been found in other organic (albeit hydrogen-bonding) liquids.
I. Introduction Chlorinated organic molecules play an important role in many industrial applications. Because of their widespread use, they have an environmental impact as groundwater pollutants,1 and as unwanted byproducts of chlorination, they are also present in drinking water.2 It is therefore of great interest to gain a molecular level knowledge of their properties in solution. In a first step along the systematic study of the methylchloromethane family ((CH)4-nCCln) of pure liquids (which includes such important examples as carbon tetrachloride and 1,1,1-trichloroethane), we have previously addressed their equilibrium properties.3 The force fields constructed to this end were highly successful as far as thermodynamic and structural properties are concerned, allowing a detailed explanation of subtle X-ray diffraction differences among the family.3 In this work, we mainly address the study of basic dynamic properties such as translational diffusion and rotational relaxation. In principle, this calculation might pose a challenge on the goodness of the nonpolarizable molecular models used so far, as it is commonly accepted that polarization may significantly contribute to liquid dynamics. The rationale is that a constant (effective) dipole moment might be able to reproduce equilibrium properties but not the fluctuations critical for the dynamics.4 We have completed an extensive series of simulations for all five pure liquids (with and without polarization and at different temperatures) in order to clarify the role of * To whom correspondence should be addressed. E-mail: estanislau.
[email protected]. † E-mail:
[email protected]. ‡ E-mail:
[email protected].
nonadditive forces. It is interesting to note that this family constitutes a remarkable case study because they are all molecules of almost identical size, with both polar and nonpolar species, which makes a systematic comparison possible. A considerable number of studies have addressed the issue of polarization effects. Liquid water has been the most intensively studied case5-20 with somewhat confusing results. A large effort has been devoted to the design of polarizable models with overall better performance, in the bulk liquid, than pair-additive ones. We understand that this line of research, while it has resulted in potentially useful force fields, has not provided a clear notion of which are the effects of polarization because of contradictory results. Even more disappointing is the fact that polarizable models have not proven to be clearly superior in the liquid phase for which they had been optimized. A second line of research has focused on transferability (a highly desirable goal if we wish to study complex environments from a sound basis), starting with the landmark paper of Barnes et al.5 This work highlighted how the inclusion of polarization resulted in a largely transferable model and provided a physical picture of polarization effects in different phases. Unfortunately, the complexity of the model, together with the good performance that simple nonadditive models have shown for the pure liquid, has precluded its use. Nevertheless, computationally more tractable schemes have been devised with the emphasis still put on transferability, trying to unambiguously compare with similar nonpolarizable models. To cite one example, the work of Bernardo et al.15 is illustrative: on the basis of the LennardJones parameters of the popular SPC model, charges are adjusted to the gas-phase dipole moment and atomic polarizabilities are added. A careful treatment of polarizability results in a suc-
10.1021/jp010390r CCC: $20.00 © 2001 American Chemical Society Published on Web 07/18/2001
7784 J. Phys. Chem. B, Vol. 105, No. 32, 2001 cessful model in both the gas and liquid phases, not requiring a fine-tuning of non-Coulomb interactions for each phase. The effects of polarizability are physically meaningful: a nonnegligible slowing down of dynamic properties (translation and rotation), together with a remarkable induced dipole moment. Moreover, effects that would not be present in nonpolarizable models arise; we have, for instance, the effect that induced and permanent dipole moments are not parallel.15 The sole inclusion of polarizability does not guarantee that accurate transferable models shall be obtained, as the study of phase coexistence has shown. Other parts of the potentials, like the short-range interactions, might have to be improved as well.21-24 Comparably, much less attention has been paid to organic liquids, with most of the work devoted to the optimization of liquid-phase potentials25-29 or to the study of inhomogeneous systems.30-35 The only systematic studies that we know of, similar to the ones just described for water, have focused on liquid alcohols.36,37 Gao et al.36 developed a polarizable force field for simulation of five different alcohols which was the first systematic study of a class of organic liquids with nonadditive potentials. They paid special attention to the contribution of polarizability, finding it was significant for energy (10%-20% of the total energy), with a dramatic increase in the molecular dipole moment (with induced dipoles of about 0.5 D). These effects were validated36 with quantum mechanical and molecular mechanical (QM/MM) methods. Here, we have intended to follow a similar route for another family of organic liquids. Our models have the experimental gas-phase dipole moments and anisotropic polarizabilities (see below) and charges adjusted to reproduce the electrostatic potential surface around the isolated molecule.3 These models are used to simulate the liquid phase, reproducing a variety of liquid properties quite satisfactorily, which should help with understanding the effects of polarization. It is interesting to note that the liquids studied here have some peculiarities which make it difficult to predict the net effect of polarization. First, electrostatic energy shows a remarkably low contribution in the corresponding fixed charge models, amounting to a maximum of a 8% of the total energy,3 in contrast with almost the opposite situation in liquid water or with an intermediate situation for some alcohols.36 In this sense, the present case constitutes an extreme scenario, as far as the balance of electrostatic vs nonelectrostatic energies is concerned. These considerations certainly suggest that polarizability should barely affect liquid properties. Actually, we have found that the very values of the (fixed) charges in the nonpolar cases (n ) 0, n ) 4) have no effect whatsoever on structural properties.3 However, it is also the case that the polar molecules in this family show substantial dipole moments (about 2 D) and for polar and nonpolar molecules alike, polarizabilities are of roughly 10 Å3, approximately a factor of 7 larger than that of water. From the latter arguments, one is led to think that polarization might be relevant after all. Moreover, the globular shape of these molecular types might help to increase polarization effects, as rotation is not hindered by steric effects but by the strength of interactions. Several methods have been described to include polarization in molecular dynamics simulations.5,6,8,10 One of the most appealing from the computational and physical standpoint is the chemical potential equalization method (CPE),38 in which partial charges within the molecule are allowed to fluctuate in order to adapt to the changing electrostatic environment. This approach has also been formulated in terms of an extended Lagrangian in order to speed up the computation.16 The method
Llanta et al. has an extreme computational simplicity as no dipole-dipole (or higher multipole) interactions have to be included and can be implemented within a conventional MD code without any significant modifications, which has facilitated its use in a wide variety of problems.16,32,34,39-41 Moreover, it allows for an easy fine-tuning of parameters to molecular properties of interest. In this connection, a substantial part of this work has aimed to develop a set of parameters for the methyl group (in the united atom approximation) which are transferable over the whole family and which therefore might also be useful in other molecular environments. Furthermore, the method is specially appropriate for the present family because of the globular shape of the molecules which allows for the generation of any induced dipole, a possibility which is precluded in other cases. We have, for instance, the fact that in the important case of water, because of the planar molecular geometry, induced moments out of the molecular plane are not possible.16 The paper is organized as follows. In the next section, we describe the chemical equalization method with a detailed account of the parametrization for the methyl group. In section III, the details of the simulations are given. Section IV contains the results obtained and their discussion. We summarize the main conclusions in section V. II. Implementation of Polarization Forces A. Chemical Potential Equalization. In the CPE method, partial charges are determined for a given liquid configuration by minimization of the electrostatic energy. In this way, the usual discrete point charge models used in MD simulation with additive interactions are readily adapted to include nonadditive forces. The function minimized for an isolated molecule is the total molecular energy expanded to second order in the partial charges:
Emolec ) E0(r) +
1
1
Jij(rij)qiqj ∑i χ0i qi + 2∑i J0i qi2 + 2∑i ∑ j*i
(1)
where E0(r) denotes the charge-independent contribution, χ0i (“atomic electronegativity”) and J0i (“atomic hardness”) are parameters characteristic of the atomic site i, and Jij(r) is a screening function, which is usually computed as the Coulomb integral of Slater ns atomic orbitals. In practice, the atomic parameters might depend on the molecular environment. Partial charges are found by imposing the condition that atom electronegativities within the molecule (χi ≡ ∂Emol/∂qi) equalize (χ1 ) χ2, ..., χ1 ) χM) while maintaining overall neutrality ( M qi ) 0). The following set of M linear equations results: ∑i)1
JQ ) ∆K
(2)
where Q stands for the atomic charge vector (q1, ..., qM), ∆K contains differences between atomic electronegativities (∆K ) (0, χ01 - χ02, ...,χ01 - χ0M)), and J is the matrix defined as
{
1 for i ) 1 Jij ) J - J ij 1j for i * 1
(3)
(having introduced the notation Jii ≡ J0i ). It is interesting to note that the polarizability tensor is easily obtained from
Rγβ ) rβJ-1∆rγ
(4)
where rβ ) (rβ1, ..., rβM) and ∆rγ ) (0, rγ2 - rγ1, ..., rγM rγ1), rβi being the β Cartesian coordinate of atom i.
Polarization Effects in Chlorinated Organic Liquids
J. Phys. Chem. B, Vol. 105, No. 32, 2001 7785
TABLE 1: Polarizability Tensor Components (in Å3) of (CH3)4-nCCln from MP2/aug-cc-pVDZ Calculations
TABLE 2: Partial Charges Computed ab Initio for the Isolated Molecule
n
4
3
2
1
0
n
qMe
symmetrya Rxx Ryy Rzz R jb R j (exptl)c
Td 9.96 9.96 9.96 9.96 11.2
C3V 10.29 10.29 8.99 9.86 10.7
C2V 10.67 8.96 9.66 9.76 10.9
C3V 9.24 9.24 10.54 9.67 12.5
Td 9.54 9.54 9.54 9.54 10.2
0 1 2 3 4
-0.156 -0.047 0.07 0.22
Symmetry of the nuclear geometry. For n ) 1-3, the symmetry rotation axis is along the z-axis. For n ) 2, the vector connecting the two Cl atoms is parallel to the x-axis. b Calculated mean polarizability. c Experimental mean polarizability. a
To generalize to liquid-state simulations, one must consider the total potential energy for N molecules N
E)
Eimol + ECoulomb + Enon-Coulomb ∑ i)1
(5)
where Eimol denotes the internal energy for each molecule (see eq 1), ECoulomb represents all the intermolecular Coulomb interactions, and Enon-Coulomb represents the short-range (typically Lennard-Jones) intermolecular interactions. If we apply the equalization principle (to the total energy and imposing electroneutrality for each single molecule), the following set of N-coupled equations is obtained:39
JQl ) ∆K + ∆Vl, l ) 1, ..., N
(6)
Note that in the rigid molecule approach that we are using, the J matrix is common for all molecules. Only the contribution ∆Vl ) (0, V1l - V2l, ..., V1l - VMl) is molecule dependent (Vil is the electrostatic potential due to all other molecules acting on atomic site i of molecule l). We have chosen to solve the N sets of equations iteratively39 at each MD time step: the charges from a previous configuration (or the gas-phase charges for the initial configuration) are used to compute ∆Vl, which produce N new charge sets, when the set of eq 6 is solved. This procedure is iterated until an accuracy of 10-6 for all charges has been reached, which requires about 12 iterations. Using the CPE method, it is relatively easy (see below) to adjust the parameters to molecular properties of interest. In particular, as we are interested in polarization effects, it is capital to have the right anisotropic polarizabilities (see eq 4). While the mean polarizabilities for the molecules considered in this work are known experimentally,47 the polar species have anisotropic polarizabilities that, to our knowledge, have not been determined experimentally. To estimate their value, we have undertaken ab initio calculations which we now describe. B. Ab Initio Polarizabilities. The polarizability tensors, R, for C(CH3)4-nCln are computed42 at ab initio MP2/aug-cc-pVDZ level (i.e., the second-order Møller-Plesset perturbation theory43 with correlation consistent basis sets of valence double-ζ plus polarization quality augmented with diffuse functions44). The electronic energies are calculated45 under uniform electrostatic fields of (0.001 and (0.002 au, followed by a numerical differentiation to determine each component of R. The nuclear geometries are from our previous optimization calculations3 at RHF/6-31G* level. The results are summarized in Table 1 together with experimentally suggested values for comparison. The calculated results show a systematic increase along n which is not observed in the experiments.46 Nonetheless, overall agreement seems quite reasonable, so we proceed to determine a scaling factor of R j (expt)/R j (calc), which is then applied to each
qCl
qC
-0.296 -0.177 -0.03 0.174
0.624 0.437 0.214 -0.13 -0.696
of the calculated components Rxx, Ryy, and Rzz. The polarizability tensors thus obtained are used as a reference in the parameter determination of the charge equilibration method. C. Methyl Group Parametrization. From eq 1, we see that there are several parameters to be fitted in the CPE method: atomic electronegativities (χ0i ) and hardness (J0i ), together with the overlap integrals, Jij(r). The latter are usually computed38 from Slater functions (Nnrn-1 e-ξir). Therefore, we have three parameters for each atomic site (χ0i , J0i , ξi) that are to be fitted to several relevant molecular electrostatic properties: basically, gas-phase partial charges and molecular anisotropic polarizabilities. In our previous work on this family of liquids,3 we detailed the procedure followed to determine the charges (see Table 2) from ab initio calculations (using the united atom approximation for the methyl group). Basically, the charges were fitted to reproduce dipole moments (with excellent accord between computed and experimental values) and electrostatic potential surfaces. As described in the previous section, (scaled) anisotropic polarizabilities have also been computed. The analysis of the n ) 4 case (CCl4) highlights the issues involved in the process of parameter fitting to the computed electrostatic properties. In this case, no parametrization should be required in principle, as Rappe´ et al.38 have determined a broad set of parameters that include chlorine and carbon, which we reproduce here for completeness:
χ0C ) 123.153 kcal/(mol e) χ0Cl ) 197.396 kcal/(mol e) (7) J0C ) 233.399 kcal/(mol e2) J0Cl ) 228.005 kcal/(mol e2) (8) ξC ) 0.8563 au of inverse length
ξCl ) 0.9154 au (9)
First of all, we can compute the polarizability from eq 4 (in this case, because of the molecular symmetry, all components are identical). We obtain a value of 11.1 Å3, which is remarkably close to the experimentally determined values of 10.5-11.5.47 This suggests that the values of J0i and ξi are perfectly acceptable (note that the polarizability tensor does not depend on the electronegativities, χ0i , or equivalently, on the partial charges). Now, we can proceed to the computation of the partial charges from minimization of eq 1, subject to the additional constraint of charge neutrality. In this case, we get simple analytical formulas, which we particularize for the chlorine atomic site:
qCl )
χ0Cl - χ0C 8JC-Cl(r) - 4J0C - J0Cl - 3JCl-Cl(r)
(10)
Surprisingly, we get the wrong sign for the partial charges. This problem was also found by Iskowitz and Berkowitz48 for CH4 and should be attributed to the molecular environment dependence of the CPE parameters, particularly of the electronegativities. It is easy to demonstrate that the denominator
7786 J. Phys. Chem. B, Vol. 105, No. 32, 2001
Llanta et al.
TABLE 3: Polarizability Tensor Components (in kcal/(mol e2))a n
0
1
2
3
method
scaled ab initio
fitted
scaled ab initio
fitted
scaled ab initio
fitted
scaled ab initio
fitted
Rxx
10.2
11.93
Rzz
10.2
R j
10.2
11.24 (5%) 10.43 (4%) 10.6 (1.7%) 10.75 (1.3%)
11.16
10.2
10.33 (13%) 10.33 (13%) 10.98 (19%) 10.54 (15%)
11.91
Ryy
10.24 (0.4%) 10.24 (0.4%) 10.24 (0.4%) 10.24 (0.4%)
11.16 (0) 11.16 (0%) 10.48 (7%) 10.93 (2%)
11.93 13.62 12.5
10. 10.78 10.9
11.16 9.76 10.7
a
Absolute PerCent Difference between Fitted and Scaled ab Initio Values are given in Parentheses.
in eq 10 is negative as a consequence of the minimum condition. Consequently, only the numerator can be fixed to get the right partial charges. In particular, we need to set χ0Cl - χ0C ) -63.94 kcal/(mol e), in comparison with a value of 74.24 kcal/ (mol e) obtained with the parameters in ref 38. To summarize, we have been able to obtain the right partial charges and polarizabilities in a two-step process: first, from the values of ξi and J0i , we get the correct polarizability tensor; second, we adjust the electronegativity differences to get the right partial charges. To address the rest of the molecules, we note that no parameters are available (in the united atom approximation) for the methyl group. Following the hints from the CCl4 case, we will proceed in two stages. First, we will develop a set of methyl parameters (ξMe and J0Me), valid for all molecules containing methyl, from a best fit to all anisotropic polarizabilities (keeping the same parameters for C and Cl as used for CCl4). Second, electronegativity differences, now dependent on the molecule, will be determined to get the gas-phase partial charges. As discussed in the previous section, we scale all polarizabilities computed ab initio to the experimental mean polarizabilities; in Table 3, we display the values resulting for each molecule. Now, we optimize the parameters for the methyl group minimizing the function
χ2 )
2 ∑n ∑i (Rfitn,i - Rscaled n,i )
(11)
where Rfit n,i denotes the i component of the polarizability tensor (with i ) xx, yy, zz, or mean) of molecule n (n ) 1-4) computed from eq 4, using the fitted atomic parameters. With standard minimization techniques,49 a satisfactory set of polarizabilities is obtained (see Table 3), with the deviation between optimized and scaled polarizabilities being below 5% in all cases, except for n ) 1. For the latter case, there is a deviation of 15% between the computed and experimental mean polarizabilities (19% for the zz component). It is to be noted that the ab initio results do not support the larger experimental polarizability for n ) 1, which does not follow the computed trend (see Table 1). This is the origin of the slightly poorer agreement for n ) 1 and, in our opinion, casts some doubts on the accuracy of the experimental value. Altogether, we regard these results as highly satisfactory and thus assign the following optimized parameters to the methyl group (independently of the molecular environment): J0Me ) 208 kcal/(mol e2) and ξMe ) 0.6 au. This set is physically consistent as it implies a smaller hardness for Me than for C and a substantially larger radius. We now turn to the computation of electronegativity differences. Only the partial charges of choice (see Table 2) and the parameters previously computed are required, with no additional optimization. With use of eq 2, we directly obtain by matrix multiplication the differences χ0i - χ0j , which are summarized
TABLE 4: Electronegativity Differences for Each Molecule (in kcal/(mol e))a
an
n
χMe - χC
0 1 2 3 4
54.97 32.21 6.25 -30.5
χCl - χC 64.82 35.34 -4.32 -63.94
denotes the molecule in the series (CH3)4-nCCln.
in Table 4. We highlight the fact that these differences are not transferable but molecule dependent, although as stressed, their calculation is trivial and thus facilitates its use in different molecular environments. III. Computational Details Geometries and Lennard-Jones parameters were not varied from our previous work, and we refer to ref 3 for a complete summary of the parameters. In all cases, the results refer to a system of 216 molecules enclosed in a cubic box. The equations of motion have been integrated with the leapfrog algorithm50 and a time step of 2 fs for both the equilibration and production runs. Because we presently neglect internal motions, the models consist of rigid molecules. The corresponding geometrical constraints were satisfied with the “shake” algorithm,51 with a relative accuracy of 10-5. For the intermolecular interaction, a spherical truncation scheme between molecular centers was taken with a cutoff value of half the box length. For electrostatic forces, the Ewald summation method52 was used, with R ) 5/L and nmax2 ) 16. NVT conditions were maintained during the simulations. To this end, the temperature was regulated through the weak coupling algorithm of Postma et al.50 with a coupling constant of 0.1 ps (we have checked in sample NVE runs that this has no effect on the computed correlation functions). For each pure liquid, we have simulated three different temperatures, which expand the full liquid-phase range at 1 bar of pressure. The experimental densities have been used in all cases. Actually, the nonpolarizable models reproduce these experimental values within statistical indeterminacy for the whole temperature range studied.3 We have also checked in sample cases that the inclusion of polarizability does not produce significant changes of the density values (a 0.5% increase for the nonpolar case n ) 4 and a 2.7% increase for the most polar case n ) 2). For each state point, an independent configuration was generated and equilibrated for 50 ps, after which it was continued during a production run of 200 ps. This same procedure was followed both without considering polarization and with including polarization effects. IV. Results The behavior found depends critically on the polarity of the molecules. The most important effects are found for polar
Polarization Effects in Chlorinated Organic Liquids
J. Phys. Chem. B, Vol. 105, No. 32, 2001 7787 TABLE 5: Gas-Phase Charges and the Corresponding Mean Values in the Liquid Statea n
T (K)
center
gas-phase value
mean liquid-phase value
4
298
C Cl
-0.696 0.174
-0.698 (0.3%) 0.175 (0.6%)
3
298
C Me Cl
-0.13 0.22 -0.03
-0.13 (0%) 0.252 (14%) -0.04 (33%)
2
293
C Me Cl
0.214 0.07 -0.177
0.215 (0.5%) 0.099 (43%) -0.207 (17%)
1
288
C Cl Me
0.437 -0.296 -0.047
0.439 (0.5%) -0.341 (15%) -0.033 (30%)
0
273
C Me
0.624 -0.156
0.624 (0%) -0.156 (0%)
a
Figure 1. Charge distribution obtained with the polarizable model (a) for n ) 2 at T ) 293 K and (b) for n ) 4 at T ) 298 K: dashed line, carbon; thin solid line, methyl; thick solid line, chlorine. The bars denote the corresponding gas-phase charges. .
molecules (n ) 1-3), while considerably smaller effects are found for the nonpolar cases (n ) 0, 4). Within each group, the trends are rather similar, and thus, in the discussion to follow, we will focus on n ) 2 and n ) 0 as representative of polar and nonpolar molecules, respectively. We start discussing the results for charges and dipole moments, which we feel provide considerable insight into the nature of the changes that polarization brings about. A. Induced Dipoles and Charge Distributions. As already emphasized, during the simulations that include polarization, charges are allowed to fluctuate. In Figure 1a, we display the charge distributions obtained for C(CH3)2Cl2. With the bars denoting the charges in the corresponding nonpolarizable models (gas-phase charges), we note that the carbon center is barely affected by the environment, with a sharp and narrow distribution of charges. For chlorine and methyl, we get broader distributions, with a full width of roughly 0.2 e. It should be noted that the hardness is similar for all centers, which suggests that the substantially more peaked result for carbon is due to electrostatic shielding from the outer centers. Basically, charge is transferred from the methyl group to the chlorine centers when the molecule is solvated. In this way, the methyl group turns slightly more positive while the chlorine atom gets more negative, reflecting the higher electronegativity of chlorine (see Table 4). This behavior is considerably less marked for the nonpolar molecules (see Figure 1b for n ) 0).
The percent deviation in absolute value is given in parentheses.
In Table 5, we summarize the mean values obtained from the charge distributions. Regarding the carbon center, the difference between gas-phase and liquid-phase charges does not exceed 0.5% for any molecule, confirming the small effect that we anticipated from the distributions. For the outer centers, the situation is completely different. For n ) 2, which we have just discussed, the methyl charge is increased by a remarkable 43%, while the chlorine charge increases (in absolute value) by 17%. Nevertheless, one cannot draw the conclusion that the methyl group is subject in all cases to the largest variations. While it is true that for n ) 1, the previous numbers only change to 30% and 15%, respectively, for n ) 3, the situation is reversed and the numbers for methyl and chlorine are 14% and 33%, respectively. For the nonpolar molecules, the qualitative behavior is different, as the charge spreading is minimal and largely centered on the gas-phase values. The maximum deviation between gas-phase and mean liquid-phase charge in the nonpolar molecules is 0.6% (for the chlorine center in n ) 4). We now turn to the analysis of dipole moment. In Figure 2a, we display the distribution obtained for the modulus of the total dipole moment for n ) 2, which again represents quite faithfully the general behavior. There is a substantial spreading of dipole moments, with a full width of the distribution of roughly 2 D (an order of magnitude smaller for the nonpolar molecules, see Figure 2b). Furthermore, and as previously discussed for the charges, there is a noticeable shift from the gas-phase value. From Table 6, this shift is about 30% for the polar molecules. If we look at the absolute value of this increase, it roughly spans the range 0.4-0.7 D for the polar molecules and is considerably smaller for the nonpolar molecules (0.06 D for n ) 0 and 0.16 D for n ) 4). The large increases in dipole moment found for n ) 2-3 exceed what has been found for alcohols by Gao et al.,36 increases of roughly 0.4 D (it is to be considered that in alcohols there is hydrogen bonding, an effect that is absent here). Moreover, the computed values contrast with the experimental estimations in the liquid phase which, except for carbon tetrachloride, produce dipole moments that do not differ significantly from the gas-phase values (see Table 6). The present calculations suggest that the liquid-phase dipole moments had been previously underestimated. Actually, this sort of mismatch is not uncommon; a similar behavior was recently found from ab initio calculations for crystalline dimethylnitramine.53 An embedded cluster model yielded a dipole moment 40% greater than the gas-phase value and 15% greater than the estimated dipole moment in the liquid phase. The case of liquid
7788 J. Phys. Chem. B, Vol. 105, No. 32, 2001
Figure 2. Distribution of the dipole moment modulus (a) for n ) 2 at T ) 293 K and (b) for n ) 4 at T ) 298 K. The bar corresponds to the gas-phase value.
TABLE 6: Gas-Phase Dipole Moments (ab Initio3 and Experimental) and Liquid-Phase Dipole Moments (Experimental and Computed Mean Values) (in debye)e
n
T (K)
gas-phase value (exptl)
4 3 2 1 0
298 298 293 288 273
0 1.79b 2.25b 2.14c 0
gas-phase value (ab initio)
liquid-phase value (exptl)
mean liquid-phase value (computed)
0 1.861 2.347 2.294 0
0.15a 1.78b 2.33b 2.16d
0.164 2.22 (24%) 2.918 (30%) 2.81 (31%) 0.064
a Reference 54. b Reference 55. c Reference 56. d Reference 57. e The percent deviation, in absolute value, between computed and gas-phase values is given in parentheses.
water is another interesting example. While the liquid-phase dipole moment estimated experimentally is ≈2.6 D58 and the most popular fixed charge model (SPC/E59) uses a dipole moment of 2.35 D, recent ab initio molecular dynamics simulations60 have shown that it is probably close to 3.0 D (60% larger than the gas-phase value and 15% larger than the estimated liquid-phase value). A result which suggested this mismatch (2.8 D) had previously been obtained with a fluctuating charge model,16 similar to the one used here. Therefore, while our results are consistent with the trends found in other systems, it is certainly surprising that for these non-hydrogenbonded liquids, the effect is even slightly more marked.
Llanta et al.
Figure 3. Radial distribution functions (a) for n ) 0 and (b) for n ) 2: thick solid line, gC-C; thin solid line, gMe-Me; dashed line, gCl-Cl. Only the lines for the polarizable model are included because they coincide almost exactly with those of the nonpolarizable model.
B. Structural Properties. We limit the structural study to radial distribution functions. As expected, mostly negligible effects are found in all cases, as is evident in Figure 3, where we display a few sample cases. Figure 3a shows gC-C(r) and gMe-Me for C(CH3)4 for both the polarizable and nonpolarizable models. Clearly, the inclusion of polarization does not affect structure to any significant extent (as is also the case for CCl4). A related result was already obtained by switching off the charges in the nonpolarizable models,3 which had no effect on structure. Figure 3b displays the corresponding results for n ) 2. Again, no noticeable differences are found here, in contrast with a calculation with no charges previously reported.3 There, it was found that the neglect of charges gives rise to almost identical radial distribution functions for the pairs Cl-Cl, ClMe, and Me-Me, consistent with the fact that these centers only have small differences in Lennard-Jones parameters. Consequently, for the polar molecules, while mean electrostatic effects (represented by fixed charge models) are nonnegligible for the structure, the addition of polarizability does not result in any significant change. C. Dynamic Properties. 1. Translational Diffusion. We have computed diffusion coefficients both from the integration of the center of mass velocity autocorrelation function and from the mean square displacement, with no differences found between the methods. Numerical results are summarized in Table 7, although the general trends are probably better appreciated in
Polarization Effects in Chlorinated Organic Liquids
J. Phys. Chem. B, Vol. 105, No. 32, 2001 7789
TABLE 7: Diffusion Coefficients n
T (K)
D (nonpolarizable)
D (polarizable)
experimental
4
273 298 313
1.24 1.71 2.05
1.25 1.77 1.94
1.284a (at 293 K)
3
273 298 313
1.46 2.03 2.35
1.38 1.92 2.34
2
273 293 313
1.67 2.09 2.81
1.57 2.04 2.76
1
273 288 313
1.83 2.22 3.01
1.77 2.13 2.99
260 273 283
3.40 4.08 4.58
3.36 4.07 4.53
0
a
2.2b
Reference 61. b Reference 62.
Figure 4. Diffusion coefficients as a function of temperature. n indicates the number of Cl atoms in the molecule. Circles correspond to the nonpolarizable model and triangles to the polarizable one. Continuous and dashed lines show the trend for each model.
Figure 5. Translational velocity power spectra (a) for n ) 2 at T ) 293 K and (b) for n ) 4 at T ) 298 K. Continuous (dashed) lines correspond to the nonpolarizable (polarizable) model.
Figure 4. In all cases, there is little or no difference between polarizable and nonpolarizable models, with the polarizable ones yielding slightly smaller diffusion coefficients. This behavior is common to all liquids (from Figure 4, it might be interpreted for CCl4 that there exists a faster dynamics in the polarizable case at 298 K, but we believe that the differences are not statistically significant). We have also analyzed the frequency dependence of translational dynamics, looking for the signature of polarization effects. To this end, we have computed the frequency spectrum for the velocities from the Fourier transform of the velocity autocorrelation functions
C ˆ (ω) ≡
∫0∞〈Vb(t)‚Vb(0)〉 dt
(12)
Figure 5 displays the results for n ) 2 and n ) 4 for polarizable and nonpolarizable models. Again, no noticeable differences are found; in all cases, the spectra are peaked at about 15 cm-1 with no additional features. Regarding comparison with experimental results, we have only been able to find some estimations of the diffusion coefficients for n ) 4 and n ) 1. In the first case, there is pretty good accord; the difference between computed (1.7 × 10-9 m2
Figure 6. Rotational relaxation functions for n ) 2 at T ) 293 K: solid line, no polarization; dashed line, polarization included.
s-1) and experimental (1.3 × 10-9 m2 s-1) can be partly attributed to the different temperatures. For n ) 1, the accord is even better and we have a numerical result of 2.2 × 10-9 m2 s-1 to be compared with an experimental value of also 2.2 × 10-9 m2 s-1. We conclude that both the polarizable and
7790 J. Phys. Chem. B, Vol. 105, No. 32, 2001
Llanta et al.
TABLE 8: Rotational Relaxation Constants (in ps)d n
T (K)
τ1 (nonpolarizable)
τ1 (polarizable)
4
273 298 313
6.04 4.93 4.31
6.14 4.96 4.31
273 298 313
5.76 4.43 4.05
6.25 4.77 4.06
273 293 313
5.15 4.17 3.52
5.54 4.76 3.82
1
273 288 313
4.39 3.76 3.02
0
260 273 283
1.62 1.45 1.44
3
2
a
τ1 (experimental)
τ2 (nonpolarizable)
τ2 (polarizable)
2.72 2.15 1.91
2.72 2.24 2.08
4.9b
2.71 2.01 1.87
2.87 2.19 1.93
2.05b
3.5c
2.09 1.78 1.44
2.24 1.96 1.55
1.53c
4.81 4.00 3.18
2.23 1.82 1.37
2.24 1.96 1.53
1.86 1.51 1.33
0.90 0.94 1.17
1.19 0.91 0.86
4.374a (293 K)
τ2 (experimental) 1.689a (293 K)
Reference 61. b Reference 55. c Reference 63. d Results with polarizable models are given in parentheses.
nonpolarizable models are indistinguishable regarding translational diffusion with very good agreement with experimental results. 2. Rotational Relaxation. The study of rotational relaxation is usually limited to the computation of the correlation functions
c1(t) ) 〈u b(t)‚u b(t)〉
(13)
and
c2(t) )
〈21[3(ub(t)‚ub(0)) - 1]〉 2
(14)
where b u is a unitary vector along a given molecular direction (here, we have checked that the relaxation is independent of the axis and thus have used in all cases a molecular symmetry axis). Usually, both functions decay exponentially, which is also the case to a good approximation for all systems studied here. In Figure 6, we display c1(t) and c2(t) for n ) 2 for only one sample temperature. From these functions, it is possible to compute the relaxation time by adjusting an exponential or by integration of the full correlation function. Here, we have adjusted exponential functions, skipping the short time Gaussian decay and starting the fit at 2 ps. Table 8 summarizes the relaxation times (τ1 and τ2) for each system and temperature, with and without polarization, together with experimental estimations. As was the case for the diffusion coefficient, no significant differences result from the inclusion of polarization. The main effect is a slight slowing down of the rotational relaxation, which translates into a slight increase of the correlation times. In Figure 7, we display the times obtained for n ) 2 at all temperatures. The dynamics is faster with increasing temperature, as expected, but the effect of polarization is largely temperature independent. As it was the case with diffusion, we conclude again that polarization is rather negligible for the rotational dynamics. Regarding the comparison with experimental estimations, the accord is quite satisfactory in all cases for which data are available, considering that experimental fits are also subject to some indeterminacy. Possibly, there is a slightly slower relaxation found in the molecular dynamics results compared with the experimental estimations. V. Concluding Remarks An exhaustive study of dynamics properties and the possible influence of polarization at different temperatures has been
Figure 7. Rotational relaxation times for n ) 2 at different temperatures. Circles correspond to the nonpolarizable model and triangles to the polarizable one. Continuous and dashed lines show the trend for each model.
performed for the methylchloromethane family of organic liquids. A considerable effort has been invested in the development of a largely transferable set of parameters for the methyl group to be used within the chemical equalization method. Only electronegativities are left as environment-dependent parameters, although they can be trivially computed from partial charges (which, in turn, can be obtained from a fit to electrostatic potential surfaces determined ab initio). The overall effect of polarization is minimal, judging from its effect on translation, rotation, and radial distribution functions. In contrast, regarding induction effects, numerical computations suggest that there exist substantial induced dipole moments for polar molecules, exceeding by about 30% what had been estimated up to now, in line with recent findings for other organic molecules. These induction effects result from charge flow between the outer centers within the molecule with a substantial shift of mean values from gas-phase charges. For nonpolar molecules, the induced dipoles are smaller by roughly 1 order of magnitude and the charge distributions are symmetrical around gas-phase charges. We anticipate that polarization will thus be critical for the dielectric properties, which we intend to study soon. The method has already proved to be
Polarization Effects in Chlorinated Organic Liquids successful for the computation of the absorption coefficient of carbon tetrachloride.64 Acknowledgment. This work was supported by DGICYT Project PB-96-0170-C03-02. References and Notes (1) Schwarzenbach, R. P.; Gschwend, P. M.; Imboden, D. M. EnVironmental Organic Chemistry; Wiley: New York, 1995. (2) American Water Works Association. Water Quality and Treatment, 4th ed.; McGraw-Hill: New York, 1990. (3) Rey, R.; Pardo, L. C.; Llanta, E.; Ando, K.; Lo´pez, D. O.; Tamarit, J. Ll.; Barrio, M. J. Chem. Phys. 2000, 112, 7505. (4) Sprik, M. J. Phys. Chem. 1991, 95, 2283; J. Chem. Phys. 1991, 95, 6762. (5) Barnes, P.; Finney, J. L.; Nicholas J. D.; Quinn, J. E. Nature 1979, 282, 459. (6) Sprik, M.; Klein, M. L. J. Chem. Phys. 1988, 89, 7556. (7) Rullman, J. A. C.; van Duijnen, P. T. Mol. Phys. 1988, 63, 451. (8) Ahlstro¨m, P.; Wallqvist, A.; Engstro¨m, S.; Jo¨nsson, B. Mol. Phys. 1989, 68, 563. (9) Zhu, S.-B.; Singh, S.; Robinson, G. W. J. Chem. Phys. 1991, 95, 2791. (10) van Belle, D.; Froeyen, M.; Lippens, G.; Wodak, S. J. Mol. Phys. 1992, 77, 239. (11) Caldwell, J.; Dang L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 9144. (12) Dang, L. X. J. Chem. Phys. 1992, 97, 2659. (13) Niesar, U.; Corongiu, G.; Clementi, E.; Kneller, G. R.; Bhattacharya, D. K. J. Phys. Chem. 1990, 94, 7949. (14) Wallqvist, A.; Berne, B. J. J. Phys. Chem. 1993, 97, 13841. (15) Bernardo, D. N.; Ding, Y.; Krogh-Jespersen, K.; Levy, R. M. J. Phys. Chem. 1994, 98, 4180. (16) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (17) Borgis, D.; Staib, A. Chem. Phys. Lett. 1995, 238, 187. (18) Chialvo, A. A.; Cummings, P. T. J. Chem. Phys. 1996, 105, 8274. (19) Dang, L. X.; Chang, T.-S. J. Chem. Phys. 1997, 106, 8149. (20) Jedlovszky, P.; Richardi, J. J. Chem. Phys. 1999, 110, 8019. (21) Kiyohara, K.; Gubbins, K. E.; Panagiotopoulos, Z. Mol. Phys. 1998, 94, 803. (22) Yezdimer, E. M.; Cummings, P. T. Mol. Phys. 1999, 97, 993. (23) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2378. (24) Chen, B.; Xing, J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2391. (25) Warshel, A. J. Phys. Chem. 1979, 83, 1640. (26) van Belle, D.; Couplet, I.; Prevost, M.; Wodak, S. J. J. Mol. Biol. 1987, 198, 721. (27) Kuwajima, S.; Warshel, A. J. Phys. Chem. 1990, 94, 460. (28) Straatsma, T. P.; McCammon, J. A. Chem. Phys. Lett. 1991, 177, 433. (29) Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. 1995, 99, 6208. (30) Bursulaya, B. D.; Zichi, D. A.; Kim, H. J. J. Phys. Chem. 1996, 100, 1392. (31) Meng, E. C.; Kollman, P. A. J. Phys. Chem. 1996, 100, 11460. (32) Stuart, S. J.; Berne, B. J. J. Phys. Chem. 1996, 100, 11934.
J. Phys. Chem. B, Vol. 105, No. 32, 2001 7791 (33) Young, W. S.; Brooks, C. L., III. J. Chem. Phys. 1997, 106, 9265. (34) Rick, S. W.; Berne, B. J. J. Phys. Chem. B 1997, 101, 10488. (35) Bryce, R. A.; Vincent, M. A.; Malcolm, N. O. J.; Hillier, I. H.; Burton, N. A. J. Chem. Phys. 1998, 109, 3077. (36) Gao, J.; Habibollazadeh, D.; Shao, L. J. Phys. Chem. 1995, 99, 16460. (37) Gonza´lez, M. A.; Enciso, E.; Bermejo, F. J.; Be´e, M. J. Chem. Phys. 1999, 110, 8045. (38) Rappe´, A. K.; Goddard, W. A. J. Phys. Chem. 1991, 95, 3358. (39) Chelli, R.; Ciabatti, S.; Cardini, G.; Righini, R.; Procacci, P. J. Chem. Phys. 1999, 111, 4218. (40) Stockelmann, E.; Hentschke, R. J. Chem. Phys. 1999, 110, 12097. (41) Ribeiro, M. C. C.; Almeida, L. C. J. J. Chem. Phys. 1999, 110, 11445. (42) We used the program GAMESS for the quantum chemistry calculations: Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. J.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (43) Møller, C.; Plesset, M. S. Phys. ReV. 1934, 45, 618. Krishnan, R.; Frisch, M. J.; Pople, J. A. J. Chem. Phys. 1980, 72, 4244. (44) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1995, 103, 4572. (45) The accuracy of the electronic wave functions has been increased by reducing the numerical thresholds (from defaults) as follows: the primitive cutoff factor is 10-30 au, the integral cutoff factor is 10-20 au, and the electron density convergence is 10-10 au. (46) We should note that the experimental methods are not common among n ) 0-4. In particular, the experimental data for n ) 1 and 2 are evaluated from the molar refraction index which thus corresponds to the polarizability at optical frequencies. (47) Lide, D. R. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1997. (48) Iskowitz, P.; Berkowitz, M. L. J. Phys. Chem. A 1997, 101, 5687. (49) Press, W. H.; Flannery, B. P.; Teukolski, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992. (50) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (51) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (52) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (53) Johnson, M. A.; Truong, T. N. J. Phys. Chem. B 1999, 103, 9392. (54) Vij, J. K. J. Chem. Phys. 1987, 87, 3357. (55) Clemett, C.; Davies, M. Trans. Faraday Soc. 1962, 58, 1705. (56) Morgan, S. O.; Yager, W. A. Ind. Eng. Chem. 1940, 32, 1519. (57) Chen, J.; Bartell, L. S. J. Phys. Chem. 1993, 97, 10645. (58) Coulson, C. A.; Eisenberg, D. Proc. R. Soc. London, Ser. A 1966, 291, 445. (59) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Chem. Phys. 1987, 91, 6269. (60) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 3572. (61) Chahid, A.; Bermejo, F. J.; Enciso, E.; Garcia-Herna´ndez, M.; Martınez, J. L. J. Phys.: Condens. Matter 1992, 4, 1213. (62) Kessler, D.; Weiss, A.; Witte, H. Ber. Bunsen-Ges. Phys. Chem. 1967, 71, 3. (63) Antony, A. A.; Smyth, C. P. J. Am. Chem. Soc. 1964, 86, 152. (64) Llanta, E.; Rey, R. Chem. Phys. Lett., in press.
