**1. Introduction**

398 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Yoshikado, A.; Ohachi, T.; Taniguchi, I.; Onoda, Y.; Watanabe, M. & Fujiki, Y. (1982). AC

Ziman, J.M. (1960). *Electron and Phonons. The theory of Transport Phenomena in Solids,* pp. 161-

*State Ionics,* Vol.7, pp. 335-344

162, Clarendon Press, Oxford, UK

Ionic Conductivity of Hollandite Type Compounds from 100 Hz to 37.0 GHz. *Solid* 

Sodium, calcium and magnesium ions are essential for the biological activity of many polyelectrolytes. This activity depends on a condensation of the metal ions. There are several clues, which suggest that interactions between the polyelectrolyte and metal ions depend on a hydration of the ion. This accounts for the great interest in the hydration of metal ions, particularly in the systems containing hydrophobic groups. In aqueous solutions the hydration of Na+, Mg2+ and Ca2+ differs. The first hydration shells of Na+ and Mg2+ consists of six water molecules and have octahedral symmetry (Dietz et al, 1982; Hawlicka & Swiatla-Wojcik, 1995). The first shell of Ca2+ is larger; it contains eight or more water molecules and does not show any regularity (Owczarek et al., 2007). X-ray diffraction studies (Tamura et al., 1992, Megyes et al., 2004) have suggested that all these cations are sixcoordinated in methanolic solutions thus their shells are octahedral.

Various experimental techniques can be employed to gain insight into a coordination shell of the ion, but a lack of theory renders even a term 'preferential solvation' misleading. The concept of the preferential solvation has been introduced to explain non-linear changes of solution properties, but now this term is commonly used to emphasise a difference of the compositions of the coordination shell and the bulk solvent. Preferential solvation is usually expected if the ion interacts stronger with one of the solvent components. There are however experimental clues that the selective solvation might be due to a microheterogeneity of the binary solvent.

Methanol-water mixture is a suitable model to study structural aspects of solvation in binary systems, particularly when hydrophobic effects may occur. Both net components are highly associated liquids, but their hydrogen-bonded networks are inconsistent. Water molecules form a 3-dimensional, tetrahedrally coordinated structure, where cavities are filled with monomers (Soper & Phillips, 1986). Extension of the hydrogen bonds over 1 nm causes that liquid water, even at room temperature, behaves like a gel (Dore et al., 2000). Hydrogen-bonded molecules of methanol form zig-zag polymer chains (Narten & Habenschuss, 1984). Though in binary mixture the molecules of methanol and water may form a common hydrogen-bonded network the bulky methyl group causes that methanol molecule cannot simply replace the water molecule in the tetrahedral structure. In consequence the methanol-water mixture may become heterogeneous on the molecular level. Neutron diffraction (Dughan et al., 2004) and X-ray spectroscopy (Guo et al. 2004)

MD Simulation of the Ion Solvation in Methanol-Water Mixtures 401

Several potentials have been proposed to describe the interaction of the water molecules, but this molecule still remains 'a challenge to model, because it is polar, polarizable, has light H atoms and is flexible' (Heyes, 1998). Only a few of the models consider the water molecule as a flexible body and permit internal vibrations. Most of them are a modification of the simple point charge potential of Berendsen *(*Berendsen et al., 1981) They use a potential of the spectroscopic type (Zhu & Wong, 1993; Ferguson, 1995) to describe an intramolecular interaction. Unfortunately SPC model neglects non-Coulomb interactions of hydrogen

In simulations presented here the interactions of the water molecules has been described by the BJH potential (Bopp et al., 1983; Jancso & Bopp, 1983) which treats the intermolecular interactions of oxygens and hydrogens by mean of the central force model (CF) (Stillinger & Rahman, 1978) and uses a three body description of intermolecular interactions. This model has been frequently employed to simulate aqueous solutions of electrolytes, both in classical (Dietz et al., 1982; Jancso et al. 1985; Probst et al., 1985, Probst et al., 1991; Hawlicka & Swiatla-Wojcik, 1995, Lavenstein et al. 2000; Ibuki & Bopp, 2009) and in QM/MM MD (Tongraar & Rode, 2003, Rode et al., 2004, Öhrn & Karlström, 2004; Tongraar & Rode, 2005, Payaka et al., 2009; Tongraar et al., 2010) simulations. The BJH potential is appropriate to simulate the methanol-water mixtures, because it is fully consistent with the PHH flexible model (Palinkas et al. 1987) of the methanol molecule. The BJH and PHH, potentials reproduce properly the structure, energies and dynamic properties of the methanol-water mixtures (Palinkas et al., 1991a; Palinkas et al., 1991a, Hawlicka & Swiatla-Wojcik, 2000). An advantage of the flexible models is, that they permit a distortion of the solvent molecules from their equilibrium geometry. In consequence a molecular polarizability is incorporated. The BJH and PHH potentials consist of two parts, which describe the inter- and

intra inter (, ) () ( ) *Vr V V r*

The intermolecular parts are the sum of Coulombic and non-Coulombic terms. The Coulombic terms result from the partial charges of the interacting sites. In water molecule the partial charges are located on oxygen (-0.66 e) and hydrogen (+0.33 e) atoms. Methanol molecule consists of the charged oxygen (-0.6 e), hydroxyl hydrogen (+0.35 e) and the methyl group (+0.25 e), considered as the pseudo-atom. Non-Coulombic intermolecular O-O, O-H and H-H interactions of the water and methanol molecules are the same as in CF2 model for water (Stillinger & Rahman, 1978) The non-Coulomb interaction of the methyl group with the hydroxyl hydrogens has been neglected and that with oxygens and methyl

Intramolecular potentials for water and methanol are based on the formulation proposed by Carney et al. (Carney et al, 1976). They are expressed as power series of the internal

> 

Usually in classical MD simulations the ion potentials are represented by the Coulomb and the Lennard-Jones terms. These potentials overestimate, however, the number of the solvent

 *i i* 

groups has been represented by the Lennard-Jones potential (Jorgensen, 1981).

coordinates, i, 'stretch' and 'bend' and the three-body interactions are included:

*VL L L* ( )

*i i <sup>j</sup>*

 

 *<sup>i</sup> <sup>j</sup> ijk i <sup>j</sup> k ijkl i <sup>j</sup> k l* (2)

(1)

atoms, which seem to be important in hydrogen-bonded systems.

**2. Effective pair potentials** 

intramolecular interactions respectively:

have confirmed that supposition. Despite apparent miscibility of both components, methanol and water clusters are observed over whole concentration range. At particular concentration, near 25-27 mol% of methanol, where transport and thermodynamic properties exhibit extrema, water and methanol form separate, percolating structures (Dughan et al. 2004).

Several experimental techniques were employed to investigate solvation of ions in methanol-water mixtures. Results are, however, inconsistent and lead to contradict conclusions. Therefore both a preferential hydration (Convington & Dunn, 1989; Hawlicka, 1995), as well as a lack of preferences (Holtz et al., 1977) have been postulated for alkali and halide ions in methanol-water mixtures. Though X-ray and neutron scattering measurements should provide a direct insight into the ion coordination shell, their results cannot be decisive for methanol-water mixture, because distances between these ions and the oxygens of either water (Neilson & Enderby, 1979; Licheri et al, 1975) or methanol (Megyes, 2004) are almost the same. Moreover a direct correlation between the alkali earth cations and the methyl group is lacking (Radnai et al., 1995). In such case the scattering techniques are not enough sensitive to investigate the preferential solvation in methanolwater mixture. Thus a molecular dynamics simulation seems to be a useful tool to provide additional information concerning the structure of the ion shell.

A quality of the simulation results depends on the methods used to describe all interactions in the solutions. *Ab initio* quantum mechanics would be the most accurate method, but its application to systems containing ions and a few hundred water molecules could not be expected for the near future. Therefore QM/MM MD simulation seems to be an elegant approach for investigating the aqueous solutions of electrolytes. The ion and its nearest water molecules are treated quantum mechanically. Such approach includes many-body interactions between particles within the solvation shell. QM/MM MD simulations were carried out for aqueous solutions of various ions (Tongraar & Rode, 2003, Rode et al., 2004, Öhrn & Karlström, 2004; Tongraar & Rode, 2005, Payaka et al., 2009; Tongraar et al., 2010). Their results evidenced a significant role of the many-body interactions on structural and dynamical properties of the hydrated ions. These simulations concerned, however, the systems, which contain only one ion, either cation or anion, and a few water molecules. Thus this technique is useless for studies more concentrated solutions, where an association of the ions occurs. The formation of the various types of the ionic pairs, solvent separated, solvent shared and contact pairs, is frequently observed in binary solvents. This phenomenon may affect significantly the solvation of ions.

Classical MD simulations are usually carried out for NVE or NVT ensembles. The volume is fixed and it depends on a number of the particles, temperature, composition of the system. Particles are placed into a periodic cube. The size of the periodic box results from the experimental density of the simulated system. Initial coordinates of particles are frequently chosen from the crystal lattice (Heyes, 1998), however the random distribution of the particles in the cube is better, because it reduces a time of equilibration.

In classical simulations the interactions between molecules are represented by a sum of the pair potentials and many-body interactions are neglected. Usually the pair potential consists of Coulomb term, for which the Ewald summation is applied, and of short-range parts, for which shifted-force potential method (Allen & Tidesly, 1987) is used.
