Molecular Dynamics Simulations of Aqueous Solutions of

Jul 13, 2006 - Department of Chemistry and Courant Institute of Mathematical Sciences, New York UniVersity,. New York, New York 10003. ReceiVed: March...
2 downloads 0 Views 158KB Size
Download PDF
14652

J. Phys. Chem. B 2006, 110, 14652-14658

Molecular Dynamics Simulations of Aqueous Solutions of Ethanolamines Roberto Lo´ pez-Rendo´ n, Marco A. Mora, and Jose´ Alejandre* Departamento de Quı`mica, UniVersidad Auto´ noma Metropolitana-Iztapalapa, AVenida San Rafael Atlixco 186 09340, Me´ xico D.F., Mexico

Mark E. Tuckerman Department of Chemistry and Courant Institute of Mathematical Sciences, New York UniVersity, New York, New York 10003 ReceiVed: March 29, 2006; In Final Form: June 1, 2006

We report on molecular dynamics simulations performed at constant temperature and pressure to study ethanolamines as pure components and in aqueous solutions. A new geometric integration algorithm that preserves the correct phase space volume is employed to study molecules having up to three ethanol chains. The most stable geometry, rotational barriers, and atomic charges were obtained by ab initio calculations in the gas phase. The calculated dipole moments agree well with available experimental data. The most stable conformation, due to intramolecular hydrogen bonding interactions, has a ringlike structure in one of the ethanol chains, leading to high molecular stability. All molecular dynamics simulations were performed in the liquid phase. The interaction parameters are the same for the atoms in the ethanol chains, reducing the number of variables in the potential model. Intermolecular hydrogen bonding is also analyzed, and it is shown that water associates at low water mole fractions. The force field reproduced (within 1%) the experimental liquid densities at different temperatures of pure components and aqueous solutions at 313 K. The excess and partial molar volumes are analyzed as a function of ethanolamine concentration.

1. Introduction Ethanolamines in aqueous solutions are widely used in the manufacture of cosmetic, pharmaceuticals, surface-active agents, and insecticides. They are also of great interest as solvents used in natural gas sweetening to separate acid gases such as carbon dioxide and sulfide acid from natural gas.1 The knowledge of the thermodynamic properties of these solutions is required for engineering design and subsequent operation. Ethanolamines are highly polar molecules that form strong hydrogen bonds. Knowledge of the properties of these mixtures is useful in connection with industrial applications. Derived excess molar volumes can be used as a basis for understanding some of the molecular interactions in water-organic mixtures. The solubility of carbon dioxide and hydrogen sulfide have been measured experimentally in aqueous solutions of highly polar molecules having different structure.2,3 The smallest molecule is 2-amino ethanol (MEA; NH2-(CH2)2-OH), which has one ethanol and one amine group. Other ethanolamines such as diethanolamine (DEA) and triethanolamine (TEA) are obtained by replacing an ethanol chain by one or both hydrogens in the nitrogen atom, respectively. Methyldiethanolamine (MDEA) contains two ethanol chains and one methyl group bonded to the nitrogen atom. Figure 1 shows the basic structure of the class of molecules studied in this work; the definitions of R1, R2, and R3 are given in Table 1. Theoretical calculations have mainly been performed on MEA.4-9 Quantum chemical calculations of basicities and proton affinities have also been recently reported.7 However, because these calculations yielded poor agreement with experimental * Corresponding author. E-mail: [email protected].

Figure 1. Schematic representation of ethanolamines.

TABLE 1: Notation Used to Distinguish Different Ethanolamines molecule

R1

R2

R3

MEA DEA TEA MDEA

-C2H5OH -C2H5OH -C2H5OH -C2H5OH

-H -H -C2H5OH -C2H5OH

-H -C2H5OH -C2H5OH -CH3

data, the differences were ascribed to the failure of the level of theory employed to properly treat the hydrogen-bonding interactions. There have been very few reported studies aimed at analyzing the properties of larger ethanolamines.10 The structural effects of the hydrogen bond in MEA have been analyzed by microwave spectroscopy,11 where it was established that the gauche form of the molecule is stabilized by a hydrogen-bonding interaction of the type O-H‚‚‚N. On the theoretical side, density functional calculations on monomers and dimers of MEA have been reported.9 These calculations, which used natural bond

10.1021/jp0619540 CCC: $33.50 © 2006 American Chemical Society Published on Web 07/13/2006

Computer Simulations of Ethanolamines orbitals, yielded geometries and hydrogen-bond formation energies that were found to be in good agreement with experimental data. A conformational analysis using ab initio methods was performed on MEA and other molecules having intramolecular hydrogen-bonding configurations.12-15 It is generally accepted that the intramolecular hydrogen bond in MEA plays an important role in its conformational landscape. The molecular geometries of DEA, TEA, and MDEA were reported by da Silva et al.10 using highly accurate ab initio methods. Button, et al.16,17 were the first to simulate MEA in the liquid phase using a simplified model to calculate the electrostatic interactions. A force field of MEA was obtained by Alejandre, et al.18 using ab initio calculations and molecular dynamics simulations at the liquid-vapor interface. In that work, the electrostatic interactions were treated more accurately via Ewald summation. The calculated coexistence densities and surface tension were found to be in good agreement with experimental results. Simulations of MEA in the liquid phase have also been performed by Gubskaya and Kusalik19 using different force fields to calculate the structure of the fluid, the heat of vaporization, and the self-diffusion coefficient. Recently Gubskaya and Kusalik20 also performed simulations of MEA in aqueous solutions to analyze the hydrogen-bonding interactions by calculating the radial distribution function (rdf). They found that water associates in water-deficient solutions. To the best of our knowledge, no simulation studies have been performed for molecules having more ethanol groups. The main interest of this work is to develop a force field for ethanolamine molecules having up to three ethanol chains that is able to reproduce their liquid densities as a function of temperature for both pure components and aqueous solutions at different ethanolamine mole fractions. The temperature range is chosen to correspond to those of certain industrial applications such as natural gas sweetening. Molecular dynamics calculations performed using this force field revealed intramolecular hydrogenbonding interactions in all the molecules leading to stable ring structures, with at least one of the ethanol chains being in the gauche structure. In molecules with two or more ethanol groups, the chains are found to have similar geometry. The remainder of the paper is organized as follows: The ab initio calculations used to develop the force field parameters are described in Section 2. Section 3 describes the force field used in the molecular dynamics simulations. The simulation details in the isothermal-isobaric ensemble are presented in Section 4. The simulation results of liquid densities, excess molar volumes, and the rdf’s of pure components and aqueous solutions are discussed in Section 5. Conclusions and references are given in Section 6. 2. Ab Initio Calculations Ab initio molecular orbital calculations were carried out with the Gaussian 98 suite.21 In what follows, we denote the molecular structure according to the notation given in Table 1. One or both hydrogens bonded to the nitrogen atom in MEA are replaced by ethanol or methyl groups to build up larger ethanolamines. Because the molecules of MEA, DEA, TEA, and MDEA can adopt more than one stable conformation or local minimum on the potential energy surface, a systematic optimization process was adopted beginning at the AM1 level of calculation. From different initial geometries, all the possible conformers were generated until the global minimum energy conformation was obtained. It is well-known that semiempirical methods are good enough to reproduce experimental geom-

J. Phys. Chem. B, Vol. 110, No. 30, 2006 14653

Figure 2. Most stable geometry of ethanolamines. These structures were optimized at the HF/6-311G++(2d,2p) level of theory. The distances shown between H‚‚‚N atoms and H‚‚‚O correspond to hydrogen bonds.

TABLE 2: Parameters of the Most Stable Geometry of Ethanolamines Obtained Using ab Initio Calculations at the HF/6-311G++(2d,2p) Levela expb

MEA

DEA

TEA

MDEA

H-O O-C C-C C-N N-R2 N-R3

1.000 1.396 1.526 1.475 1.017 1.017

bond distance 0.942 0.941 1.396 1.397 1.517 1.516 1.457 1.451 0.998 0.998 0.996 1.448

H-O-C O-C-C C-C-N C-N-R2 C-N-R3 R2-N-R3

108.0 112.1 108.1 110.4 111.3 109.9

bond angles 108.1 108.1 111.9 112.0 109.8 110.2 111.3 109.8 111.7 114.1 107.5 108.5

108.6 112.3 113.0 113.5 113.5 113.5

108.6 111.9 112.8 111.9 111.2 112.0

H-O-C-C O-C-C-N C-C-N-R2 C-C-N-R3

-27.0 55.4 78.2 -159.5

dihedral angles -45.20 -44.76 59.82 60.78 76.01 66.75 -163.78 -171.20

-74.62 54.00 76.02 -152.31

-52.85 64.05 82.56 -151.42

0.940 1.400 1.522 1.452 1.452 1.452

0.940 1.400 1.518 1.460 1.457 1.452

a The bond distances are in Å, the bond and torsional angles are given in degrees. b Experimental data of MEA taken from ref 11.

etries;22 however, a reoptimization of all the geometrical parameters with the Hartree-Fock method at the 6-31G* level of theory was performed to improve the electronic structure. The final geometry, atomic charges, and dipole moments were calculated with the HF/6-311++G(2d,2p) level of theory for all molecules. Figure 2 shows the most stable geometry of the molecules studied in this work. The optimized geometries of MEA are in good agreement with those reported by other authors.4-9 The calculated H-O-C-C dihedral angle of MEA, -45.2, differs from the experimental value, -27 ( 6, probably because of the low precision in the experimental determination of this parameter.9 For larger ethanolamines our optimized geometries agree well with those of da Silva, et al.,7,10 who employed higher levels of theory. The similarity of the results justify the more approximate level of theory used in this work. The geometry of the most stable conformer of all ethanolamines is given in Table 2, where a comparison with the experimental data of MEA11 is also given. The main differences lie in the

14654 J. Phys. Chem. B, Vol. 110, No. 30, 2006

Lo´pez-Rendo´n et al.

TABLE 3: Nonbonding Distances in Å of Ethanolamines at the HF/6-311G++(2d,2p) Level MEA

DEA

TEA

MDEA

distance

R1

R1

R2

R1)R2)R3

R1

R2

r(O- - -N) r(OH- - -N) r(NH- - -O) r(OH- - -O)

2.86 2.40a 3.37 4.85

2.87 2.42a 3.27 6.17

2.85 3.66 2.46a 2.90

2.88 2.67a

2.96 2.56a

2.93 2.53a

3.52

3.60

a

These distances correspond to hydrogen bond formation according to ref 23.

Figure 4. Rotational energy barriers of ethanolamines on the OCCN bond of the ethanol chain denoted as R1.

Figure 3. Absolute entropy as a function of the temperature of the ethanolamines of one molecule in the ideal gas approximation obtained using the package Gaussian 98 at 1 atm. The results correspond to the molecular geometry given in Figure 2.

bond distances N-R2 and N-R3, where R2 and R3 denote the ethanol or CH3 groups. The bond and torsional angles have approximately the same values as those of MEA. The torsional angle OC-CN in all ethanolamines is approximately 60° and favors a ringlike structure in the HO-CC-N group due to an intramolecular hydrogen bond. The criterion used in this work for the existence of hydrogen bonds is that the distance r(H‚‚ ‚B) should be smaller than the sum of the van der Waals radii of H and B atoms. Using the tabulated values by Bondi23 (rH ) 1.2 Å, rO ) 1.4 Å, and rN ) 1.5 Å, where ri is the radius of atom i), a hydrogen bond is defined when r(H‚‚‚N) < 2.7 Å and r(H‚‚‚O) < 2.6 Å. Table 3 gives the intramolecular bond distances for specific atom pairs. It is shown that, according to ref 23, the terminal hydrogen of the ethanol chain forms an intramolecular hydrogen bond with the nitrogen atom. Figure 2 also shows these distances for MEA, DEA, and MDEA molecules. The parameters given in Table 2 for the most stable geometry suggest that ethanol chains are equivalent in different molecules. This conclusion can also be drawn from Figure 3, where it is seen that the change in absolute entropy between DEA and MEA, ∆S° ) S° (DEA) - S° (MEA), at a given temperature is the same as that between TEA and DEA, that is, ∆S° ) S° (TEA) - S° (DEA). It can be seen, for example, that ∆S° ) 20 cal mol-1 K-1 and 45 cal mol-1 K-1 at 300 and 1000 K, respectively. The values of S° in MDEA are larger than those in DEA because the methyl group has larger contributions to translations and vibrations than the hydrogen in DEA. The absolute entropy S° of MEA from this work is in good agreement with those reported by Chang et al.15 The charges were obtained using the electrostatic potential, yielding calculated dipole moments for MEA, DEA, and TEA of 3.04, 3.00, and 3.50 D, respectively. The corre-

sponding experimental values11 are 3.05 D for MEA, and 2.80 and 3.60 D24 for DEA and TEA, respectively. The largest difference is 7% in DEA. The dipole moment of MDEA obtained in this work is 3.32 D, for which we were unable to find an experimental value. By taking the global energy minimum conformation as the initial point, all rotational barriers (HO-CC, OC-CN, and CCNR, where R denotes H, CH3, or an ethanol chain) were calculated at the HF/6-311++G(2d,2p) level for all target molecules. In each case, the rotation angle was varied by increments of 30° from 0 to 360°. The rotational barriers for the OC-CN obtained using the present ab initio protocol are shown in Figure 4 for all molecules studied in this work. The most stable conformer was found at 60°. These results also support the idea that, in the ethanolamines, the ethanol chains are equivalent. It should be noted that our rotational barrier NC-CO of MEA is in good agreement with that reported by Chang, et al.,15 although these authors performed their calculations at the more demanding MP2 level of theory. This suggests that it is not necessary to use such high-level calculations to have a reasonable description of the rotational barriers in molecules having strong intramolecular interactions. Consequently, the barriers in other molecules were calculated using the same protocol as that of MEA. 3. Force Field The force field used in this work to perform the molecular dynamics simulations has intra- and intermolecular contributions. The bond bending and bond stretching interactions are modeled as harmonic potentials

U(rij) )

kr (r - r0)2 2 ij

U(θijk) )

kθ (θ - θ0)2 (3.1) 2 ijk

where rij ) |ri - rj| is the distance between atoms i and j, ri is the position of atom i, and θijk is the bond angle between the vectors ri - rj and rk - rj; the subscript 0 is used to denote an equilibrium value, and kr and kθ are the corresponding spring constants. The torsion potential between four atoms is written as in ref 25

Computer Simulations of Ethanolamines

J. Phys. Chem. B, Vol. 110, No. 30, 2006 14655

5

U(φijkl) )

∑ Cn cosn φijkl n)0

(3.2)

where i, j, k, and l indicate the four atoms in terms of which the dihedral angle is defined, and Cn are constants obtained by fitting the torsional barriers to the ab initio calculations described in the previous section. The trans conformation is defined when φijkl ) 180°. The Cn coefficients for the dihedral angles HOCC, OC-CN, and CC-NR2 were obtained using results of DEA and those of MDEA were used for CC-NR3. Figure 5 shows the rotational barriers used to simulate the ethanolamine molecules. The rotational barrier of OC-CN has stable conformations at (60° and 180°, similar to linear hydrocarbons. The improper torsional angles that involved nitrogen were calculated using a harmonic function,

U(φijkl) )

kφ (φ - φ0)2 2 ijkl

(3.3)

where kφ is the spring constant and φ0 is the torsional angle at equilibrium. The interaction between any two atoms separated by more than three bonds is given by Lennard-Jones and Coulomb forms:

U(rij) )

[( ) ( ) ]

σij 1 qiqj + 4ij 4π0 rij rij

12

-

σij rij

6

(3.4)

where qi is the charge of atom i, and 0 is the vacuum permittivity. The diameter of the united atom is σ, whereas  measures the strength of the interaction. Interactions between unlike species are determined by the Lorentz-Berthelot mixing rules: σij ) (σii + σjj)/2 and ij ) (iijj)1/2. To simplify the force field, we considered that the geometry and torsional barriers of ethanol chains are the same in all molecules, thereby reducing the number of parameters. The nitrogen and the atoms bonded to it have different charge and Lennard-Jones parameters. Tables 4 and 5 give the potential parameters for the intra- and intermolecular interactions, respectively. The parameters for interactions between atoms separated by more than three bonds in a molecule were not scaled, that is, they retain the same value as those used to calculate the intermolecular forces. To simulate water, a flexible model with an oxygen diameter slightly larger than that of the SPC/E model was used, following Tironi, et al.26 This model gives a density 0.6% greater than the experimental value at 1 atm and 313 K. The parameters for the intramolecular interactions are kr ) 557733.9 K Å-2, kθ ) 46067.4 K rad-2, r0 ) 1.0 Å, and θ0 ) 109.4° and are the same as that of Tironi, et al.26 used to simulate flexible water using the SPC/F model. The Lennard-Jones parameters and charges for water are also given in Table 4. 4. Simulation Details The force field was validated by computing properties in the liquid phase. Molecular dynamics simulations in the isothermalisobaric ensemble were carried out to obtain the liquid density in one-component systems at different temperatures as well as in aqueous solution. A new integrator,27 which preserves the phase space volume, was used with two Nose´-Hoover chains28 coupled to the system: one to the particles and the other to the barostat. The method is based on a new approach applied to the original Martyna, Tobias, and Klein equations of motion.29 The total number of molecules was 500 for simulations of both the pure liquids and aqueous solutions. The molecules were

Figure 5. Rotational energy barriers of ethanolamines obtained by fitting the ab initio calculations to the function given by eq 3.2. We show the four most important barriers used in the simulations.

TABLE 4: Intramolecular Force Field Parameters for Ethanolamines Used in the Molecular Dynamics Simulationsa bond distance

r0 (Å)

kr (K Å-2)

H-O O-C C-C C-N H-N

0.96 1.42 1.54 1.44 1.00

102608 102608 102608 102608 102608

θ0

kθ (K rad-2)

108.0 110.5 114.5 111.2 112.8 108.48

71869.9 131550.9 121762.0 76164.4 64796.6 64796.6

bond angle H-O-C O-C-C C-C-N C-N-H C-N-C H-N-H torsion

C0

C1

C2

C3

C4

C5

H-O-C-C O-C-C-N C-C-N-R2 C-C-N-R3

503.0 1775.2 1385.2 1822.5

-275.6 -3746.4 -478.9 -1746.9

524.5 -1081.6 -2136.1 -258.1

1352.6 1866.7 593.6 2871.5

196.2 1261.3 4015.3 891.1

115.0 2298.7 -50.1 995.7

molecule DEA MDEA TEA

improper torsion

kφ (K rad-2)

φ0

C-C-N-R3 C-C-N-R2 C-C-N-R2 C-C-N-R3

92608 92608 92608 92608

-176.6 -158.9 -139.4 72.71

a The angles are given in degrees, and the constants for regular torsions are given in K.

arranged in a face-centered cubic (fcc) crystalline array at low density (F ) 0.2 g/cm3) to avoid high energy configurations due to overlaps between molecules. The Lennard-Jones and Coulomb interactions were handled via Ewald summation30-32 to properly evaluate the long-range interactions. The cutoff distance Rc ) 9.5 Å was used to calculate the real part of the potential for all the simulations. The R parameters in the Ewald sums were 0.27 and 0.29 Å-1 for the real part of the LennardJones and Coulomb potentials, respectively. The reciprocal-space term in both models was calculated using the smooth particle mesh Ewald method.33 For the electrostatic interactions, the spline order was set at 6, and the grid density was set at 0.5 Å in every direction, while, for the Lennard-Jones interactions,

14656 J. Phys. Chem. B, Vol. 110, No. 30, 2006

Lo´pez-Rendo´n et al.

TABLE 5: Intermolecular Force Field Lennard-Jones Parameters and Charges for Sites in Ethanolamines and Watera,b site

H2O

MEA

DEA

MDEA

TEA

Hw 0.4238 Ow -0.8476 H-(O) 0.369 0.369 0.369 0.369 O -0.660 -0.660 -0.660 -0.66 -CH2-(O) 0.291 0.291 0.291 0.291 -CH2-(N) 0.250 0.114 0.030 -0.170 -N-0.938 -N-0.572 -N-0.09 -N0.510 H-(N) 0.344 0.344 0.030 CH3

σ/Å

/kB

0.0 3.182 0.0 3.07 3.905 3.905 3.250 4.600 3.25 3.050 0.980 4.050

0.00 78.2 3.90 85.51 59.42 59.42 85.51 85.51 85.51 85.51 3.90 59.42

a Reference 26. b The charges are in electron units. The subscript “w” is for atoms in water.

Figure 6. Densities of liquid ethanolamines and water as a function of temperature at 1 atm. The lines correspond to experimental values taken from ref 39 for ethanolamines and from ref 40 for water. The symbols are results from the simulations of this work.

the parameters were 5 and 1 Å, respectively. The external pressure was 1 atm. The thermostat and barostat parameters were τp ) 20∆t and τb ) 5∆t in all the simulations. A multiple timestep algorithm27,28,34,35 was used with the smallest time step δt ) 0.3 fs. The largest time step ∆t ) 6δt was used to evaluate the reciprocal-space interactions. Thus, each large time step involves six intramolecular force evaluations, two real-space force evaluations, and only one reciprocal-space evaluation. The equilibrium was reached after 105 of the largest time steps. The liquid densities and rdf’s were obtained by averaging 40 independent blocks of 5000∆t for pure components and binary mixtures. The calculated average temperature and pressure were the same as the imposed values, and no drift in the conserved quantity was observed. 5. Simulation Results The liquid densities as functions of temperature for the ethanolamines and water are shown in Figure 6. The temperature range has been chosen so as to correspond with specific industrial applications, such as pharmaceutical manufacture and natural gas sweetening.2 The calculated values are in good agreement with experimental data39 at all temperatures, where the difference between both sets of data is less than 1%. The largest difference is found in TEA at high temperatures. The force fields of DEA and MDEA described herein were used to perform simulations at 313 K on aqueous solutions at

Figure 7. Densities of aqueous solutions of DEA and MDEA at 313K and 1 atm. The lines correspond to experimental values taken from ref 41 for DEA and ref 38 for MDEA. The filled circles and filled squares are results from the simulations of this work for DEA and MDEA, respectively. The error bars in densities are shown.

Figure 8. Radial distribution function of selective hydrogen bonds between MDEA and water: (a) Hw-Ow, (b) H-O, (c) Hw-O, (d) H-Ow. The symbols with subscript “w” are atoms that belong to water; otherwise they are atoms of MDEA. The dotted and long dashed lines are for mixtures having X2 ) 0.25 and X2 ) 0.75 of MDEA, respectively. The solid lines in panels a and b are for pure components.

different ethanolamine mole fractions, X2. Figure 7 shows the calculated liquid densities compared with the experimental data.41,38 The difference between both set of data is less than 1% at all ethanolamine concentrations. Ethanolamines are molecules that form strong intramolecular and intermolecular hydrogen bonds. The intramolecular interactions were described previously in section 2. In pure components, the important hydrogen bonds are between the hydrogen and oxygen atoms. In mixtures, additional hydrogen bonds are found between ethanolamines and water. Figure 8 shows typical rdf’s of aqueous solutions of MDEA at different concentrations. In addition to results for the pure components, results for the mixtures having compositions of X2 ) 0.25 and X2 ) 0.75 are also shown. The results for the other ethanolamines follow the same trend as that found in MDEA. The hydrogen bond between water molecules is shown in Figure 8a. It is seen that, as the concentration of water decreases, the first two peaks in the rdf increase as a result of association between water molecules. It is well-known that, at low density, water associates, forming large voids.36 In this case, no voids are found; however, water seems to behave in qualitatively the same manner. To check

Computer Simulations of Ethanolamines

J. Phys. Chem. B, Vol. 110, No. 30, 2006 14657

whether the water associates, we performed a simulation of pure water at constant temperature with the initial configuration generated by removing the MDEA molecules from the solution at X2 ) 0.75. For this water simulation, the first two peaks of the rdf of hydrogen-oxygen interactions were found in the same position as those shown in Figure 8a; however, the maxima in the rdf were found to be 42 and 29, respectively, well above that found in solution at the same water mole fraction. The arrangement between water molecules in solution seems to follow the same behavior as that found in pure water at low concentrations. This finding is in agreement with simulation results of aqueous solutions of MEA20 at low concentrations of water. The rdf’s between the hydrogen of one molecule and the oxygen from another in the ethanol chains in MDEA are shown in Figure 8b. The strength in the hydrogen bond as the concentration of MDEA increases is opposite that found in water; in this case, the intensity is higher in the pure liquid than in the mixture. Figure 8c shows the rdf’s between the oxygen in water and the hydrogen in MDEA, and Figure 8d shows the correspondent results for the rdf’s between the oxygen of MDEA and the hydrogen of water for two aqueous solutions. In both cases, the strength of the hydrogen bond increases with the composition of MDEA. The excess molar volume of the mixtures VEXC was calculated using the standard definition

VEXC ) Vm - (X1V /m1 + X2V /m2)

(5.5)

where Vm ) (X1M1 + X2M2)/F is the molar volume of the / mixture, Xi is the molar fraction, V mi is the molar volume of the pure component, and Mi is the molar mass of component i. The density of the mixture is denoted by F and was obtained during the simulation in the isothermal-isobaric ensemble. The calculated results of DEA and MDEA are compared with experimental data from ref 41 and refs 37 and 38, respectively. In the latter case, the reported results are different at moderate and high concentrations of MDEA, showing the difficulty of measuring VEXC for this highly polar molecule. Figure 9a shows the excess molar volume of aqueous solutions studied in this work. The calculated VEXC shows negative deviations from ideality at all concentrations of DEA and MDEA, in agreement with experimental data. Except at high concentrations of ethanolamine, the simulation results are overestimated. We found that very small changes in liquid densities result in large deviations in VEXC. The excess molar volumes from simulations of this work were fitted to a Redlich-Kister equation:38 5

VEXC ) X2(1 - X2)

An(1 - 2X2)n-1 ∑ n)0

(5.6)

The values of the fitted parameters An for DEA were A0 ) -1.445, A1 ) -0.334, A2 ) 1.198, A3 ) 0.603, A4 ) -1.677, and A5 ) -1.229, and, for MDEA, they were A0 ) -5.333, A1 ) -0.686, A2 ) 7.647, A3 ) 0.699, A4 ) -7.373, and A5 ) -0.464). We also calculated the partial molar volumes,38 V h m2, of DEA and MDEA by using

∑n An (1 - 2X2)n-1 + 2X2(1 - X2)2∑An (n - 1)(1 - 2X2)n-2 n

V h m2 ) V /m2 + (1 - X2)2

(5.7)

Figure 9. Properties of aqueous solutions of ethanolamines: (a) excess molar volume and (b) partial molar volume. The dot-dashed line is for the experimental results (ref 41) of DEA, and the solid and dashed lines are the experimental results (refs 37 and 38) of MDEA. The symbols are the simulation results obtained in this work, and they have the same meaning as those in Figure 7. The unit of volume is cm3 mol-1.

The results of V h m2 for DEA and MDEA are shown in Figure 9b and compared with the experimental results.37,38,41 The VEXC values from ref 37 were fitted to a Redlich-Kister equation of order 5, instead of taking their values of An fitted up to order 3. Except at moderate concentrations of ethanolamine (0.1 < X2 < 0.3), the agreement is satisfactory for the other concentrations. The major difference is found in MDEA at X2 ) 0.25, where V h m2 is 2.4% greater than the experimental value. The calculated partial molar volumes of ethanolamines at infinite dilution, V h 02 3 -1 are 93.89 and 110.58 cm mol for DEA and MDEA, respectively. The corresponding experimental results41,38 are 94.6 and 110.7 cm3mol-1. 6. Conclusions By using ab initio calculations and molecular dynamics simulations at constant temperature and pressure, we developed a force field of ethanolamines having up to three ethanol chains. The results were obtained in a temperature range relevant for industrial applications. In all the molecules, the gauche conformation of one ethanol chain was obtained and attributed to a strong intramolecular hydrogen bond between hydrogen and nitrogen atoms. The interaction parameters of the ethanol chains are transferable from one chain to another. The main difference between ethanolamine molecules is found in the nitrogen atom where the charge changes from negative to positive values from MEA to TEA, respectively. All ethanolamine molecules have intermolecular hydrogen bonding higher than water, both as pure components and in aqueous solution. Water associates in waterdeficient mixtures in a fashion similar to that found in pure water at low density. The force field of pure ethanolamines and their mixture with water produces liquid densities and partial molar volumes in good agreement with experimental data. This suggests that a polarizable force field model might not be

14658 J. Phys. Chem. B, Vol. 110, No. 30, 2006 required to simulate this class of complex polar molecules in the liquid phase. The excess molar volume calculated in this work is slightly overestimated in DEA at all concentrations and in MDEA at low concentrations of ethanolamine. The difference is of the same order of magnitude as that found in experiments at high concentrations of MDEA. Simulations at the liquidvapor interface are currently being performed for DEA and MDEA molecules. Acknowledgment. R.L.R. thanks Instituto Mexicano del Petro´leo (IMP) and Conacyt for a scholarship. J.A. and M.A.M. thank Conacyt for financial support. J.A. also thanks Tom Darden for helpful discussions about the SPME method. References and Notes (1) Amararene, F.; Balz, P.; Bouallou, C.; Cadours, R.; Lecomte, F.; Mougin, P.; Richon, D. J. Chem. Eng. Data 2003, 48, 1565. (2) Silkenbaumer, D.; Rumpf, B.; Lichtenthaler, R. N. Ind. Eng. Chem. Res. 1998, 37, 3133. (3) Kamps, A. P.-S.; Balaban, A.; Jo¨decke, M.; Kuranov, G.; Smirnova, N. A.; Maurer, G. J. Eng. Chem. Res. 2001, 40, 696. (4) Raudino, A.; Millefiori, S.; Zuccarello, F.; Millefiori, A. J. Mol. Struct. (THEOCHEM) 1979, 51, 295. (5) Vanquickenborne, L. G.; Coussens, B.; Verlinde, C.; de Ranter, C. J. Mol. Struct. (THEOCHEM) 1989, 201, 1. (6) Kelterer, A.-M.; Ramek, M. J. Mol. Struct. (THEOCHEM) 1991, 232, 189. (7) da Silva, E. F. J. Phys. Chem. A 2005, 109, 1603. (8) Silva, C.; Duarte, M. L.; Fausto, R. J. Mol. Struct. 1999, 482, 591. (9) Vorobyov, I.; Yappert, M. C.; Dupre´ D. B. J. Phys. Chem. A 2002, 106, 668. (10) da Silva, E. F.; Svendsen, H. F. Ind. Eng. Chem. Res. 2003, 42, 4414. (11) Penn, R. E.; Curl, R. F., Jr. J. Phys. Chem. 1971, 55, 651. (12) Buemi, G. Int. J. Quantum Chem. 1996, 59, 227. (13) Kelterer, A.-M.; Ramek, M.; Regina, F. F.; Cao, M.; Scha¨fer, L. J. Mol. Struct. (THEOCHEM) 1994, 310, 45. (14) Chang, Y.-P.; Su T.-M. J. Mol. Struct. (THEOCHEM) 1996, 365, 183. (15) Chang, Y.-P.; Su, T.-M.; Li, T.-W.; Chao, I. J. Phys. Chem. A 1997, 101, 6107. (16) Button, J. K.; Gubbins, K. E.; Tanaka, H.; Nakanishi, K. Fluid Phase Equilib. 1996, 116, 320. (17) Button, J. K.; Gubbins, K. E. Fluid Phase Equilib. 1999, 158, 175. (18) Alejandre, J.; Rivera, J. L.; Mora, M. A.; de la Garza V. J. Phys. Chem. B 2000, 104, 1332. (19) Gubskaya, A. V.; Kusalik, P. G. J. Phys. Chem. A 2004, 108, 7151. (20) Gubskaya, A. V.; Kusalik, P. G. J. Phys. Chem. A 2004, 108, 7165.

Lo´pez-Rendo´n et al. (21) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; AlLaham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A7; Gaussian, Inc.: Pittsburgh, PA, 1998. (22) Dewar, M. J.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. P J. Am. Chem. Soc. 1985, 107, 3902. (23) Bondi, A. J. Phys. Chem. 1964, 68, 441. (24) Lide, D. R. Handbook of Chemistry and Physcis; Springer-Verlag: New York, 1997; p 15. (25) Ryckaert, J.-P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1978, 66, 95. (26) Tironi, I. G.; Brunne, F. M.; van Gunsteren, W. F. Chem. Phys. Lett. 1996, 250, 19. (27) Tuckerman, M. E.; Alejandre, J.; Lo´pez-Rendo´n, R.; Jochim, A. L.; Martyna, G. J. J. Phys. A: Math. Gen. 2006, 39, 5629. (28) Martyna, G. J.; Tuckerman, M. E.; Klein, M. L. J. Chem. Phys. 1992, 97, 2635. (29) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177. (30) Karasawa, N.; Goddard, W. A., III J. Phys. Chem. 1989, 93, 7320. (31) Lo´pez-Lemus, J.; Alejandre, J. Mol. Phys. 2002, 100, 2983. (32) Lo´pez-Lemus, J.; Romero-Bastida, M.; Darden, T. A.; Alejandre, J. Mol. Phys., in press. (33) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (34) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117. (35) Tuckerman, M. E.; Martyna, G. J.; Berne, B. J. J. Chem. Phys. 1992, 97, 1990. (36) Mountain, R. D. J. Chem. Phys. 1999, 110, 2109. (37) Bernal-Garcia, J. M.; Ramos-Estrada, M.; Iglesias-Silva, G. A.; Hall, K. R. J. Chem. Eng. Data 2003, 48, 1442. (38) Maham, Y.; Teng, T. T.; Mather, A. E.; Hepler, L. G. Can. J. Chem. 1995, 73, 1514. (39) Diguillo, M. R.; Lee, R.-J.; Schaeffer, S. T.; Brasher, L. L.; Amyn, S. T. J. Chem. Eng. Data 1992, 37, 239. (40) Lemmon, E. W.; McLinden, M. O.; Friend, D. G. Thermophysical Properties of Fluid Systems, NIST Chemistry WebBook, NIST Standard Reference Database; Linstrom, P. J., Mallard, W. G., Eds; U.S. Department of Commerce: Washington, DC, 2005. http://webbook.nist.gov. (41) Maham, Y.; Teng, T. T.; Hepler, L. G.; Mather, A. E. J. Solution Chem. 1994, 23, 195.