Calculation of Coulombic Forces

For the simulations of ioics, the treatment of the Coulombic force of long range is included. Traditionally, Ewald methods [7] and/or related methods such as particle mesh Ewald [8] are used for long- range Coulombic force. Multipolar expansion [20] is also applicable even when the system contains defects. Recently, Wolf method [21, 22] is also used in ionic systems in which the Coulomb force is truncated at a certain length because the force was considered not necessarily long scale due to the screening by other ions. The method is useful to save the cost of the calculation in several situations in ionic systems [23]. Charge neutrality of the system within a basic cell of MD, which is usually imposed to the system, allows to use the Wolf method in disordered system [23]; however, the charge neutrality imposed to the cell (equal number of cation and anion in the case of monovalent systems.) itself is one of the causes of the finite- size effects of the MD simulations, since it causes at least some structural constraints. Therefore, a careful choice of the system size is also required. The check of the applicability of each method in each case seems to be the responsibility of researchers. In the case of molten salt and/or ionic liquids as well as network systems and more complex systems, the Coulombic force is not necessarily screened within a short length scale and such length scale is a strong function of temperature and/or compositions as shown later. In the following sections, the existence of several length scales in some ionic systems including IL will be shown.

In nanostructured materials without PBC, the Ewald method and related one assuming periodicity cannot be used. A careful choice of the method of calculation of the Coulombic force is necessary.

Treatment and Tools of Nanostructured Materials in MD Simulations

MD simulations often treat the nanoscaled systems with and without periodic boundary conditions, and therefore general structural analyses mentioned in the following sections can be applied to nanostructured materials. Tools explained here are pair correlation functions, running coordination numbers, from which one can determine several length scales and fractal dimensions. Fractal dimension analyses [24, 25] are expanded to multifractal structures [26-28] as well as those for dynamics. Analyses of network structures, by using Q_n distributions [29-31], cluster analyses using a dendrogram [32] are also introduced and some of them will be applied to data obtained by MD simulations in Chapter 9.

In nanoparticles and/or porous systems, the roles and characteristics of the surface area are often targets of debates. Porosity and connectivity of pores are also used to characterize porous materials.

Some methods, calculation of the accessible surface area and geometric pore size distribution, analysis of the structural connectivity and percolation analysis of the porous space are explained in Ref. [33].

Results for the relations among porosity, density, and dynamics in porous lithium silicate systems can be found in Chapters 6 and 7.

Some of these values have correlated with each other. For example, as already discussed in Chapter 2, porosity is a function of the fractal dimension of holes. Since the porosity is a function of density (see Chapter 6), it means that the fractal dimension of the hole is a function of the density as schematically shown in Fig. 2.12 in Chapter 2. The similar relations should hold among density, packing fraction and fractal dimension of the structure, at least for the monofractal case.

Structural Information Obtained by Molecular Dynamics Simulations—With Examples of NaCl,Silicates, and Ionic Liquids

Pair Correlation Functions

Liquids have homogeneous structure at longer length scales and characteristics of them can be well represented by the function based on the radial distance r.

Pair correlation function [1, 2] is a fundamental function of the theory of liquids [34, 35] and can also be used for characterizing supercooled liquids and glasses.

From MD simulation, the pair correlation function among different kind of particles i and j is obtained by counting the number of j particles within the shell with width r at distance rfrom particle пДг- Ar/2, r + Ar/2).

where V is a volume of the system, /V, and Nj are the numbers of the species i and j and the term 4лг² is for the surface area of the sphere.

When the species i and j are the same one, the function is

where now n,(r - Ar/2, r + Ar/2) is the number of i particles within the shell with width Ar at distance r. The (N_t - 1) term is the number of surrounding particles, which does not include the central particle /.

Practically, the function can be calculated during the MD run and it is useful to check the status of the calculation.

Pair correlation functions can be calculated among ions, atoms, molecules. For example, in the case of ILs with internal structures, the function g(r) for each pair can be determined for the center of mass position (or center of charges) of each ion. These functions become 1 when the mean density is attained.

If one calculates the function during MD runs, distances among i and j particles will randomly appear, and this may prevent the effective coding for parallel computations. A coding with random access may be better to put it out of the loop of the calculation of force, if necessary.

Figure 4.1 An example of the correlation function, g(r). Data are taken from the result of molecular dynamics simulation for molten NaCl at 1500 К using the same potential parameters as in Chapters 8 and 9. Dashed (red) curve: Na-Cl pair, Dotted (black) curve: Na-Na, Solid (blue) curve: Cl-Cl pair.

The structure thus examined can be used to obtain the more coarse-grained model for the simulations with keeping the characteristics of the structures in the atomistic-level. Of course, it is also possible to calculate £f(r) for each atomic species or for partial structures if necessary. Starting from pioneer works by Fumi and Tosi [36, 37], nowadays many kinds of potential parameters are available for alkali halides including polarizable ion models [38]. Researchers can select suitable ones by the required quality, calculation cost, etc., to satisfy the objective of each work. In Fig. 4.1, an example of pair correlation functions for molten NaCl obtained by classical MD simulations is shown.

Here we used the potential parameters in [39], for which NaCl-water interactions in Chapters 8 and 9 are calculated so that results can be compared with colloidal solutions in these chapters. As usual, molten NaCl, Na-Cl pair shows a sharp peak with the shortest distance, while Na-Na and Cl-Cl pairs show quite similar functions to compensate for the charges in Na-Cl pairs.

Running Coordination Number and Fractal Dimension Analysis

The running coordination number N[r) is related to the pair correlation function g(r) as follows when / and j are different species [2, 3].

Here, n is chosen such that nAr becomes r, and < > indicates an ensemble average.

The coordination number of j atoms around / atom can be obtained as a function of the distance r and the first coordination shell is often defined by the value at the position of the first minimum of the pair correlation function of the i-j pair.

From the running coordination number, one can determine the fractal dimension, d_{, of the structure as follows:

If the system consists of several domains or hierarchy structures, more than one fractal dimension will be obtained as shown for several cases of gels in Chapter 3. The coexistence of several exponents can also be analyzed by the multifractal analyses as shown later.

The fractal dimension of the structure affects the motions on this structure, and it can be determined from the trajectories of any units of structures such as atoms, ions, and molecules.

Charge Radial Distribution Function

An alternative arrangement of cation and anion in the ionic system tends to neutralize structure at a certain length scale. This behavior is regarded as screening in the strongly coupled system such as NaCl [40, 41] and charge radial distribution in the following form is suggested.

where A/r is the amplitude, f is the phase shift, d is the period of the oscillations the Я corresponds to the Debye length.

In ionic structures, and one can consider the charge distribution function Q(r) [41-43] by assuming that the charge is simply on the center of mass (or charge) position,

For monovalent systems such as NaCl with formal charges, e (a weight by a unit charge) is 1. The function is related to layers of charges (charge density wave, CDW). In this function, the mean value for the structure around cation and that around anion was taken. This is reasonable for example in molten NaCl, for which g(r) for Na-Na pair is comparable to that for the Cl-Cl pair. However, if such structures are modified by other components, sometimes each central ionic species need to be distinguished. This function decays to 0 when the charge neutralization is attained at a certain distance. Function concerned with the density wave (DW) regardless of charge distribution can be considered even for the ionic system. The density distribution function G(r) is defined by the following equations.

Typical length scales for these functions are informative to learn the difference of the main factors controlling them. This function becomes 1 when the value becomes the mean density. If we compared the Q(r) and G(r) of the molten NaCl, as in Fig. 4.2, one can see the CDW continues longer than that for the

DW. The first peak in CDW is affected by the DW and not necessarily determined by charge ordering alone.

Figure 4.2 (a) Examples of charge distribution function,

Determination of Debye Length from Charge Density Wave (CDW)

If the plots of ln|Q(r)r| against r have a straight line when the maxima of the peaks (envelopes) are connected, the characteristic length Aq, which corresponds to the Debye length for the screening of the Coulombic term in the simple dilute ionic systems, can be determined from the slope, -1/Aq. Such analyses have been used to examine the length scales of ILs near the glass transition temperature [43] to show the contribution of different length scales. The method was applied to learn the length scale associated with the structure of salts in water (and water in salt) with different concentrations [44]. Figure 4.3 is an example of the function in the molten NaCl, where Coulombic interaction was calculated by the Ewald method. In this case, oscillation with the comparable widths and intervals continues for a long length scale. Caution should be paid for the physical meaning of the Debye length, since the values in dense ionic systems are not necessarily the same one as in a dilute system. The interaction observed is a renormalized one by

Figure 4.3 An example of the plot of ln|r(?(r)| against r. Data are for molten NaCl at 1500 K. From the slope connecting the maxima of peaks (envelopes), one can determine a Debye length, although the physical meaning of the dense ionic systems is not necessarily the same as in that in dilute solution.

the interactions from other ions. The Debye length determined from the slope is found to be 5.94 A and the value in the form of Eq. 4.3 means the length at which the |r(?(r)| value is reduced to 1/e (~0.36788) times of the initial value. Of course, it does not mean the correlation is lost completely at this length. Correlation in Fig. 4.3 continued to cover whole observed regions ~24 A or longer. Since the function is rapidly decaying, good statistics are required to discuss this kind of length scales.

Screening Length in Interfacial Systems

Screening at ionic surfaces is an interesting problem concerning with this structure of the CDW. Desai [40] has examined the system of NaCl confined between two slabs (neutral and charged ones). The function Q[z) = p(z)⁺ + p(z)~, where z is a distance measured from the surface and following one-dimensional form is successfully applied for both octopolar (111) and dipolar (110) surfaces.

where A is amplitude, Aq is the characteristic length corresponds to the Debye length. Here d and ф are the period of oscillation and phase shift, respectively.

Large effects of surfaces are reported, and their results suggest that the role of the apparatus used in the experiments cannot be neglected.