J. Phys. Chem. B 2002, 106, 861-869
861
Solvation of N-Methylformamide by Ethanol: A Comparison of Molecular Dynamics Calculations with the Experimental Data Jan Zielkiewicz* and Jan Mazerski Technical UniVersity of Gdan˜ sk, Department of Chemistry, Narutowicza 11/12, 80-952 Gdan˜ sk, Poland ReceiVed: July 24, 2001; In Final Form: October 23, 2001
The molecular dynamics calculations for {N-methylformamide+ethanol} binary mixture of six compositions and for pure components were performed using various sets of parameters. The results obtained were compared with the thermodynamic data for this system and with deductions derived from experimental data using the Kirkwood-Buff theory of solutions. Analyzing both molecule-molecule and atom-atom radial distribution functions, a good agreement was found between our calculations, and the deductions derived from thermodynamic results. From calculated radial distribution functions the radius of solvation shell was estimated, and values of the local mole fractions were evaluated. Moreover, the formation of complexes was investigated; the lifetime of created NMF-ethanol and NMF-NMF complexes was estimated. The general picture obtained from these calculations is consistent with the thermodynamic results and complements the “thermodynamic” point of view on the solvation of amides by alcohol molecules.
Introduction This work is a part of the systematic study of molecular interactions between the hydroxy and amide groups in the solutions. The amide group can serve as a model of the peptide bond, and the mentioned interactions play an important role in the solvation of the peptides in aqueous solutions. In the previous papers (see refs 1-3, and references given therein) the thermodynamic quantities (the excess Gibbs energy and the molar volumes of mixing) were reported for the binary and ternary mixtures containing amide (N-methylacetamide,1 N,Ndimethylformamide,2 and N-methylformamide3), water, and aliphatic alcohol (methanol, ethanol, and propanol). These quantities were exploited for the description the solvation process of amides by water and/or alcohol: determination of the local mole fractions and orientational effects within the solvation shell. These results were obtained by using the Kirkwood-Buff theory of solutions. The Kirkwood-Buff theory of solutions4 is the exact statistical theory, describing thermodynamic properties such as chemical potential, partial molar volumes, and compressibility factor over the whole concentration range using the quantities
Gij )
∫0∞ (gij(r) - 1)4πr2 dr
(1)
where gij(r) is the radial distribution function and Gij are called the Kirkwood-Buff integrals or fluctuation integrals. This theory links macroscopic (thermodynamic) properties of solution with microscopic description, given by the radial distribution functions, gij(r), and the link is provided by eq 1. In other words, if radial distribution functions, gij, are known (and then the Gij quantities are known, too), we can determine the thermodynamic properties of solution using the equations of the KirkwoodBuff theory.4 On the other hand, as it was pointed first by Ben Naim,5 the original Kirkwood-Buff theory can be inversed, and the Gij parameters can be determined from the experimental * Corresponding author: e-mail:
[email protected].
values of thermodynamic quantities (chemical potential, partial molar volumes, and isothermal compressibility factor). Next, the ciGij product (where ci ) Ni/V is the concentration of species i in the mixture) describes the total average excess (or deficiency) of i molecules in the entire neighborhood of a j molecule.5,6 Then, exploiting the idea presented previously in the literature,1-3,5,6 it is possible to estimate the local mole fractions of the components of the solution. The possibility of evaluating these quantities was pointed out by many authors6-10 and seems a valuable and convenient tool to describe the solvation process. As it was stated above, the Kirkwood-Buff theory is an exact one, and therefore the determination of the Gij values from thermodynamic data does not involve any assumptions. That is why the description of the solvation phenomena obtained on the basis of this approach is fully reliable, and the Kirkwood-Buff integrals formalism is a widely applied method in the solution investigations.11,12 It is quite obvious that the pure thermodynamic data cannot give a detailed insight into the structure of the solvation shell. The information on the local structure of this shell is contained in the radial distribution functions, but these details are lost in the integration process (eq 1). On the other hand, the molecular dynamics calculations may provide a complete information on the solution structure (the radial distribution functions: moleculemolecule and atom-atom, formation of complexes, and orientational effects). We can say, that the results of MD calculations may complement and confirm the conclusions obtained from the thermodynamic measurements. This is why this work has been initiated. For correct reproduction of the thermodynamic properties of the mixture, a proper and reliable set of intermolecular interactions parameters is required. Many papers report the molecular dynamics or Monte Carlo simulations, but a relatively small number of them concern simulations of mixtures, especially over the whole composition range. Although the pure component properties are described correctly, it is not quite obvious that thermodynamic properties of mixture will be correctly described,
10.1021/jp012833x CCC: $22.00 © 2002 American Chemical Society Published on Web 12/28/2001
862 J. Phys. Chem. B, Vol. 106, No. 4, 2002
Zielkiewicz and Mazerski and B molecules is given in the form
TABLE 1: Geometrical Parameters of NMF and Ethanol Molecules bond lengths (nm)
bond angles (deg)
EAB )
O-H CH2-O CH2-CH3
Ethanol 0.100 C1-O-H 0.143 C2-C1-O 0.153
109.5 111
CAdO CA-H1 CA-N NA-H2 N-CH3
N-Methylformamide 0.123 H1-CA-O 0.109 H1-CA-N 0.134 CA-N-H2 0.100 CA-N-CH3 0.147
120 120 116 123
∑i ∑j
(
qiqje2
+
rij2
(C12)ij
-
)
(C6)ij
rij12
rij6
(2)
Standard combining rules were used:
(C6)ij ) {(C6)ii(C6)jj}1/2 (C12)ij ) {(C12)ii(C12)jj}1/2
(3)
Two sets of parameters for NMF and three sets for ethanol molecule were used. For the NMF(1) molecule the LennardJones parameters were taken from the GROMOS 96 force field, and the partial charges were taken from Jorgensen and Gao.14 For the NMF(2) molecule all the parameters were taken from Jorgensen and Gao14 (the OPLS parameters). For the ethanol(1) molecule all the parameters were taken from the GROMOS 96 force field. For the ethanol(2) molecule the Lennard-Jones parameters were taken from the GROMOS 96 force field, and partial charges on atoms were taken from Jorgensen.15 For the ethanol(3) molecule all the parameters were taken from Jorgensen15,16 (OPLS parameters). Geometrical parameters of NMF and ethanol molecules are given in Table 1 and values of the parameters that were used are summarized in Table 2. To mimic correctly the hydrogen bond (using the GROMOS 96 force field), GROMOS 96 increases the Lennard-Jones repulsion term, C12, for the polar atoms if they interact with the hydrogen. That is why two values of C12 parameters are given in Table 2 for the oxygen and nitrogen atoms. The initial configurations of the investigated mixtures were created by arranging the appropriate number of NMF and ethanol molecules in a cubic box. Then, the system was equilibrated by at least 1500 ps or more with the time step of a single iteration equal to 1 fs. Finally, the 600-750 ps run was performed for data collection. The atomic coordinates were
too. Then, we use various sets of parameters in our simulations, and results of this calculations are compared with experiment. The aim of this work is a direct comparison of the results of the thermodynamic measurements and conclusions derived from them with the results obtained from molecular dynamics calculations. The investigated system is the binary mixture {Nmethylform-amide(NMF)+ethanol}, at the temperature 313.15 K. Calculations The molecular dynamics calculations were carried out by using the GROMOS 96 package.13 Each model system consisted of 729 molecules in a cubic box with the periodic boundary conditions. The calculations were carried out at the NPT ensemble, assuming a temperature of T ) 313 K and a pressure p equal to 1 atm. The molecules were represented as a collection of interacting sites: a four-site set for ethanol and a six-site (or five-site for OPLS model) set for N-methylformamide. The CH3 and CH2 groups were taken as a single sites centered on the carbon atom (the united atom approximation). There are static, nonpolarizable models whose site-site interactions are sums of the Coulombic potentials and the short-range Lennard-Jones term. The energy of intermolecular interactions between the A
TABLE 2: Parameters of eq 2 for the ith site of the amide and ethanol molecules (C6 and C12 given in kJ mol-1 nm6 and kJ mol-1 nm12 units) N-Methylformamide cis C6 × 103
atom
trans
C12 × 106
CA N
2.341 2.436
O
2.262
H1 H2 CH3
8.464‚10-2 0.0 5.513
3.374 1.693 3.775 0.7415 1.266 0.0151 0.0 15.116
CA N O H2 CH3
5.791 3.352 2.362 0.0 8.556
1.994 3.952 1.590 0.0 25.796
C6 × 103
q
C12 × 106
q 0.58 -0.55
8.464‚10-2 0.0 5.513
3.374 1.693 3.775 0.7415 1.266 0.0151 0.0 15.116
5.791 3.352 2.362 0.0 8.556
1.994 3.952 1.590 0.0 25.796
0.58 -0.55 -0.53 0.30 0.20
NMF(1) 0.53 -0.55
2.341 2.436
-0.53
2.262
0.0 0.35 0.20 NMF(2) 0.53 -0.55 -0.53 0.35 0.20
-0.53 0.0 0.30 0.20
Ethanol atom
C6 × 103
H O
0.0 2.262
CH2 CH3
7.105 9.916
C12 × 106 Ethanol(1) 0.0 1.266 1.505 25.776 33.547
q
C6 × 103
0.398 -0.548
0.0 2.262
0.150 0.0
7.105 9.916
C12 × 106 Ethanol(2) 0.0 1.266 1.505 25.776 33.547
q
C6 × 103
C12 × 106
q
0.435 -0.700
0.0 2.381
Ethanol(3) 0.0 1.994
0.435 -0.700
0.265 0.0
7.005 10.000
24.840 29.009
0.265 0.0
Solvation of N-Methylformamide
J. Phys. Chem. B, Vol. 106, No. 4, 2002 863
TABLE 3: Details of Simulations No. molecules NMF ethanol
simulation time (ps)
box length l (nm)
box volume V (cm3/mol)
750 750 750 750 600 750 810
4.192 ( 0.006 4.175 ( 0.006 4.154 ( 0.005 4.143 ( 0.005 4.132 ( 0.005 4.124 ( 0.005 4.106 ( 0.005
60.87 ( 0.24 60.15 ( 0.24 59.25 ( 0.23 58.77 ( 0.21 58.30 ( 0.22 57.94 ( 0.20 57.18 ( 0.22
cis-NMF(1) + ethanol(1) -57.01 ( 0.23 -75.63 ( 0.24 -99.27 ( 0.26 -113.26 ( 0.25 -127.42 ( 0.27 -136.50 ( 0.27 -158.19 ( 0.34
135 243 378 459 540 594 729
594 486 351 270 189 135 0
0.185 0.333 0.518 0.630 0.741 0.815 1.000
750 600 750 750 600 600 510
4.195 ( 0.006 4.180 ( 0.005 4.162 ( 0.005 4.149 ( 0.006 4.139 ( 0.006 4.133 ( 0.005 4.117 ( 0.005
trans-NMF(1) + ethanol(1) 60.99 ( 0.29 -53.04 ( 0.25 60.64 ( 0.21 -67.08 ( 0.24 59.57 ( 0.24 -88.15 ( 0.24 59.02 ( 0.24 -96.58 ( 0.25 58.61 ( 0.24 -107.46 ( 0.25 58.33 ( 0.20 -114.98 ( 0.24 57.64 ( 0.20 -134.24 ( 0.25
135 243 378 459 540 594 729
594 486 351 270 189 135 0
0.185 0.333 0.518 0.630 0.741 0.815 1.000
600
4.126 ( 0.006
61.91 ( 0.26
Ethanol(1) -34.14 ( 0.25
0
729
0.000
750 750 750 750 600 750 810
4.139 (0.006 4.127 (0.005 4.117 (0.005 4.114 (0.005 4.113 (0.005 4.110 (0.005 4.106 (0.005
58.58 (0.25 58.06 (0.23 57.67 (0.21 57.53 (0.23 57.49 (0.21 57.35 (0.22 57.18 (0.22
cis-NMF(1) + ethanol(2) -68.90 (0.24 -85.03 (0.24 -105.51 (0.25 -117.79 (0.26 -130.31 (0.26 -138.49 (0.31 -158.19 (0.34
135 243 378 459 540 594 729
594 486 351 270 189 135 0
0.185 0.333 0.518 0.630 0.741 0.815 1.000
750 600 600 750 600 600 510
4.141 (0.006 4.136 (0.007 4.126 (0.006 4.128 (0.005 4.125 (0.005 4.120 (0.005 4.117 (0.005
58.68 (0.26 58.27 (0.27 57.92 (0.23 57.70 (0.36 57.91 (0.21 57.83 (0.21 57.64 (0.20
trans-NMF(1) + ethanol(2) -64.68 (0.21 -76.75 (0.43 -94.35 (0.30 -101.78 (0.70 -110.28 (0.22 -117.03 (0.25 -134.24 (0.25
135 243 378 459 540 594 729
594 486 351 270 189 135 0
0.185 0.333 0.518 0.630 0.741 0.815 1.000
600
4.155 (0.006
59.27 (0.26
Ethanol(2) -49.32 (0.20
0
729
0.000
100
Large system: trans-NMF + ethanol(2), RCUTP)2.0 nm, RCUTL)3.5 nm 8.248 (0.004 57.94 (0.08 -98.01 (0.08 3024
2808
0.518
total potential energy Epot (kJ/mol)
mole fraction of amide xA
750 750 750 750 600 750 750
4.152 (0.006 4.164 (0.005 4.179 (0.005 4.189 (0.005 4.202 (0.005 4.208 (0.005 4.236 (0.006
59.15 (0.25 59.64 (0.21 60.30 (0.21 60.75 (0.23 61.31 (0.20 61.55 (0.22 62.81 (0.25
cis-NMF(2) + ethanol(3) -64.28 (0.23 -80.88 (0.25 -101.93 (0.25 -114.57 (0.26 -127.38 (0.25 -135.93 (0.26 -157.28 (0.25
750 750 750 750 750 750 750
4.160 (0.005 4.174 (0.005 4.194 (0.005 4.204 (0.005 4.214 (0.005 4.223 (0.005 4.239 (0.005
59.50 (0.22 60.11 (0.22 60.93 (0.24 61.39 (0.23 61.83 (0.24 62.23 (0.23 62.95 (0.23
trans-NMF(2) + ethanol(3) -58.92 (0.21 -71.11 (0.22 -86.48 (0.23 -95.84 (0.23 -105.27 (0.23 -111.53 (0.23 -127.60 (0.23
135 243 378 459 540 594 729
594 486 351 270 189 135 0
0.185 0.333 0.518 0.630 0.741 0.815 1.000
600
4.140 (0.006
58.60 (0.25
Ethanol(3) -44.02 (0.20
0
729
0.000
collected in the time interval equal to 2-3 ps for calculations of various quantities. The large system, obtained by multiplication of equilibrated “small” system, contains 729 molecules. The obtained structure was subsequent equilibrated by 250 ps, and finally the 100 ps run for data collection was performed. Details of our simulations are given in Table 3. The GROMOS 96 molecular dynamics program offers the twin-range method for evaluating the long-range nonbonding interactions.13 The short-range cutoff distance, RCUTP, was assumed equal to 1.2 nm, and the long-range cutoff distance,
135 243 378 459 540 594 729
594 486 351 270 189 135 0
0.185 0.333 0.518 0.630 0.741 0.815 1.000
RCUTL, was equal to 1.9 nm. The list of interacting pairs at distances between RCUTP and RCUTL was updated at every fifth simulation step. For the large system these values are RCUTP ) 2.0 nm and RCUTL ) 3.5 nm. In this case the calculations are very time expensive: the pair list contains 2.3 million elements approximately. Results and Discussion The NMF molecule can exist in the solution in two forms: the cis and trans conformers (see Figure 1). The barrier to
864 J. Phys. Chem. B, Vol. 106, No. 4, 2002
Zielkiewicz and Mazerski
Figure 1. The cis and trans conformers of NMF molecule. For the OPLS model the H-C group is approximated by one site (the united atom approximation).
internal rotation in monosubstituted amides (N-methylformamide and N-methylacetamide) was determined experimentally17 and is equal to 22.6 kcal/mol for trans-to-cis conversion, and 19.8 kcal/mol for cis-to-trans conversion. These values were determined from the NMR measurements, using 1,2-dichloroethane and water as solvents. It is easy to calculate that the trans conformer is 2.8 kcal/mol more stable. The error of these measurements was estimated17 to be 1.8 kcal/mol, and therefore the measurement is no possible evaluation of the solvent effect on the relative stability of both conformers. On the other hand, we can expect that these conformers may be different in their ability for hydrogen bond formation. Dixon et al.18 have found that for N-methylacetamide, (NMA), the cis isomer is predicted to be 2.3 kcal/mol less stable than the trans one for the uncomplexed NMA. In the water solution, the cis-NMA‚H2O complex is only 0.5 kcal/mol less stable than the trans-NMA‚ H2O complex. Taking into investigation the dimeric molecule, (NMA)2, formed by two amide molecules, they found that the cis isomer forms a dimeric complex 6.7 kcal/mol more stable than the complex formed by the trans isomer. This indicates that the cis and trans conformers interact with water or other amide molecules in a different manner. We can expect that in the alcohol solution, a similar behavior for NMF should be observed. Thus, for the binary mixture {NMF+EtOH} investigated in this work, the calculations for the cis and trans isomers of a NMF molecule were performed separately. In the discussion, we take into account the following quantities: the total potential energy, Epot, the molar volume, Vmol, the molecule-molecule and atom-atom radial distribution functions, gij, and formation of complexes between amide and alcohol molecules. (a) Total Potential Energy and Molar Volumes. The total potential energy may be split into three terms: the LennardJones interactions, the Coulomb (electrostatic) interactions, and the potential energy of distortions of molecules (bond lengths, bond angles, etc.). For the ideal mixture, the total potential energy plotted against the mole fraction should be a straight line. For real solutions the deviations from this line will be observed as heat of mixing. From the experimental data for heat of mixing19 and for excess volumes of mixing20 the for {NMF+EtOH} binary mixture, it appears that both the total potential energy and the box volume should be varied linearly versus the amide mole fraction. Observed deviations from this linear dependence are very small: at the equimolar composition the heat of mixing, HE ) +0.6 kJ mol-1, and excess volume of mixing, VE ) -0.09 cm3 mol-1. Results of our calculations of Epot and Vmol are collected in Table 3; additionally, they are presented in the Figure 2. As may be seen from this figure, the linear variation of the total potential energy, Epot, against the mole fraction of NMF is roughly fulfilled for both cis and trans conformers. Using the OPLS partial charges, there are visible differences between the two conformers: for example, the trans form deviates from the
Figure 2. Composition dependence of the total potential energy of intermolecular interactions, Epot, and the molar volume of mixture, Vmol. The dotted line represents the experimental data.
straight line. Observed deviations correspond with the endothermal mixing, but compared with experiment their absolute values are too large. For the molar volume there are no visible differences between the conformers (omitting difference in the molar volume of pure forms). For both forms, negative deviations of Vmol from the ideal linear dependence are observed. The sign of these changes agrees with experiment, but their absolute values are too large. Notice that combination of the Lennard-Jones parameters from the GROMOS 96 force field and partial charges from Jorgensen15 gives the best reproduction of the molar volumes of pure compnents. (b) Molecule-Molecule Radial Distribution Functions. For the investigated system {N-methylformamide+ethanol}, the molecule-molecule radial distribution functions (between centers of mass), rdf, g11(r), g12(r), and g22(r) are shown in Figure 3 for all the compositions. The graphs represent results for the NMF(1) and ethanol(2) set of parameters; for other sets of parameters the picture is very similar. As can be seen from this figure, for g11(r) and g12(r) rdf’s (describing solvation of the NMF molecule by NMF and ethanol, respectively) the welldefined first peak of the rdf is observed; this peak corresponds with the first coordination sphere around the amide molecule. We define here the solvation shell as a sphere of radius Rc, centered on the amide molecule, and within this sphere the mole fractions of components differ from the bulk ones.6-10 The radius of the first coordination sphere can be estimated from the position of the first minimum on the rdf; it equals to 0.65 nm approximately. Moreover, the second peak of the rdf is visible too, and it corresponds with the second coordination sphere. Additionally, it can be concluded that moleculemolecule rdf’s for both the cis and trans conformers differ only slightly. Variation of the number of ethanol and NMF molecules, along with the total number of molecules within the solvation shell, is presented in the Figure 4. The total number of molecules within this shell (at Rc < 0.65 nm) is nearly constant, equal to 38.3 approximately for the NMF(1)+ethanol(1) system, 38.8
Solvation of N-Methylformamide
J. Phys. Chem. B, Vol. 106, No. 4, 2002 865
Figure 5. Values of the local mole fractions, xijloc(r), as a function of radius r, calculated from eqs 6 and 7. Figure 3. Molecule-molecule radial distribution functions, g11(r), g22(r), and g12(r), calculated at various compositions for NMF(1) and ethanol(2) sets of parameters.
can be read from Figure 3. As it was previously found,3 the linear coefficient of preferential solvation describing solvation of the NMF molecule by ethanol is δ21 ) -5 cm3/mol at equimolar composition (see Figure 1 in ref 3). Then, the local mole fraction of ethanol around the NMF molecule can be estimated as
x21 ) x2 +
Figure 4. Concentration dependence of the number of molecules within the solvation shell around the amide molecule, calculated for the NMF(2) and ethanol(3) sets of parameters: 1: ethanol; 2: NMF; 3: total.
for the NMF(1)+ethanol(2) system, and 36.0 for the NMF(2)+ethanol(3) system. It is the same for both conformers, and, as can be seen, is nearly independent of the choice of parameter set. It was previously deduced3 from thermodynamic results that for the investigated binary mixture {NMF+ethanol} the local mole fractions of components are nearly the same as the bulk ones. Moreover, it was stated3 that the radius of solvation shell, defined as above, around the NMF molecule does not exceed the radius of the first coordination sphere. We try to examine these conclusions as follows. Using the radius of solvation shell, Rc, it is easy to calculate its volume:
4 Vsolv ) π(0.653 - 0.303) ) 1.037 nm3 ) 624 cm3/mol (4) 3 where 0.30 nm is the “hard core” radius of amide molecule; it
δ21 5 ) 0.5 ) 0.492 Vsolv 624
(5)
In other words, it is smaller than the bulk value by 2% approximately. On the other hand, we can calculate the local mole fractions from our molecular dynamics results. According to the definition of the molecular radial distribution function, gij(r), we can write
∫0r gij(x)4πx2 dx
nij(r) ) Fi
(6)
where nij(r) symbolizes the number of i molecules within the sphere of radius r around the central j molecule and Fi is the number density of i molecules (Fi ) Ni/V). Then, the local mole fractions within this sphere can be calculated as follows:
xijloc(r) )
nij(r)
∑k nkj(r)
(7)
The xijloc(r) functions calculated from this equation are displayed in Figure 5. The solvation of the amide molecule describes x11 (amide around amide) and x21 (ethanol around amide) local mole fractions. As can be seen, the local mole fractions in the limit of the first coordination sphere (at Rc ) 0.65 nm) differ from the bulk ones. Especially, the local mole fraction of ethanol around the NMF molecule is smaller than the bulk one. The difference, however, exceeds the estimated above 2% at
866 J. Phys. Chem. B, Vol. 106, No. 4, 2002
Zielkiewicz and Mazerski
Figure 6. Kirkwood-Buff integrals, Gij(r), as a function of radius r, calculated from molecule-molecule radial distribution functions and using eq 8, performed for the mixture at the amide mole fraction equal to 0.518 and using the NMF(1) and ethanol(2) sets of parameters.
equimolar composition. Then, we can conclude that our results agree with experimental data only qualitatively. As may be seen from Figure 5, in the limit of the second solvation shell (at r > 1 nm) the local mole fractions do not differ from the bulk ones. That behavior was found for aqueous amide solutions, too; the Marcus thermodynamic investigations21 show that correlation in the second solvation shell is very weak and disappears altogether from the third shell onward. Similar conclusions were derived for the (formamide+water) mixture and using molecular dynamics calculations. Puhovski and Rode22 have found that the local mole fractions of components around the formamide molecule are nearly equal to the bulk ones. This is for water solutions; we found previously, nevertheless, that solvation of amides by water and alcohols is very similar.1-3 As it was described in the introduction, the Kirkwood-Buff integrals can be calculated from GE and VE data given in the literature.23,20 For example, at xA ) 0.518 these values are G11 ) -45.6 cm3 mol-1, G12 ) -69.6 cm3 mol-1, and G22 ) -43.5 cm3 mol-1. Then we try to interpret the molecule-molecule radial distribution functions in terms of the Kirkwood-Buff integrals (defined by eq 1). In our calculations we have chosen a relatively large number of molecules (729 molecules), because we expected that the Gij functions
Gij(r) )
∫0r (gij(x) - 1)4πx2 dx
(8)
should be converged asymptotically to the finite value if the radius r increases infinitely. In this limit, the Gij value should be equal to the experimental one. Such asymptotic behavior was observed by Kojima et al.24 for the calculated Gij(r), obtained by solving the Percus-Yevick equation (see Figure 1 in ref 24). However, for the investigated system {NMF+EtOH}, the Gij(r) functions do not converge. Figure 6 shows the results of these calculations at exemplary composition xA ) 0.518. For other compositions the results are very similar. It is easy to
Figure 7. Kirkwood-Buff integrals, Gij(r), as a function of radius r, calculated for the large system at the amide mole fraction equal to 0.518 and using the NMF(1) and ethanol(2) sets of parameters.
understand this result: since the [gij(x) - 1] function (eq 8) is multiplied by the factor x2 prior to the integration, even small deviations of the gij(x) from the unit value (visible on the graphs in Figure 3) give, at large x, large contributions to the Gij value. As a consequence, the divergence of Gij(r) is observed. The results obtained indicate, therefore, that even 729 molecules are insufficient to calculate the Kirkwood-Buff integrals by integration of the molecular rdf. As we suppose, similar difficulties were probably observed by Chitra and Smith25 for the binary {2,2,2-trifluoro-ethanol+water} mixture. The authors have calculated the Kirkwood-Buff integrals using the upper limit of integration equal to 1 nm; the calculated values differ, however, from experiment even by factor 10 or more. Since in the cited paper25 the functions Gij(r) are not displayed, then direct comparison with our results is impossible. We may suppose that for large systems the approximation of the Gij values should be better. Then, we perform the simulation run for the large system containing 3024 t-NMF molecules and 2808 ethanol molecules (5832 molecules in total, xA ) 0.518). The set of parameters used in the calculations is: NMF(1) and ethanol(2). Total simulation time was equal to 100 ps, coordinates were collected after each 2 ps. Results of integration are visible in Figure 7. As may be seen, the integration procedure in this case fails too. The possible explanations of this fact are that the investigated system is still too small or simulation time is too short. (c) Atom-Atom Radial Distribution Functions. We take into account the following ethanol-amide atom-atom radial
Solvation of N-Methylformamide
J. Phys. Chem. B, Vol. 106, No. 4, 2002 867
Figure 8. Atom-atom radial distribution functions, calculated for the NMF(1) and ethanol(2) sets of parameters.
distribution functions: Hethanol-Oamide, Oethanol-Oamide, HethanolNamide, Oethanol-Namide, and HN-amide-Oethanol radial distribution. They are displayed in Figure 8 for the NMF(1) and ethanol(2) sets of parameters, and at composition xA ) 0.518. For other sets of parameters and other compositions the picture is very similar. As can be seen, the strongly orientational effects of ethanol molecules around the oxygen atom of amide are observed; as we can expect, the orientation of ethanol molecules around the amide nitrogen is less visible. These effects are clearly visible within the solvation shell; within the second coordination sphere they are nearly invisible (Figure 8). Similar orientational effects of water around the formamide molecule, resulting from hydrogen bond creation, were observed by Puhovski and Rode.22 They found that there are significant orientational effects for both water-water and water-formamide interactions, observed at short ranges corresponding to the first solvent shell. The most intense peaks indicate the formation of the strong hydrogen bonds Hethanol-Oamide, Hethanol-Oethanol, and HN-amideOethanol. Integration of these peaks allows to estimate the mean number of hydrogen bonds created by alcohol or amide molecules. Especially interesting is the estimation of the alcohol structure in the solution. Experimental results (small deviations from the Raoult’s law, small deviations of the local mole fractions from the bulk ones, and small values of the heat of mixing) indicate that ethanol retains its structure over the whole composition range. The simple quantitative measure of this structure in the mixture is the total number of hydrogen bonds created by the alcohol molecule. Results of these calculations are presented in the Figure 9. In this figure the differences between various sets of parameters are clearly visible. In our opinion, the GROMOS 96 parameters for ethanol (ethanol(1)) lead to the disagreement with the experiment. Using this set of parameters, the mean number of hydrogen bonds created by the ethanol molecule strongly decreases if the amide mole fraction increases. Only for the OPLS partial charges does ethanol retain its structure over the whole composition range. Figure 10 presents the mean number of hydrogen bonds created by the oxygen atom of the amide molecule. As can be seen, the oxygen atom forms 1.7-1.8 hydrogen bonds approximately, and this value nearly does not depend on the mole fraction of amide. The OPLS parameters were successfully used for simulations of pure amides26,27 and ethanol.28-31 It can be concluded now that OPLS partial charges correctly describe the binary mixture investigated in this work. (d) Formation and Mean Lifetime of the HydrogenBonded Complexes. The clearly visible first peak of the ethanol-amide Hethanol-Oamide and HN-amide-Oethanol atomatom radial distribution functions indicates that a hydrogen bond between NMF and ethanol molecules exists. This hydrogen bond is responsible for formation of complexes between these molecules. Creation of these complexes was reported previously
Figure 9. Mean number of hydrogen bonds created by one ethanol molecule. 1: ethanol-NMF; 2: ethanol-ethanol; 3: total. The cutoff radius Rc ) 0.25 nm.
Figure 10. Mean number of hydrogen bonds created by one NMF molecule. 1: NMF-ethanol; 2: NMF-NMF; 3: total. The cutoff radius Rc ) 0.275 nm.
in the literature.19,32-34 Thus, we may consider the complexes appearing in the resulting files obtained from our molecular dynamics simulations. We assume that a complex is formed if the hydrogen bond between molecules is created. The hydrogen bond in the molecular dynamics calculations can be defined using either energetic15 or geometric29,30 criteria. In our calculations we use for simplicity only the cutoff RH for the donor
868 J. Phys. Chem. B, Vol. 106, No. 4, 2002
Zielkiewicz and Mazerski TABLE 4: Fitted N0 and τ Values (eq 9) for the Investigated Systems Hamide-Oethanol complexes, binary mixture, xA ) 0.518 NMF(1)+ethanol(1) NMF(1)+ethanol(2) NMF(2)+ethanol(3) N0 τ/ps
cis
trans
cis
trans
cis
trans
95 2.86
90 2.04
107 7.79
129 5.74
104 5.24
136 4.63
Hamide-Oamide complexes, binary mixture, xA ) 0.518 NMF(1)+ethanol(1) NMF(1)+ethanol(2) NMF(2)+ethanol(3) N0 τ/ps
cis
trans
cis
trans
cis
trans
110 5.57
83 3.01
112 7.40
101 4.67
108 6.08
92 4.45
Hamide-Oamide complexes, pure amide NMF(1) N0 τ/ps
Figure 11. Mean number of complexes formed by the Hethanol and Oamide atoms as a function of the elapsed time. Points represent results of calculations, and the solid line is calculated from eq 9, using the fitted N0 and τ values.
oxygen-acceptor hydrogen distance. According to this criterion, the complex is formed if the distance H-O does not exceed some limiting value RH. We take into account both HethanolOamide and HN-amide-Oethanol hydrogen bonds, and assume RH ) 0.25 nm and RH ) 0.27 nm, respectively. The calculations were carried out as follows. The trajectory of the system in the configuration space was written at the time interval ∆t (usually 3 ps). First, we look the complexes in the starting coordinate file, written at the moment t0, for which we find, say, N(t0) complexes. Next, we look how many of these N(t0) complexes exist in the second file, written at the time t0+∆t, and we find, say, N(t0 + ∆t) complexes. Then, we look at how many of these N(t0 + ∆t) complexes exist in the subsequent file written at t0 + 2∆t, and so on. Finally, this procedure was averaged over the whole simulation run (usually 750 ps). The mean lifetime of complexes, τ, can be calculated from the equation
N ) N0e-t/τ
(9)
where N symbolizes number of complexes observed at the time t. In this equation N0 and τ are adjustable parameters. Results of our calculations for complexes formed by the Hethanol and Oamide atoms are given in Figure 11 at the composition xA ) 0.518 and for various sets of parameters; fitted N0 and τ values are given on the graphs. As may be seen from this figure, eq 9 quite correctly describes our results. Differences between various sets of parameters are clearly visible. Moreover, differences between both conformers are clearly visible for ethanol(2) and ethanol(3) sets of parameters: the trans conformer always creates more stable complexes than does the cis conformer. It should be stated here that calculated lifetime, τ, strongly depends on the set of parameters chosen. The same calculations were carried out for complexes formed by Hamide and Oethanol atoms too. Results of this investigation
NMF(2)
cis
trans
cis
trans
345 6.85
338 3.73
324 4.87
328 3.77
are reported in Table 4. In this case the cis form creates more stable complexes (the lifetime is longer). It may be explained as an effect of formation the “cyclic” complex using two hydrogen bonds: Hamide-Oethanol and Hethanol-Oamide. Results of our calculations can be compared with literature data. Recently Woutersen et al.38 have found, using 2D-IR spectroscopy, that the lifetime of the hydrogen bond formed by the amide oxygen and methanol is equal to 10-15 ps. It may be supposed that for ethanol complexes this lifetime should be longer; this may be concluded from the NMR relaxation data for alcohols.36 Then, we can say that the lifetime value obtained by us for NMF-ethanol complexes agrees with the available literature data. The procedure described above can be adopted to estimate the lifetime for the Hamide-Oamide hydrogen bond in the amideamide complex. These calculations were carried out in the pure amide (xA ) 1.0) and in the binary mixture at composition xA ) 0.518. The obtained results are included in the Table 4, too. As may be seen, the cis conformer forms more stable complexes (the lifetime is longer). This can be explained as effect of formation of the cyclic dimers by the cis form of amide (see Figure 1). From the NMR data the correlation time, τc, can be determined; this quantity is coupled to the lifetime of the hydrogen bonded state.36,37 In general, the correlation time for H-bonded molecules is longer than that for “homologic” molecules without hydrogen bonds, and the lifetime of the hydrogen bond is longer than the correlation time τc.36 Unfortunately, there is no direct link between the correlation time τc and the lifetime of the bonded state. Recently Seipelt and Zeidler have found that correlation time in the pure (liquid) NMF is equal to 4.7-4.9 ps.38 It is impossible, nevertheless, to estimate the lifetime of the hydrogen bond from these results. It may be concluded only that the lifetime values obtained by us are not in conflict with the available data. Conclusions 1. Comparison of Various Sets of Parameters. In our opinion, the best description is given by the combination of partial charges from OPLS model and the Lennard-Jones parameters from the GROMOS 96 force field. This combination gives the best reproduction of the molar volumes of pure liquids. The full OPLS set of parameters (partial charges and Cij
Solvation of N-Methylformamide coefficients) gives, generally, very similar results. The partial charges for ethanol taken from the GROMOS 96 force field give, in our opinion, unreliable concentration dependence of mean number of hydrogen bonds per alcohol molecule. 2. Analysis of Various Quantities. By analyzing various quantities (density of mixture, total potential energy, values of the local mole fraction, hydrogen bond, and creation of complexes) it may be concluded that the general picture obtained from analysis of the molecular dynamics results is consistent with thermodynamic data. The agreement between our calculations and experiment is, however, qualitative only. Moreover, it was impossible to calculate very important quantities, the Kirkwood-Buff integrals, even for the large system. Our results complement the experimental data: for example, the differences between both conformers are visible and the lifetime of hydrogen bonded complexes was estimated. Acknowledgment. Calculations were carried out at the Academic Computer Center in Gdan˜sk. References and Notes (1) Zielkiewicz, J. Phys. Chem. Chem. Phys. 2000, 2, 2925-2932. (2) Zielkiewicz, J. J. Phys. Chem. 1995, 99, 4787-4793. (3) Zielkiewicz, J. J. Chem. Soc., Faraday Trans. 1998, 94, 17131719. (4) Kirkwood, J. G.; Buff F. P. J. Chem. Phys. 1951, 19, 774-777. (5) Ben Naim, A. J. Chem. Phys. 1977, 67, 4884-4890. (6) Ben Naim, A. J. Chem. Phys. 1989, 93, 3809-3813. (7) Ben Naim, A. Pure Appl. Chem. 1990, 62, 25-34. (8) Mansoori, A. G.; Ely J. F. Fluid Phase Equilib. 1985, 22, 253275. (9) Rubio, R. G.; Prolongo, U.; Cabrerizo, M.; Diaz Pena, M.; Renuncio J. A. R. Fluid Phase Equilib. 1986, 26, 1-13. (10) Rubio, R. G.; Prolongo, U.; Diaz Pena, M.; Renuncio J. A. R. J. Phys. Chem. 1987, 91, 1177-1184. (11) Fluctuation Theory of Mixtures; Matteoli, E., Mansoori, G. A., Eds.; Taylor and Francis: New York, 1990. (12) Newman, K. E. Chem. Soc. ReV. 1994, 23, 31-40.
J. Phys. Chem. B, Vol. 106, No. 4, 2002 869 (13) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R.; Tironi, L. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide, Biomos b.v., Zu¨rich, Groningen, 1996. (14) Jorgensen, W. L.; Gao, J. J. Am. Chem. Soc. 1988, 111, 42124216. (15) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 1276-1284. (16) Jorgensen, W. L. J. Am. Chem. Soc. 1984, 106, 6638-6646. (17) Drakenberg, T.; Frosen, S. J. Chem. Soc., Chem. Commun. 1971, 1404-1405. (18) Dixon, D. A.; Dobbs, K. D.; Valentini, J. J. J. Phys. Chem. 1994, 98, 13435-13439. (19) Pikkarainen, L. Thermochim. Acta 1991, 178, 311-319. (20) Zielkiewicz, J. J. Chem. Thermodyn. 1996, 28, 313-318. (21) Marcus, Y. J. Chem. Soc., Faraday Trans. 1990, 86, 2215. (22) Puhovski, Y. P.; Rode, B. M. J. Phys. Chem. 1995, 99, 15661576. (23) Zielkiewicz, J. J. Chem. Thermodyn. 1997, 29, 229-237. (24) Kojima, K.; Kato, T.; Nomura, H. J. Solution Chem. 1984, 13, 151165. (25) Chitra, R.; Smith, P. E. J. Chem. Phys. 2001, 114, 426-435. (26) Essex, J. W.; Jorgensen, W. L. J. Phys. Chem. 1995, 99, 1795617962. (27) Yashonath, S.; Rao, C. N. R. Chem. Phys. 1991, 155, 351-356. (28) van Leeuwen, M. E. Mol. Phys. 1998, 87, 87-101. (29) Saiz, L.; Padro, J. A.; Guardia, E. J. Phys. Chem. B 1997, 101, 78-86. (30) Saiz, L.; Guardia, E.; Padro, J. A. J. Chem. Phys. 2000, 113, 28142822. (31) Fidler, J.; Rodger, P. M. J. Phys. Chem. 1999, 103, 7695-7703. (32) Schmid, E. D.; Brodbek, E. Can. J. Chem. 1985, 63, 1365-1371. (33) Eaton, G.; Symons, M. C. R.; Rastogi, P. P. J. Chem. Soc., Faraday Trans. 1989, 85, 3257-3271. (34) Eberhardt, E. S.; Raines, R. T. J. Am. Chem. Soc. 1994, 116, 21492150. (35) Woutersen, S.; Mu, Y.; Stock, G.; Hamm, P. Chem. Phys. 2001, 266, 137-147. (36) Hertz, H. G.; Zeidler, A. M. Nuclear Magnetic Relaxation in Hydrogen Bond Liquids, in The Hydrogen Bond, Shuster, P., Zundel, G., Sandorfy, C., Eds.; North-Holland: Amsterdam, 1976; Vol. III. (37) E. Breitmaier, E.; Spohn, K. H.; Berger, S. Angew. Chem., Int. Ed. Engl. 1975, 14, 144-159. (38) Seipelt, C. G.; Zeidler, M. D. Ber. Bunsen-Ges. Phys. Chem. 1997, 101, 1501-1508.