17956
J. Phys. Chem. 1995, 99, 17956-17962
Dielectric Constants of Formamide and Dimethylformamide via Computer Simulation Jonathan W. Essex? and William L. Jorgensen* Department of Chemistry, Yale University, New Haven, Connecticut 06520-8107 Received: May 30, 1995; In Final Form: August 30, 199.5@
Monte Carlo simulations of liquid formamide and dimethylformamide, using the united-atom OPLS potential energy functions, have been undertaken. Basic thermodynamic properties and radial distribution functions are largely unaffected by the inclusion of long-range electrostatic interactions in the form of a reaction field, whereas the dipole-dipole correlation function is observed to be very sensitive. Static dielectric constants of the liquids were calculated by measuring the polarization of the systems as a function of an applied electric field. The dielectric constant calculated for dimethylformamide is in good agreement with experiment, whereas the computed value for formamide is too low. It is likely that treatment of explicit polarization is needed for hydrogen-bonded liquids and related systems including proteins to correctly describe dielectric behavior. The dependence of the calculated dielectric constants on system size and reaction field dielectric constant was also examined.
Introduction The dielectric constant of a medium plays a critical role in determining chemical behavior and, in particular, in moderating the strength of ion-ion interactions.' Therefore, it is important that the empirical potential energy functions used in computer simulations of fluids be tested for their ability to reproduce experimental dielectric constants. In this paper the dielectric constants of liquid formamide and dimethylformamide (DMF), as described by the united-atom OPLS potential energy functions,2 are calculated using Monte Carlo simulations. These liquids are of particular interest for two reasons. First, they represent two important classes of solvents, hydrogen-bonded and dipolar aprotic, and although the gas-phase dipole moments of the molecules are very similar, their dielectric constants are markedly differ en^^,^ The dipole moments of fonnamide and DMF are 3.7 and 3.9 D, respectively,2 whereas the dielectric constants of these liquids are approximately 109 and 37.3 Reproducing this pattem constitutes a challenging test of the potential energy functions. It is unclear whether the commonly used fixed-charge models can simultaneously yield correct densities, heats of vaporization, and such a high dielectric constant as for formamide. Second, these molecules can be regarded as temnlates for protein backbones; it seems important that results for pure liquid amides should reproduce experimental dielectric behavior if one is to have confidence in the performance of such potential functions in a protein environment. This is particularly true in light of recent simulation results suggesting that local dielectric constants in proteins are ca. 10-30.'34 This study also addresses the dependence of thermodynamic and structural properties of pure liquids on the treatment of longrange electrostatic interactions.
Methodology It is accepted that the calculation of liquid properties dependent on intermolecular angular correlations, via computer simulation, requires the inclusion of long-range electrostatic interactions; the static dielectric constant is one such There are two common methods for evaluating long-range Current address: Department of Chemistry, University of Southampton, Highfield, Southampton, SO17 1BJ United Kingdom. * To whom correspondence should be addressed. Abstract published in Advance ACS Absrracrs, December 1, 1995. @
0022-365419512099-17956$09.0010
electrostatic interactions. In the Ewald summation approach, electrostatic interactions are evaluated through an infinite lattice sum over all periodic images, thereby imposing pseudocrystalline order on the l i q ~ i d . ~In, ~the . ~reaction field approa~h,5-~.'~*'' electrostatic interactions beyond the spherical cutoff, r,, are approximated by treating the system in excess of r, as a polarizable continuum of dielectric constant cd. Each molecule of the liquid experiences an electric field arising from the polarization of the continuum by the net dipole moment within the cutoff sphere centered on that molecule. This approach therefore assumes that the liquid is structureless beyond r,, so that a continuum approximation can be applied. It has been demonstrated that the Ewald sum and reaction field methods yield indistinguishable dielectric constants for SPC water, but through markedly different orientational behavior.I2 In the work presented here, the reaction field approach was adopted because of its conceptual simplicity and computational efficiency. The nonphysical replication of the instantaneous dipole for the central cell in the Ewald method reduces the attractiveness of that approach. The reaction field is applied by the addition of a single term to the Coulombic potential for the interaction between site i and all sites j in other molecules within rc.13
qi and % represent the charges on atomic sites i and j , respectively, and r,j is the distance between the atomic sites. This expression is valid provided that the molecules of the liquid are electrically neutral and that the nonbonded cutoff is applied with respect to the molecule as a whole. It is apparent from eq 1 that a value of cd has to be selected. Ideally, cd should be identical to the dielectric constant of the liquid within r,, errif discontinuities at the cutoff boundary are to be avoided. This does, however, require knowledge of the property the simulation is intended to determine. In practice, however, any choice of cd such that cr 5 5 -has been demonstrated to yield consistent values of cr for moderate to highly polar l i q ~ i d s . ' ~ , ' ~ This observation is confirmed in the work presented here. A number of methods exist for the calculation of dielectric constants through computer simulation. Arguably the most widely used approach is to calculate the average of the square 0 1995 American Chemical Society
Dielectric Constants of Formamide and Dimethylformamide of the total dipole moment of the system, (M2). For reaction field boundary conditions, eq 2 is then applied to give ere
V is the volume of the simulation system, kb is the Boltzmann constant, and T the temperature. This approach is referred to subsequently as the fluctuation method. It has been shown that markedly different values of (M2) arise from different choices of e, but that identical values of Er are obtained through the consistent application of eq 2.'33'59'6 In practice, however, the choice of e, x er yields the most rapid convergence of (M2).'63'7 A recent extension of this general method involved the application of perturbation theory to obtain values of E, for which E , = Er*'6 Another commonly used approach for the calculation of dielectric constants is to determine the polarization response of the liquid to an applied electric field, &.I8 If the average system dipole moment along the direction of the applied field, per unit volume, is (P), then cr can be cakulated using eq 3 for reaction field g e ~ m e t r y . ' ~
This approach is between 2 and 4 times more efficient than the fluctuation method for SPC water.I9 It does, however, require the perturbation of the liquid structure through the application of an electric field. In practice, therefore, care must be taken to ensure that the observed polarization, (P), is linear in the applied field strength, EQ;i.e., dielectric saturation is a~oided.'~,~O A perturbation theory based on the change in free energy of the liquid on application of an applied field has recently been derived. However, the method is less efficient than the conventional polarization response approach outlined above.2 Other methods for calculating the static dielectric constants of liquids include the use of ion-ion potentials of mean force,22 and an umbrella sampling approach whereby the complete probability distribution of the net dipole moment, P ( M ) , is ~ a l c u l a t e d .A ~ ~related approach, which also yields the necessary (M2) for eq 2, has been reported and involves computing P(M2) explicitly for only small ranges of M2.24 For a softsphere plus imbedded dipole fluid, this procedure is reported to be 5 times more efficient than the standard fluctuation method. Comparisons with the efficiency of the polarization method are desirable for a molecular system. It may also be noted that the full frequency dependence of the dielectric constant can be computed via molecular dynamics simulations with a fluctuation-based procedure, as done recently for a polarizable model of ~ a t e r . ~ 5 In the work presented here, the polarization response method was adopted because of its reported effi~iency.'~ Furthermore, this work could, in principle, be extended to include the dielectric saturation behavior of the amides in an identical manner to that described in ref 20. Liquid Simulations Monte Carlo statistical mechanics simulations were initially performed on the pure liquids to determine the sensitivity of their thermodynamic properties and radial distribution functions to system size, cfi, and simulation ensemble. Very long simulations (150-300 million configurations) in conjunction with an applied field, Eo, were then performed to determine cr. An applied field strength of 1.5 x lo8V m-l was selected for these studies, and all simulations were performed under constant NVT conditions. This value of &,has been shown to lie within the linear region of polarization response for SPC ~ a t e r , 'and ~ . it ~~
J. Phys. Chem., Vol. 99, No. 51, 1995 17957 is demonstrated in this paper that this is also true for OPLS formamide and DMF. The OPLS force field and potential energy parameters for formamide and DMF, as described in ref 2, were implemented in the pure-liquid Monte Carlo program MCLIQ.26 Cubic simulation systems consisting of either 128 or 256 molecules were constructed and extensively equilibrated for at least 5 million configurations; the magnitudes of molecular translation, rotation, and dihedral angle movement were selected to yield an acceptance probability in the course of each Monte Carlo simulation of approximately 40%. All simulations were performed at 25OC with Metropolis ampl ling.^ Where appropriate, a pressure of 1.0 atm-was applied through attempted changes in system volume of up to 400 A3 every 1200 configurations. For simulations performed under constant NVT conditions, the system density was set equal to the experimental value.2 A spherical nonbonded cutoff was applied to truncate intermolecular interactions at distances of 10.0 and 12.0 A for formamide systems of 128 and 256 molecules, respectively, and 12.0 and 15.0 A for DMF systems of 128 and 256 molecules, respectively. The cutoff is calculated based on the separation of carbonyl carbon atoms for formamide and on the nitrogen atom separation for DMF; this convention is identical to that adopted in ref 2. A correction is applied to the interaction energy in the course of the simulations to take account of Lennard-Jones interactions neglected beyond the cutoff.5
Results and Discussion Unpolarized Liquids. The thermodynamic properties computed for pure liquid formamide and DMF with e, = 1.0 (no reaction field) and E , = 00 (conducting boundary conditions) are given in Tables 1 and 2, together with the data reported in the original paper.2 For simulations performed under constant NVT conditions, enthalpies of vaporization (AHvBp) and constantpressure heat capacities (C,) are not given. The statistical uncertainties (flo) reported here were obtained from the fluctuations in separate averages over batches of 0.5 million configuration^.^ The statistical uncertainties reported in ref 2 were estimated as approximately f 0 . 2 A3 molecule-' for the molecular volume and f0.025 kcal mol-' for the enthalpy of vaporization. However, the uncertainties are probably underestimated because of the small block size (0.05 million configurations) used in those calculations. Consider first the results obtained without the presence of a reaction field, for which e, is unity. Comparison of the data reported in ref 2 with that obtained under identical conditions, but over a longer simulation (for formamide, column 1 with column 6, and for DMF column 1 with column 5 ) , reveals differences of only about 1% in the energy and density. For both formamide and DMF, increasing the system size together with the nonbonded cutoff results in a 1-3% decrease in the enthalpy or energy of vaporization (columns 1 and 2). However, the enthalpies of vaporization and the molecular volumes are all close to the experimental data. Comparison of the results obtained under identical simulation conditions, but in the NPT and NVT ensembles (columns 2 and 3), demonstrates that changing the ensemble has no statistically significant effect on the intermolecular interaction energy, Ei. This follows from the very small errors in the computed volumes. Consider now the simulations performed using a reaction field with e, = 00. Application of a reaction field term reduces the magnitude of the intermolecular interaction energy for formamide (columns l and 4 of Table l); this is consistent with previously reported studies on liquid water.I9 Furthermore, the
Essex and Jorgensen
17958 J. Phys. Chem., Vol. 99, No. 51, 1995
TABLE 1: Thermodynamic Properties of Liquid Formamid@ simulation NMOL Erf
ensembJe cutoff1A configs MC TIT -Ed) Emtrd) Eintra(g) AEwp AHvap AH,,, exptd CoP(igld
CP(U Cp(l)exptd V V exptd
1
2
3
128 1.o NPT 10.0 5M 25.0 14.39 f 0.02 0.31 f 0.00 0.30 14.38 f 0.02 14.97 f 0.02 14.7 11.1 20.2 f 1.2 25.7 66.3 i 0.1 66.3
256 1.o NPT 12.0 10 M 25 .O 13.98 f 0.01 0.31 f0.00 0.30 13.97 & 0.01 14.56 f 0.01 14.7 11.1 19.9 f 0.8 25.7 66.7 i 0.1 66.3
256 1.o NVT 12.0 10 M 25.0 14.00 f 0.02 0.31 f 0.00 0.30 14.00 f 0.02 NIA 14.1 11.1 N/A 25.7 (66.3) 66.3
5
6'
128
256
m
W
NVT 10.0 50 M 25.0 14.22 f 0.01' 0.31 f 0.00 0.30 14.21 f 0.01' NIA 14.7 11.1 N/A 25.7 (66.3) 66.3
NVT 12.0 10 M 25.0 13.86 f 0.02c 0.31 f 0.00 0.30 13.85 f 0.02' NIA 14.7 11.1 NIA 25.7 (66.3) 66.3
128 1.o NPT 10.0 1M 25.0 14.20 0.31 0.30 14.19 14.78 14.7 11.1 19.2 25.7 66.8 66.3
4
Energies are given in kcal mol-I; molecular volumes, V , are given in A3 molecule-); specific heat capacities are given in cal mol-' K-I. Simulation data reproduced from ref 2. Includes the self-interaction energy of an individual dipole in its own reaction field. Experimental data reproduced from ref 2.
TABLE 2: Thermodynamic Properties of Liquid Dimethylformamid@ simulation
1
2
3
128 1.o NPT 12.0 5M 25.0 10.78 f 0.02 0.30 f 0.00 0.30 10.77 f 0.02 1 1.36 f 0.02 11.1 21.3 30.4 f 0.9 35.8 128.9 f 0.2 128.6
256 1.o NPT 15.0 10 M 25.0 10.69 f 0.01 0.30 f 0.00 0.30 10.69 f 0.01 11.28 f 0.01 11.1 21.3 31.5 f 1.1 35.8 129.0 f 0.1 128.6
256 1.o NVT 15.0 10 M 25.0 10.70 f 0.01 0.30 f 0.00 0.30 10.70 f 0.01 NIA 11.1 21.3 N/A 35.8 (128.6) 128.6
4
56
NVT 12.0 50 M 25.0 10.80 5 0.00' 0.30 f.0.00 0.30 10.80 f 0.00' NIA 11.1 21.3 NIA 35.8 (128.6) 128.6
128 1.o NPT 12.0 1M 25.0 10.85 0.30 0.30 10.85 11.44 11.1 21.3 32.0 35.8 128.1 128.6
128 W
Energies in kcal mol-'; molecular volumes, V, in A3 molecule-'; specific heat capacities in cal mol-' K-I. Simulation data reproduced from ref 2. Includes the self-interaction enerev of an individual diDole in its own reaction field. Experimental data reproduced from ref 2. e Estimated using the additivity rules given in ref 35.-
interaction energy still exhibits a cutoff dependency (columns 4 and 5 ) . In contrast, the influence of the reaction field on the energetics for DMF is negligible (columns 1 and 4 of Table 2). It should be noted that intermolecular interaction energies reported here include the self-interaction term arising from the interaction of a dipole with its own reaction field,27as given in eq 4. This is just a constant which amounts to -0.078 kcal/ mol for formamide (p = 4.33 D) and to -0.083 kcal/mol for dimethylformamide (p = 4.45 D) when r, = 12 A. (4) In Figures 1 and 2, the computed N-X radial distribution functions (rdfs) are presented for formamide and DMF, respectively; there are only marginal differences in the liquid structure as represented by the rdfs. Indeed, in Figure 3, one region of the N - 0 rdf is shown which suggests that changing the system size, together with the associated cutoff, has a greater effect on the rdf than the application of the reaction field. This observation of rdf invariance to the value of is consistent with other s t ~ d i e s . ' ~ ,Inspection ' ~ ~ ~ ~ ~of~ additional ~ results (data not presented) c o n f i i s that changing simulation conditions through nonbonded cutoff, system size, ensemble, and cd causes only small changes to the rdfs.
In Figures 4 and 5 the dipole-dipole correlation functions, defined in eq 5 , are presented for formamide and DMF. (5) p~ and pug are the dipole moments of molecules A and B
separated by a distance r, and (cos &r)) is the mean value of cos &r), where O(r) is the angle between the dipole moment vectors of molecules A and B. The distance r is calculated on the same basis as the nonbonded cutoff. The data were accumulated in histogram bins of width 0.1 A. Unlike the rdfs, the dipole-dipole correlation functions show significant dependency on simulation conditions. This is consistent with the observed sensitivity of orientational correlation functions in water simulation^.'^*^^.*^ However, it is notable that the results obtained for 256 molecules, in the absence of the reaction field term but in different ensembles, are almost identical, suggesting that changing the simulation ensemble from NET to NVT has little effect on intermolecular orientational structure. The dipole-dipole correlation functions obtained in the presence of the reaction field show discontinuities at the value of the nonbonded cutoff, demonstrating the role of long-range electrostatic interactions in influencing molecular ordering. In each case with the reaction field, the dipole-dipole correlation
Dielectric Constants of Formamide and Dimethylformamide
J. Phys. Chem., Vol. 99, No. 51, 1995 17959
Formamide
Formamide
N-X radial distribution h~~~ctions I
-
2.4
"
'
b
"
"
"
1.8
M
1.6
-
1.4
-
1.2
-
.
1.0
-
0.8
-
ti
N-N A 0
g
V
q ,11, ,
I .o
-03-
-
-04
-
2.0
3.0
,
4.0
,
,
,
,
6.0
5.0
, 7.0
,I 8.0
ji
- 128 molecules. NPT. E~ = 1 0
'i
-- -
___
256 molecules. NFT, E,, = 1 0 256 molecules. NVT, E,, = I 0 128 molecules, NVT, E,, = 256 molecules, NVT, E,, =
--
:
-05
-
0.6
0.0
I
-
2.0
v
"
"
128 molecules. NPT, E,,= 1.0 .......... 256 molecules, NPT, E,, = 1.0 256 molecules, NVT, E,, = 1 .O 128 molecules, NVT, E, = 256 molecules, NVT, E,, = w
2.2
-
Dipolc-dipole correlation function
- 09 -08
00
20
40
60
80
120
LOO
140
r/ A
rlA
Figure 1. Rdfs for liquid formamide. H(C)and H(T) are the hydrogen atoms cis and trans to the oxygen atom.
Figure 4. Dipole-dipole correlation functions for liquid formamide. Dimethylfomamide
Dimethylformamide N-X radial distribution functions 2.4
- 128 molecules, NPT, E,, = 1.0 256 molecules, NPT, E,,= 1.0
-
___ 256.molecules. NVT. E, = 1.0 2.0
128 molecules, NVT. E,, =
1.6 -
- 128 molecules,NPT.
. M 1.2 1.4
v
1.0
-
4 7 46
0.8 1 4 9 4,8
0.6
-1.1
0.2
I
0.0
= 1.0 256 molecules. NPT. E,, = 1.0 256 molecules. NVT. E,, = 1.0 128 molecules, NVT, E, =
-
i i
"
2.0
"
4.0
"
6.0
"
8.0
"
10.0
"
12.0
"
14.0
16.0
r/A 1 .o
2.0
3.0
4.0
5.0
6.0
7.0
8.0
rlA
Figure 2. Rdfs for liquid dimethylfonnamide. C(C) and C(T) are the methyl groups cis and trans to the oxygen atom. Formamide N - 0 radial distribution functions 128 molecules, NPT, E,, = 1.0 256 molecules, NFT, E# = 1 0 256 molecules., NVT.. E-" = 1.0 128 molecules, NVT. &, = w 256 molecules, NVT. &, =
-
~
i
-1.0
0.4
0.0
....
0.6 I 3.0
4.0
5.0
6.0
7.0
8.0
r/A
Figure 3. N-0 rdfs for formamide. function beyond r, mimics the results for the larger system size in the absence of the reaction field.
Figure 5. Dipole-dipole correlation functions for liquid dimethylformamide. Thus, the basic thermodynamic properties and rdfs are largely insensitive to changes in system size, nonbonded cutoff, and ~d. However, orientational behavior is very sensitive to system size and inclusion of the reaction field. Polarized Liquids. As described in the Methodology section, the dielectric constants of the fluids are calculated from the observed polarization as a function of the applied electric field strength. Initially, ed was given a value of corresponding to conducting boundary conditions. In Figure 6, the accumulated average dielectric constant of formamide is presented as a function of simulation length, for a system of 128 molecules, together with the instantaneous value of the dielectric constant obtained in batch averages of 0.5 million configurations. In Table 3, the values of the observed polarization along the direction of the applied field are presented in batch averages of 50 million configurations. eI was calculated from these data using eq 3. The final converged value for the dielectric constant from this simulation was 56.2 f 1.8. The quoted uncertainty corresponds to the standard error obtained over the batch averages of 50 million configurations. To determine the consistency of the results as function of ~ d 300 , million configurations were then executed with ed set to 56.0. In Figure 7 the average dielectric constant is presented as a function of simulation length, together with the instantaneous value obtained over batch averages of 0.5 million configurations. The associ-
a
17960 J. Phys. Chem., Vol. 99, No. 51, I995
Essex and Jorgensen Formamide dielectric constant
Formamide dielectric constant E- = '-'-
i
m8tantaneous e, accumulated e,
-
I 1
8001
-. 128 molecules
e,, = 56.0.256 molecules 1500
,
I
.
-
140.0
l"ne0Us
I
.
I
.
,
.
I
,
I
E,
accumulatcdc,
130.0
I
I
.
120.0
110.0 100.0 900
80 0 w-
70 0
600 50 0 400
100 -
300
300
1
200
i
I
20.0 0.0
100.0
150.0
00 00
200.0
50.0
Monte Carlo configurations x lod
Figure 6. Accumulated and instantaneous dielectric constants, el,as a function of simulation length for a liquid formamide system of 128 molecules, ~d = m. TABLE 3: Observed Polarization, (P), and Calculated Static Dielectric Constant., cr, for Formamide" NMOL 128 128 256 m 56.0 56.0 6i batch 1 7.121 3.746 5.082 batch 2 6.200 3.915 3.925 batch 3 6.747 5.242 4.716 batch 4 6.289 4.222 4.955 batch 5 batch 6 4.397 6.589 f 0.214 56.2 f 1.8 109f 1
(P) Er
4.413 f 0.239 55.9 f 4.6 109% 1
4.574 f 0.342 59.0 f 6.8 109 f 1
(exptI3 Each batch is 50 million configurations of Monte Carlo sampling.
150.01
,
,
Formamide dielectric constant e,, = 56.0,128 molecules ,
I
,
I
,
I
,
,
,
I
I
130.0 120.0 100 1100 0
t
0.0 0.0 ~
50.0
100.0
150.0
200.0
250.0
300.0
Monte Carlo configurations x lod
Figure 7. Accumulated and instantaneous dielectric constants, er, as a function of simulation length for a liquid formamide system of 128 molecules, ~d = 56.0. ated polarization data are presented in Table 3. A final value for the dielectric constant of formamide of 55.9 f 4.6 was obtained from this simulation. The close agreement of this result with that obtained using conducting boundary conditions suggests that the observed dielectric constant is independent of ed, provided that er 5 ed 5 00. This is in agreement with data
'
'
200
'
"
400
"
600
"
800
"
1000 Monte Carlo configurauons x IOd
'
1200
1400
1600
Figure 8. Accumulated and instantaneous dielectnc constants, el,as a function of simulation length for a liquid formamide system of 256 molecules, = 56.0. reported on other liquids.I3,l4 Furthermore, as the two studies give an identical result to within simulation error, but through markedly different observed polarizations, (P),the results suggest that an applied field of 1.5 x lo8 V m-l is not sufficient to cause dielectric saturation. It is also apparent from a comparison of Figures 6 and 7, together with the statistical uncertainties reported above, that use of a large dielectric constant for the reaction field term results in a larger absolute polarization of the medium and more rapid convergence. This is in direct opposition to the behavior observed using the fluctuation method; there, more rapid convergence is obtained if ed = e r . I 6 , l 7 To determine the sensitivity of the calculated dielectric constant to system size, the calculation was repeated using a box of 256 formamide molecules and a value of 56.0 for E,+ A Monte Carlo simulation covering 150 million configurations was performed on this system, and the polarization and dielectric constant data are presented in Figure 8 and Table 3. A value for the dielectric constant of 59.0 i~ 6.8 was obtained. This agrees, within statistical accuracy, with the values obtained using the smaller system. Further computation to obtain a more precise value for the dielectric constant was not undertaken; however, the results provide no evidence to suggest that use of only 128 molecules is inadequate. Identical conclusions have been reached for the SPC and TIP4P water model^^^.'^ and for a Stockmeyer fluid.17 In Figure 9, the accumulated average dielectric constant for the DMF is presented as a function of simulation length, with an elf value of c-. Also presented are the instantaneous values of Er averaged over the preceding 0.5 million configurations. In Table 4, the values of the observed polarization are presented in batch averages of 50 million configurations, giving a final converged value of 32.1 f 0.3 for er. A second simulation of 200 million configurations was then performed using a value of 32.0 for erf. The polarization and dielectric constant data are given in Table 4 and Figure 10. A final value for the dielectric constant of 31.9 f 1.9 was obtained. The computed er is clearly insensitive to the magnitude of the reaction field dielectric constant, although large values of ed yield more rapid convergence. As observed for formamide, the insensitivity of the results to ed suggests that dielectric saturation is not occumng. Since the static dielectric constant of formamide was shown to be independent of system size, the DMF calculations were
Dielectric Constants of Formamide and Dimethylformamide
J. Phys. Chem., Vol. 99,No. 51, 1995 17961
Dimethylformamide dielectric constant
tion and experiment. Such discrepancies are not unexpected for force fields of this generation since the dielectric constant was not used in the parametrization procedure (see, for example, refs 13, 19, 27, and 30). Not only were the available computational resources insufficient to calculate crreliably, but also the theoretical basis of Er calculations has only become well understood in the past 10 years. One would expect the computed value to be smaller than the experimental one since contributions to the dielectric constant from electronic polarization and anglebond distortion have not been evaluated. However, it is unclear whether these contributions could account for the difference in er of some 50 units.31 Furthermore, the close agreement between various experimental values suggests that the experimental data are reliable.3 For DMF, the average experimental value of cr is approximately 37;3 the computed result with the OPLS unitedatom potential functions of 32 f 2 is in close agreement. It seems likely that the weaker interactions and lack of hydrogen bonding in liquid DMF allow potential functions without explicit polarization to yield a reasonable dielectric constant.
e, = -, 128 molcculcs 50.0
45.0
15.0
c
I
........... instantaneous E,
-accumulated
-
E,
NMOL
128
€4
m
128 32.0
batch 1 batch 2 batch 3 batch 4
3.786 3.705 3.729 3.626
2.612 2.333 2.332 2.726
(P)
3.712 f 0.033 32.1 f 0.3 37f 1
2.501 f 0.100 31.9 f 1.9 37f 1
€r
(exptI3
Each batch is 50 million configurations of Monte Carlo sampling.
Dimethylformamide dielectric constant E"
= 32.0, 128 molecules
80.0
I
75.0
65.0 /
10.0
'
0.0
I
i
I
50.0
1m.o
150.0
200.0
Monte Carlo configurations x 10"
Figure 10. Accumulated and instantaneous dielectric constants, er,as a function of simulation length for a liquid dimethylformamide system of 128 molecules, = 32.0.
not extended to 256 molecules; the value of er for formamide is likely to be more sensitive to simulation protocol owing to cr being larger and the molecular volume being smaller. The experimental value for the dielectric constant of formamide is approximately 109.3.3 The computed result with the OPLS potential-energy parameters is 56 zk 2. Though both values are large, there is clearly disagreement between simula-
Conclusion Monte Carlo simulations of liquid formamide and DMF, as described by the united-atom OPLS potential energy functions, have been performed in both the presence and absence of a reaction field treatment of long-range electrostatic interactions. The thermodynamic properties and radial distribution functions of the liquids are virtually unaffected by changes in system size, simulation ensemble, and the reaction field dielectric constant, E,+ However, angular correlations as represented by the dipoledipole correlation function are very sensitive to simulation conditions. The static dielectric constants of the liquids were calculated by monitoring the system polarization, (P),as a function of an applied electric field of strength, Eo. The consistency of the results as a function of erf and system size suggests that the values of cr are reliable. Furthermore, the insensitivity of the calculated cr to the value of erf supports not only the condition that for polar fluids any value of crf is appropriate provided cr I erf I but also that the applied electric field, Eo, was insufficient to cause dielectric saturation. It is also noticeable that more rapid convergence of (P)is observed for crf= 00, in direct contrast to the result$ observed using the fluctuation method.I6-l7 For formamide, the simulation result for the dielectric constant of 56 f 2 substantially underestimates the experimental value of 109.3. For DMF, however, the simulation result of 32 f 2 is in reasonable agreement with the experimental value of 37, especially given that contributions to er arising from electronic polarization and bondangle deformation of the molecules in the liquid are not evaluated. The poor performance for formamide is perhaps not surprising; force field models with fixed charges tend to underestimate the static dielectric constant for hydrogen-bonded The results suggest that the formamide model should be used with caution in simulations involving ion screening effects but that the DMF model can be used with confidence. Substantial errors are also expected for secondary amides, and by extension, caution is also warranted for the description of dielectric effects in simulations of proteins with such nonpolarizable force fields.
-
Acknowledgment. The authors are grateful for support from the National Institutes of Health and NATOEPSRC and for discussions with Professor K. B. Wiberg. References and Notes (1) Warshel, A.; Aqvist, J. Annu. Rev. Biophys. Biophys. Chem. 1991, 20, 267.
17962 J. Phys. Chem., Vol. 99, No. 51, 1995
Essex and Jorgensen
(2) Jorgensen, W. L.; Swenson, C. J. J . Am. Chem. Soc. 1985, 107, 569. (3) Wohlfarth, C. Landolt-Bornstein: Numerical Data and Functional Relationships in Science and Technology, Group 4, Volume 6, Static Dielectric Constants of Pure Liauids and B i n a n Liauid Mixtures; Springer-Verlag: Berlin, 1991. (4) Smith, P. E.; Brunne, R. M.; Mark, A. E.; van Gunsteren, W. F. J. Phys. Chem. 1993, 97, 2009. (5) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. ( 6 ) de Leeuw, S . W.; Perram, J. W.; Smith, E. R. Annu. Rev. Phys. Chem. 1986, 37, 245. (7) Barker, J. A.; Watts, R. 0. Mol. Phys. 1973, 26, 789. (81 Ewald, P. Ann. Phvs. 1921, 64, 253. (9) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London Ser. A 1980, 373, 27. (10) Onsager, L. J. Am. Chem. Soc. 1936, 58, 1486. (11) Kirkwood, J. G.J. Chem. Phys. 1939, 7, 911. (12) Belhadj, M.; Alper, H. E.; Levy, R. M. Chem. Phys. Lett. 1991, 179, 13. (13) Neumann, M. J. Chem. Phys. 1985, 82, 5663. (14) Zhu, S . B.; Wong, C. F. J. Chem. Phys. 1993, 98, 8892. (15) Neumann, M. Mol. Phys. 1983, 50, 841. I
.
(16) Smith, P. E.; van Gunsteren, W. F. J. Chem. Phys. 1994,100,3169. (17) Kusalik, P. G.J. Chem. Phys. 1990, 93, 3520. (18) Watts, R. 0. Chem. Phys. 1981, 57, 185. (19) Alper, H. E.; Levy, R. M. J. Chem. Phys. 1989, 91, 1242. (20) Alper, H. E.; Levy, R. M. J. Phys. Chem. 1990, 94, 8401. (21) Ni, X.; Fine, R. M. J. Phys. Chem. 1992, 96, 2718. (22) Dang, L. X.; Pettitt, B. M. J. Phys. Chem. 1990, 94, 4303. (23) Kurtovic, Z.; Marchi, M.; Chandler, D. Mol. Phys. 1993, 78, 1155. (24) Kusalik, P. G.; Mandy, M. E.; Svishchev, I. M. J. Chem. Phys. 1994, 100, 7654. (25) Rick, S . W.; Stuart, S. J.; Beme, B. J. J. Chem. Phys. 1994, 101, 6141. (26) Jorgensen, W. L. MCLIQ, Yale University, New Haven, CT, 1994. (27) Neumann, M. J. Chem. Phys. 1986, 85, 1567. (28) Klein, M. L.; McDonald, I. R.; Beme, B. J.; Rao, M.; Beveridge, D. L.; Mehrotra, P. K. J. Chem. Phys. 1979, 71, 3889. (29) Pangali, C.; Rao, M.; Beme, B. J. Mol. Phys. 1980, 40, 661. (30) Fonseca, T.; Ladanyi, B. M. J. Chem. Phys. 1990, 93, 8148. (31) Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1992, 96, 8450. (32) Benson, S . W. Thermochemical Kinetics; Wiley: New York, 1976.
JP95 1486M
