Molecular dynamics simulations of the aqueous solution of tetramethylammonium chloride

18 April 1997

CHEMICAL PHYSICS LETTERS

ELSEVIER

Chemical Physics Letters 268 (1997) 325-330

Molecular dynamics simulations of the aqueous solution of tetramethylammonium chloride E. Hawlicka, T.

Dlugoborski

Institute of Applied Radiation Chemistry, Technical University, Zeromskiego 116, 90-924 Lodz, Poland

Received 28 June 1996; in final form 25 February 1997

Abstract A molecular dynamics simulation of 0.555 molal tetramethylammonium chloride in aqueous solution has been performed at 293 K. The tetramethylammonium cations (TMA) are treated as charged 'pseudo-atoms' and a flexible 3-site model employed for the water molecules. Structural properties of the system are discussed on the basis of radial distribution functions, orientation of the molecules in the hydration shells of the ions and their geometrical arrangement. The power spectra and average number of H-bonds of a water molecule in the hydration shells of the ions and in bulk solvent were computed. The TMA exhibits 'apolar' and 'ionic' properties and its influence on the structure of bulk water is minor. Association of ions occurs even for moderate concentration and solvent-separated ion pairs are formed.

1. Introduction Hydrophobic phenomena occurring in aqueous solutions of non-polar solutes are believed to be induced by the structural peculiarities of H-bonded networks of water [1]. An important questions is how far the solute perturbs the water structure. Though this information is important for understanding the role of hydration in various processes like protein folding or micelle formation, suggestions have been controversial and perturbations ranging from one or two hydration layers up to hundreds of ~tngstr~Sm have been postulated [2]. X-ray and neutron scattering techniques may provide direct information about the solution structure, but the total pair correlation function obtained from an experiment is the weighted sum of the partial pair correlation functions of atoms, and even for small molecules their number becomes too large to be

interpreted unambiguously. Therefore, perturbations of the water structure due to hydrophobic effects have to be investigated in model systems. Some of the best model systems are aqueous solutions of tetraalkylammonium halides. These salts are highly soluble in water and it is possible to follow the changes in the hydrophobic effects with increasing alkyl chain. Recently, neutron diffraction experiments involving isotopic substitution have shown [3] that the tetramethylammonium cation is more 'apolar' than ' i o n i c ' ; thus it can serve as a model system for studying hydrophobic hydration. This Letter presents the results of a MD simulation of an aqueous solution of tetramethylammonium chloride (TMACI). The aim was to gain insight into structures of the hydration shells of T M A cations and bulk water and to investigate the dynamic properties of water in various environments: the solvation shells of ions and the bulk.

000%2614/97/$17.00 Copyright © 1997 Elsevier Science B.V. All fights reserved. PII S0009-261 4(97)00229-7

E. Hawlicka, T. Dlugoborski / Chemical Physics Letters 268 (1997) 325-330

326

2. Potential model

Rigid models are commonly employed to describe intermolecular interactions and they correctly reproduce many measurable quantities of liquid water such as radial distribution functions, internal energy, pressure and relative permittivity. Modifications, which introduce flexibility and polarizability of the water molecule are not commonly adopted and their necessity is in dispute [4]. Flexible models may, however, be better suited to simulating the properties of water in the hydration layers. The flexible models include intramolecular interactions and they permit changes in the molecular geometry, bond length and valence angle. Thus intramolecular vibrations can be computed from the MD simulation results and the power spectra can be compared with the spectroscopic data, which are frequently used to investigate hydration phenomena. Most flexible models of water are based on SPC [5] or MCY [6] potentials. They are employed to describe intermolecular interactions whereas various spectroscopic potentials [4,7-10] are added for intramolecular interactions. Recently, Ferguson developed a simple flexible model [11] based on the SPC potential, which reproduces many properties of water including intramolecular vibrations. In the simulation presented here, however, the water molecules have been described via the flexible BJH model [12], though it is more complicated than Ferguson's model. The BJH model reproduces fairly well the structural and dynamical properties of liquid water and of aqueous solutions of electrolytes [13]. This model was used previously to investigate solutions of NaCI Table 1 Parameters (in k J / m o l ) for i o n - w a t e r and i o n - i o n potentials (for distances expressed in A.). The potential takes the form:

V(ri)) = --÷QiJ Ai) ri) ri7 10

TMA TMA TMA TMA CI CI CI -

O H TMA CIO H CI

B,j r6 2

Qij

- 9.17447 4.58723 13.8957 - 13.8957 9.17447 - 4.58723 13.8957

10

8

Ai j

4.38046 = 0 246.789 19.2841 0.173801 = 0 1.08931

10-6

Bi )

7.79638 = 0 76.7566 14.8274 1.07322 --- 0 2.43543

in water and methanol-water mixtures [14-16]. Thus the comparison of MD results obtained for aqueous solutions of TMACI and NaC1 may be helpful for understanding hydrophobic effects. The TMA cations are treated as charged 'pseudoatoms'. All their interactions and those of the chloride ion have been expressed as a sum of the Coulomb and Lennard-Jones potentials. The non-Coulombic parameters for the ion-oxygen interactions are summarized in Table 1. The non-Coulombic interactions between the ions and water hydrogens have been neglected.

3. Details of the simulation

In the MD simulation of an 0.555 molal aqueous solution of TMAC1 the cubic box contained 400 molecules of water, 4 cations and 4 anions. The box size, 23.33 ,~, was calculated from the experimental density of 0.9997 g / c m 3 [17] at 293 K. An Ewald summation was employed for all the Coulombic interactions, whereas for the non-Coulombic interactions the shifted force potential method was used [18]. The time step was equal to 2.5 × 10-16 s. The system was equilibrated for 5 ps, then the simulation was extended for a further 50 ps, without any rescaling of the velocities. The average temperature of the run was 293 _+ 3 K. The stability of the total energy was better than 0.1%.

4. Results and discussion

Radial distribution functions (RDF's) obtained for the ions and the water sites, oxygen and hydrogen are shown in Fig. 1. The broad peak of the gVMAO(r) function appears between 4.2 and 5.4 ,~. Its highest value, gTMAO(rmax)= 2.5, is observed at about 4.47 A. The peak position and peak height agree quite well with data obtained in neutron diffraction experiments with nitrogen isotopic substitution [3,19-21]. A distance between othe TMA center and the water oxygen atom of 4.7 A, and a peak height, of about 2, have been reported [3]. The shape of the gVMAO(r) and the peak position also coincide with the MC results [22,23]. Integration of the gTMAO(r) function

327

E. Hawlicka, T. Dlugoborski / Chemical Physics Letters 268 (1997) 325-330 I

'

L

'

I

'

1

'

I

'

I

'

[

'

I

'

I

g

........i..

'

~

'

I

1 I

'

2

'

~

'

3

I

'

I

I

'

I

4

'

I

'

I

5

'

I

'

I

6 I

'

'

I

7 '

I

'

I

8 '

I

9 '

I

4

b

3 I

2

vectors; one vector points from the ion to the oxygen and the second vector is the dipole moment of the water molecule. The distribution of the O angle in the hydration shell of TMA is almost uniform, thus there is a lack of the preferred orientation of the water molecule. The O angle shows only a fiat maximum between 60 and 90 ° and a similar feature has been noticed for the hydration layers of the non-polar molecule of methane [25]. The average number of H-bonds per water molecule (nHB) has been computed independently for the three subsystems: hydration shells of cations, anions and bulk water. The H-bond has been defined via the following geometrical criterion: (i) the distance between the oxygens of the H-bonded molecules cannot exceed 3.4 ,~, (ii) the distance between the H atom of the proton donor and the O atom of the acceptor must be shorter than 2.6 A, and (iii) the angle between the vector connecting the oxygens and the intramolecular OH bond cannot exceed 20 °. The analysis has been done in 0.05 ps intervals over the whole simulation run. The average numbers of H-bonds obtained for a water molecule in the hydration shell of TMA cation and in the bulk solvent are (nHB}TMA = 2.43 and (nHB}bul k = 2.65, respectively. The difference is small, less than 10%, particularly when it is compared with the remarkable reduction of the numbers of H-bonds in the hydration shells of Na + ( ( n H B ) N a = 1.5 [27]). The RDFs of the chloride ion and the oxygen and hydrogen atoms are shown in Fig. lb. The gclo(r) and g c m ( r ) functions exhibit well pronounced first peaks centered at about 3.3 and 2.3 ,~, respectively. The peaks appear at distances about 3 - 5 % shorter than with those reported for a 1.1 molal solution of NaCI [13]. The difference may be due to lower concentration of TMAC1. Integration of the gclo(r) up to its first minimum, at 3.8 ,~, yields the average hydration number of CI- 6.6, whereas the average value equal 8.0 has been found for CI- shell in 1.1 m NaC1 solution [13]. The decrease is probably due to the association of ions. Though the salt concentration is moderate a broad peak of gTMACl(r) appears between 4.5 and 5.8 ,~. In order to examine the formation of ion pairs the distances between the cations and anions have been examined in 0.01 ps time intervals over the whole simulation run. Although there were no permanent ion pairs, at each of o

1

i

J ,

0 0

I

1

1

i

2

3

4

i

I

5

i

I

6

,

I

7

,

I

8

,

I

9

r [A]

Fig. 1. Radial distribution functions of ion-oxygen (full line) and ion-hydrogen (dashed line) for cations (upper part) and anions (bottom part).

up to the first minimum, at 5.9 ,~, yields 25 water molecules in the hydration shell of TMA. This agrees quite well with the experimental results, 20-I-2 molecules calculated from the neutron scattering measurements [19-21] and 25 molecules deduced from the deuterium relaxation time [24]. The gTMAH(r) function also exhibits a broad peak between 4.2 and 5.6 .~. As seen from Fig. l a the hydrogen atoms are, on average, at the same distance from the TMA cation as the oxygens. This feature is consistent with the experimental results [3] and it confirms the conclusion that the TMA cation behaves like an apolar species [25,26]. There are additional clues indicating that the TMA ions are more 'apolar' than 'cationic': (i) lack of a preferred orientation of the water molecules in the first hydration shell and (ii) almost the same average number of hydrogen bonds per water molecule in the hydration shell and the bulk. The orientation of the water molecules in the hydration shells has been deduced from the distribution of O angles. The 6) angle is defined by two

328

E. Hawlicka, T. Dlugoborski / Chemical Physics Letters 268 (1997) 325-330

the sampled times at least one ion pair was always identified. Though only one ion pair is formed, it comprises 25% of all the ions in the system and its presence affects the average hydration number of the ions. The analysis of the ion surroundings leads to the conclusions that the hydration shells of 'free' C1- consist of either 8 or 7 water molecules. The anion involved in the ion pair shares the nearest neighbors with the cations and in that case the anion hydration number is only 4. The water molecules in the hydration shells of C1- form an almost linear H-bond with the anion, but they also form two H-bonds with the water molecules from the second layer. The average number of H-bonds computed, according to the criterion defined above, for the CI- shells in the NaC1 and TMACI solutions are the same, (nHB)CI = 2.0 In order to deduce the effect of TMACI on the water structure the goo(r), goH(r) and gHH(r) functions have been calculated and are presented in Fig. 2. As might be expected the presence of TMACI does not affect the distances between the water molecules and the positions of the peaks are the same as in pure water and aqueous NaC1 [15]. The presence of TMAC1 affects, however, the shape and intensity of the goo(r), g o n ( r ) and gHH(r) functions. The first peaks of these RDF's are broader and higher compared with the NaC1 solution. Integration of the goo(r) function up to its first minimum at 3.1 A, yields a number of nearest neighbors equal to 4.1. It is less than reported for pure water, 4.7, and 1.1 molal NaC1 solution, 4.4. Though the average number of H-bonds computed for bulk water is slightly higher ( n H B ) b u t k = 2.65 as compared with the numbers found for pure water and bulk water in NaCI solution, (nHB)pure ~-- 2.5 a n d (nHB)butk=2.53, respectively, the influence of TMACI on the water structure is minor. Total spectral densities of water have been computed via Fourier transformation of the normalized velocity correlation functions. The details of the computations have been reported previously [27]. Separate calculations have been performed for the water molecules belonging to the hydration shells and the bulk solvent. The thickness of the hydration layers has been taken as the position of the first minimum of the g i o n _ o ( r ) function. The ion-oxygen distances have been computed in 0.001 ps intervals

1

2

3

4

5

6

7

8

9

4

O-H

H-H

0

1

J J 2

.

.

3

.

.

4

.

5

.

.

6

.

7

8

9

rlA

Fig. 2. Radial distribution functions for the water sites in 0.555 molal TMACI solution.

o

over the whole simulation ran. The molecule was considered as belonging to the hydration shells if the ion-oxygen distance was less than the corresponding value of rm~.. The power spectra were averaged over several blocks of 20000 time steps (z = 5 ps). Only molecules staying in one of the subsystems within the whole time interval r have been included in the calculation of the power spectra. In Table 2 the frequencies of the OH stretching and HOH bending vibrations are summarized. The bending frequency is almost independent of the surrounding of the water molecule. The peak position of the HOH band noticed for pure water and the hydration shells and bulk water in TMAC1 and NaCI solutions are the same. The OH stretching frequency is sensitive to the surrounding of the water molecule and the differences between the frequencies obtained for pure water and the various subsystems in TMAC1 and NaCI solutions are remarkable. In order to dis-

E. H a w l i c k a , T. D l u g o b o r s k i / C h e m i c a l

I

+I +1 ~l" oo

Physics Letters 2 6 8

(1997) 325-330

329

cuss the influence of TMAC1 on the OH vibrations the frequency shift (A tO)subsystemhas been calculated as follows:

I (A

to)subsystem

=

&)subsystem - - t O g a s ,

(1)

~-H

Wsubsystem and tOgas denote the peak position of the OH band for the subsystem and for the 'gas phase' water, respectively. The frequency in the gas phase was obtained from a simulation neglecting all intermolecular interactions. The results are listed in Table 2 and are compared with the data for pure water and 1.1 molal NaCI solution [27]. In TMACI solution the OH frequencies obtained for each subsystem are shifted to lower wavenumbers. The red-shift for bulk water is the same, within statistical uncertainty, as that for pure water. This shows that TMAC1 does not affect the H-bonds in bulk solvent, neither their number nor their strength. The frequencies of the molecules in the hydration shells of C1- in solutions of NaC1 and TMAC1 show almost the same red-shifts. As compared with bulk water there is an additional red-shift by about - 2 5 _+ 5 cm -~. This indicates that the H-bonds formed between the anion and the water molecule are stronger. Compared with bulk water the molecules in the TMA hydration shells have OH frequencies shifted to lower wavenumbers by about - 2 0 +_ 5 cm -z, Because the average numbers of H-bonds in the hydration shells of TMA and bulk water are similar the additional red-shift indicates that the H-bonds are stronger. The strengthening of the H-bonds of the water molecules in the vicinity of cations has been noticed for Na +, but that effect is remarkably stronger and an additional red-shift of about - 9 0 cm-~ has been found [27]. The strengthening of the H-bonds in the field of cations has been reported from spectroscopic experiments as the 'cooperative effect' [28]. Though the cooperative effect for TMA cations is small, it shows, however, that this cation exhibits some 'ion properties' and they may induce strengthening of H-bonds. where

o _= -H +1

Z

I Z

o

¢-4 +1 +1

"-d ..Z

¢~

+1 --H

< I tt~ +1 +1 0

o

z

0

.~

.=_ -H +l

o

o I -H-H 00 c~ o o

e~

..~ 0

+i +l

-7

.~_°

~

t'N

+1 +1 ~ m

©

m +1 +1

I

+1 -H

+t +1

I +1 -H

.~_

~ =©

Acknowledgements

m

Authors thank the Polish State Committee for Scientific Research for the financial support (grants

330

E. Hawlicka, T. Dlugoborski / Chemical Physics Letters 268 (1997) 325-330

no. 3 T09 010 10) and Computer Center of Technical University in Lodz and Supercomputer Center Cyfronet in Krakow (grant no. CD/079/95) for computational facilities.

References [1] A. Pohorille, L.R. Pratt, J. Am. Chem. Soc. 112 (1990) 5066. [2] P.M. Wiggins, Microbiol. Rev. 54 (1990) 432. [3] J.Z. Turner, A.K. Soper, J.L. Finney, J. Chem. Phys. 102 (1995) 5438. [4] J.L. Barrat, I.R. McDonald, Mol. Phys. 70 (1990) 535. [5] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [6] O. Matsuoka, E. Clementi, M. Yoshimine, J. Chem. Phys. 64 (1976) 1351. [7] R.J. Bartlet, I. Shavitt, G.D. Purvis, J. Chem. Phys. 71 (1979) 281. [8] K. Toukan, A. Rahman, Phys. Rev. B. 31 (1985) 2634. [9] J. Anderson, J.J. Ullo, S. Yip, J. Chem. Phys. 87 (1987) 1726.

[10] S.B. Zhu, C.F. Wong, J. Chem. Phys. 99 (1993) 9047. [11] D.M. Ferguson, J. Comput. Chem. 16 (1995) 501. [12] P. Bopp, G. Jansco, K. Heinzinger, Chem. Phys. Lett. 98 (1983) 129.

[13] K. Heinzinger, in: Computer Modelling of Fluids, Polymers and Solids, C.R.A. Catlow (Ed.), Kluwer, Dordrecht, 1990, p. 357. [14] G. Palinkas, E. Hawlicka, K. Heinzinger, Chem. Phys. 158 (1991) 65. [15] E. Hawlicka, D. Swiatla-Wojcik, Chem. Phys. 195 (1995) 221. [16] E. Hawlicka, D. Swiatla-Wojcik, Comput. Chem. in press. [17] B.E. Conway, R.E. Verrall, J.E. Desnoyers, Trans. Faraday Soc. 62 (1966) 2738. [18] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, New York, 1987. [19] J.L. Finney, J. Turner, Physica B 156 and 157 (1989) 148. [20] J. Turner, A.K. Soper, J.L. Finney, Mol. Phys. 70 (1990) 679. [21] J. Turner, K. Soper, J. Chem. Phys. 101 (1994) 6116. [22] W.L. Jorgensen, J. Gao, J. Chem. Phys. 90 (1986) 2174. [23] A.K. Soper, Chem. Phys. 202 (1996) 295. [24] P.O. Eriksson, G. Lindblom, E.E. Bumell, G.J.T. Tiddy, J. Chem. Soc. Faraday Trans. 1 84 (1988) 3129. [25] I.I. Vaisman, F.K. Brown, A. Tropsha, J. Phys. Chem. 98 (1994) 5559. [26] R.L. Mancera, A.D. Buckingham, J. Phys. Chem. 99 (1995) 14632. [27] E. Hawlicka, D. Swiatla-Wojcik, Chem. Phys. in press. [28] H. Kleeberg, in: Interactions of Water and Ionic and Nonionic Hydrates, H. Kleeberg (Ed.), Springer Verlag, Berlin, 1987, p. 89.