1 Molecular Dynamics Simulations of Liquids with Ionic
Downloaded via 79.133.106.109 on July 12, 2018 at 21:06:38 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
Interactions K. HEINZINGER, W. O. RIEDE, L. SCHAEFER, and GY. I. SZÁSZ Max-Planck-Institut fuer Chemie (Otto Hahn-Institut), Mainz, Germany
Molecular dynamics (MD) simulations of liquids with ionic interaction have so far been performed for molten salts and aqueous electrolyte solutions. The characteristic problem for this kind of simulation are the far ranging Coulombic forces. The f i r s t preliminary MD calculations for molten salts have been reported by Woodcock in 1971 (1). The large amount of work published in the meantime has been reviewed by Sangster and Dixon (2). In the case of aqueous electrolyte solutions only preliminary results for various alkali halide solutions have been published so far (3). The more advanced state of the art for the molten salts leads to a concentration of the effort on the improvement of the algorithm for the integration of the equation of motion necessary to calculate dynamical properties with still higher accuracy. One example where very high accuracy is needed is the calculation of the isotope effect on the diffusion coefficients. In the section on molten salts below one way to improve the algorithm is d i s cussed and the progress is checked on the mass dependence of the diffusion coefficients in a KCl melt. Results with a certain degree of r e l i a b i l i t y from MD simulations of aqueous solutions reported up to now are restricted to structural properties of such solutions. In the section on aqueous solutions below very preliminary velocity autocorrelation functions are calculated from an improved simulation of a 0.55 molal NaCl solution. The problem connected with the stability of the system and the different cut-off parameters for ion-ion, ionwater and water-water interactions are discussed. Necessary steps in order to achieve quantitative results for various dynamical properties of aqueous electrolyte solutions are considered. Molten Salts The Hamilton differential equation system can be solved 0-8412-0463-2/78/47-086-001$07.00/0 © 1978 American Chemical Society Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
2
COMPUTER MODELING OF MATTER
numerically for a N particle system with a simple "leap-frog" algorithm as used by Verlet (4) or by a predictor-corrector a l gorithm. The most important criterion for the choice of the a l gorithm is the numerical s t a b i l i t y . The final decision follows from necessary accuracy with which for example the energy is preserved, which in turn is determined by the desired accuracy of the transport properties, say self diffusion coefficient. In the simulated system of molten KC1 i t has been investigated how the total energy, the total momentum, and the velocity autocorrelation functions varies with different algorithms. The MD calculation for the KCl-system was carried out at the melting point (1043K). The basic cube contains N=2-108 particles and the cube length S is than calculated to be 20.6A. For the pair potential the Born-Mayer-Huggins potential ip(r) = + er"
1
+ bexp(-Br) + Cr" + Dr' 6
(1)
8
is used with the parameters given by Tosi and Fumi (!5), listed in Table I. Table I Parameters for the Born-Mayer-Huggins Potential b
+
K
+
-
K
+
- CI"
CI"'-
cr
D
C
10~ erg
A-1
1991.67768
2.967
1224.32459
2.967
-
4107.95500
2.967
- 145.5
12
K
B
10" erg A 12
-
6
10" erg A8 12
24.3
-
24.0
48.0
-
73.0
- 250.0
Three different algorithms were investigated. In the f i r s t version (I) the Coulombic energies and forces were evaluated with the use of the erfc part of the Ewald method (6) only. In the other two versions (II and III) the f u l l Ewald method was employed. In a l l versions the separation parameter T), the number of nearest neighbours n~ , and the maximum value of the reciprocal lattice vectors h] (not occuring in version I) were chosen to be: i\ = 0.175,
n^ = 6,
h
1
= 8.
The algorithm which was used in version I is a so called "leapfrog" algorithm: Ii(n+1) = r.(n-1) + 2At v.(n-1) + (At /m.) JF.(n) 2
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
(2)
1.
HEINZINGER ET AL.
MD Simulations of Liquids with Ionic Interactions
3
v.(n) = v^n-1) +[v.(n-1) - v.(n-2)] + (A^/m^F.fn) Here r.-j( )> v.i( )> n
n
a n c l
Lj( ) n
a
r
e
t
h
e
- ^(n-1)]
(3)
position, velocity and force
of particle i at time t = nAt respectively and m-j are the masses -14 of the ions. The time step length was chosen to be 0.5-10 sec. The behavior of the total energy for version I is shown in F i gure 1 where (AE/E) =6 - 10-3. A change from the "leap-frog" algorithm to a predictorcorrector algorithm has been made in version II. For the position vectors one can write: r?(n+1) = r?(n-1) + 2 At v.(n-1) + (At /m.)[ F.(n) + F-(n-l)] 2
(4)
r?(n+1) = r?(n-1) + 2 At v.(n-1) + (At /4m.)[ F.(n+1) + 4F.(n) + 3F.(n-1)] 2
(5)
Here the indices p and c indicate predicted and corrected values respectively. The formula for the velocities has been changed slightly compared to the formula (3): ^ ( n ) = v.(n-1) + (At/2m.)[F (n) + F.(n-1)] i
(6)
With this version II, two simulations with different time steps were carried out starting from the same configuration (0.5 10"1^s and 0.25 10~14s) over a time interval of 1.2 ps. The total energies for both runs are shown in Figure 1 . The value (AE/E)Q ^ has decreased compared to the one of version I by one order of magnitude. Moreover AE/E decreased further by a factor 0.5 when the time step was shortened to 0.25.10"14 . In version III the positions and velocities have been treated with the same predictor-corrector algorithm. This means that F-j ( t ) , the derivative of the forces with respect to time, appears in the formula for the velocities. s
rP(n+1) = r?(n-1) + 2 At v?(n-1) + (2At /3 ) [2F (n) + F.(n-1)] 2
mi
i
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
(7)
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
J
1 01
Etotal
1 •
1 05
.
.
3
.
.
, 1 psec
.
.
•
3
Figure 1. Total energy of version I ('upper curve) in units of 10 J/mol and of version II flower curve) with a energy scale enlarged by a factor of 4 A. The length marked by (jimm) in the upper scale corresponds to the interval of the lower scale. The time scale starts after an equilibration run over 7.2 psec from a roughly equilibrated starting configuration.
627 6
633-
1.
HEINZINGER ET AL.
MD Simulations of Liquids with Ionic Interactions
r?(n+1) = r + 1 (-w)(-1) (-2)
V
d
+
2
w w
=(
+
>
2
0=1
H )
J ^
e
I
a
+
B
/ d
a
a (-a)
d
+
V (r) = + 2 / r ++
e
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
B
16
COMPUTER MODELING OF MATTER
The switching function, S ^ r ) has been introduced to reduce unr e a l i s t i c Coulomb forces between very close water molecules ( 1 2 ) . e is the elementary charge and q = 0.23 e is the charge in the ST2 water model, d and r denote distances between point charges and LJ centres, respectively. The choice of a and 6, odd for positive and even for negative charge yields the correct sign. The LJ parameters for water-water interaction are taken from the ST2 model, and for cation-cation from the isoelectronic noble gases (Y3). The anion-anion parameters are derived from the cation ones on the basis of the Pauling radii as shown in detail in a previous paper (14). The LJ parameters between different part i c l e s are calculated by application of Kong's combination rules (J5J. A l l e and o are given in Table III. Table
III
Lennard - Jones parameters for the various pair potentials in a NaCl solution is given in A and in units of 10"*6 erg. H0
Na
3.10
2.92
4.02
54.84
30.83
2.73
3.87
5937
28.32
-
27.87
2
H0 2
Na cr
+
52.61
-
+
CI"
4.86
The classical equations of motion are integrated in time steps of At = 1.09-10-16 (22).In order to take account of the different strengths of the various interactions and to keep the computer time in reasonable limits, different cut off parameters are chosen for ion-ion, ion-water and water-water interactions. The Ewald summation (6), which gives the correct Coulombic energy for an unlimited number of image ions, is used for ion-ion interactions. As the ions residing in the f i r s t and second neighbour image boxes are included in the direct calculation of the so called error function part, a suitable choice of the separation parameter allows to neglect the Fourier part of the Ewald sum. The ion-water interaction is cut off at 9.1 A, half the sidelength of the basic periodic box. The water-water interaction is cut off at 7.1 A which means that only f i r s t and second nearest neighbour molecules are included in the calculation. The structural properties of the NaCl solution calculated from this simulation are not significantly different from preceding ones. The characteristic distances an heights of the ionwater and water-water radial pair correlation function remain a l most unchanged whether a temperature control mechanism is built s
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
1.
HEINZINGER ET AL.
MD Simulations of Liquids with Ionic Interactions 17
in or not. The radial pair correlation functions from the simulation of a 2.2 molal NaCl solution (200 water molecules and 8 ions of each kind) are compared in detail with X-ray investigations in a previous paper (14). The problem connected with the determination of hydration numbers have also been discussed previously (16). In addition i t has been concluded from the MD s i mulations of alkali halide solutions that in the f i r s t hydration shells a lone pair orbital of the water molecules is directed towards the cation and a linear hydrogen bond is formed with the anion (1?). This result is also in agreement with conclusions from experimental and other theoretical investigations. The real aim of MD simulations i s , of course, the calculation of the dynamical properties of the liquid. The velocity autocorrelation function have been calculated from this 10.000 time step run. In a l l cases the functions have been averaged over 360 starting vectors. The s t a t i s t i c a l quality is much better, of course, in the case of the water molecules where the ensemble average extends over 200 particles compared with 2 in the case of the ions. In Figure 8 the angular (solid line) and the translational (dashed line) velocity autocorrelation functions of the water molecules are shown. The correlation time for the angular velocity is as expected much shorter than for the translational velocity. In Figure 9 the velocity autocorrelation functions for the cation (solid line) and the anion (dashed line) are drawn. If these autocorrelation functions are integrated, in spite of a l l s t a t i s t i c a l uncertainties i t is found that the corresponding diffusion coefficients are of the same order of magnitude as is known or estimated from experimental investigations. Therefore i t is j u s t i f i e d to continue the investigation along this line. It is obvious that much longer simulation times are required to make quantitative statements about the dynamical properties. But not only an increase in simulation time is necessary, also the simulation i t s e l f needs further improvement. The system simulated corresponds to a microcanonical ensemble and consequently i t should have a constant total energy after an equilibration time of suitable length. In Figure 10 the total energy (solid line) and the potential energy (dashed line) for the system are shown over the 10.000 time steps run. I can be seen that the increase of the total energy is equally distributed between kinetic and potential energy. The unavoidable cumulative errors due to the algorithm, which were discussed above in detail in the case of the molten s a l t s , are not the main source of energy increase in the simulation of aqueous solutions. Responsible for the problems here are mainly the cut-off r a d i i . Figure 11 shows the variations of the total energy of the basic box E during a simulation run of 270 & t, which is essentially caused by the cut-off radii R and Ri in the evaluation of the potential energy E t ( t ) . Let us consider the expression: w
p o
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
18
COMPUTER MODELING OF MATTER
Figure 9. Normalized velocity autocorrelation functions qp(t) for cations (solid) and anions (dashed in a 0.55m NaCl solution)
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
1. HEINZINGER ET AL.
M D Simulations of Liquids with Ionic Interactions 19
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
20
COMPUTER MODELING OF MATTER
-232 AE
-233
50
100
150
200
250 At
Figure 11. Variation (cut-off jumps) of the total energy of the basic box in units of 10 erg for 200 water molecules and eight Ntf CZ. The time step is 2.18 • 16 sec and AE/E = 3 • 10 . 12
16
4
s
Lykos; Computer Modeling of Matter ACS Symposium Series; American Chemical Society: Washington, DC, 1978.
1.
E
MD Simulations of Liquids with Ionic Interactions 21
HEINZINGER ET AL.
pot
( t +
A t
>- pot E
•
tV (r (t^t)-V (r (t))}
( t ) =
i
1 j
1 j
i j
^(r^tt^t))-
1 j
^(r^t))
Z
(28)
where the summation extends over i and j such that I
r^.JtJ^R and ^ .(t+At) to S to
C to ο
H2
^
3
II .2 to
5
"•*» c Ο
3
^
II Î3 i
1 ο
