Molecular dynamics simulation of helium diffusion in

Journal of Non-Crystalline Solids 443 (2016) 47–53

Contents lists available at ScienceDirect

Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/locate/jnoncrysol

Molecular dynamics simulation of helium diffusion in vitreous silica Sergey V. Kuhtetskiy ⁎, Elena V. Fomenko, Alexander G. Anshits Institute of Chemistry and Chemical Technology, Siberian Branch of the Russian Academy of Sciences, Akademgorodok 50/24, Krasnoyarsk, 660036, Russia

a r t i c l e

i n f o

Article history: Received 29 February 2016 Received in revised form 6 April 2016 Accepted 7 April 2016 Available online xxxx Keywords: Vitreous silica Helium diffusion Molecular dynamics simulation Empirical potentials

a b s t r a c t The helium diffusion through the non-porous vitreous silica was simulated by the molecular dynamics (MD) method. Truncated Morse potential was used to create models of vitreous silica. It was found that the quality of the simulation has been determined by the two main factors: the correct model of matrix in terms of the interstitial space geometry, and the quality of the gas-matrix interaction potential parameterization. The density of silica glass was in a good accordance with the experimental value 2.2 g/cm3. The helium-silica interaction potential's parameterization has been done as close as possible to the conditions of a real migration of helium through the doorways connecting adjacent solubility sites. Taking into account of these two factors has allowed simulating the diffusion of helium through vitreous silica by molecular dynamics simulation quantitatively and calculating the diffusion coefficients which are close to the experimental values available in the literature. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Non-porous silicate glasses are considered as promising systems for the creation of highly selective gas separation membranes. Their gas mixture's separation factors are: α(Не/СН4) = 106, α(Не/Ne) = 102, and α(Не/Н2) = 10 at 400 °С [1,2]. Computer simulation of the gas diffusion predicts the membrane's diffusion properties according to its structure. It is useful to design the new highly selective membrane materials. The quality of the simulation depends on the matrix and gas-matrix interaction models. First of all due to geometry of the interstitial void space of matrix, and height of the potential barriers of the doorways connecting adjacent solubility sites where the atoms of a gas are dissolved. Unfortunately, these characteristics of the materials are quite difficult for direct experimental study of modern structural techniques. For example, even for a pure vitreous silica which structure is studied more than eight decades starting from pioneering works of Zachariasen [3], middle-order structure (from a few angstroms to nanometers) is still discussed [4,5]. A number of characteristics of the interstitial void space, which are important for practice, in particular, the size distribution function of solubility sites, were determined using the mechanical model of gas solution proposed by Shackelford and colleagues [6]. This model underlay the statistical theory of monatomic gas diffusion in vitreous silica [7], which showed good agreement with experimental data [8]. Unfortunately, the theory contains a number of parameters relating to the structure of a particular material (jump distance, vibration frequency in a doorway and activation energy for diffusion), which were determined

⁎ Corresponding author. E-mail address: [email protected] (S.V. Kuhtetskiy).

http://dx.doi.org/10.1016/j.jnoncrysol.2016.04.013 0022-3093/© 2016 Elsevier B.V. All rights reserved.

by fitting to the experimental data. Therefore, this theory does not predict reliable diffusion properties of a newly developed material. Model can be constructed by MD-simulation method that does not use parameters which are associated with the matrix structure. This method gives an opportunity to simulate the main stages of the membrane's material design: structure, synthesis of the samples, diffusion properties. However, there are few papers dealt with this method of the gas diffusion modeling through non-porous matrix in the literature. In particular, molecular dynamics simulation has been used to model diffusion of helium through the alpha-quartz and vitreous silica [9]. Unfortunately, a small size of the samples (144 molecules of SiO2) did not allow the authors to obtain quantitative data and estimate the advantage of the molecular dynamics for simulation of gas diffusion through glasses. This paper is devoted to molecular dynamics simulation of helium diffusion into vitreous silica, including (a) preparation of the matrix model, which correctly describes the geometry of the interstitial void space, (b) helium-matrix interaction, (c) the calculation of numerical values of the diffusion coefficient of helium over a wide temperature range, and (d) their comparison with the experimental values available in the literature.

2. Simulation details Molecular dynamics simulation model includes: model of vitreous silica (matrix), model of helium-matrix interaction, the preparation of matrix samples with embedded helium atoms at the given pressure and temperature, and the simulation of the migration of helium atoms in vitreous silica, followed by the calculation of their mean square displacement and diffusion coefficient.

48

S.V. Kuhtetskiy et al. / Journal of Non-Crystalline Solids 443 (2016) 47–53

2.1. Model of vitreous silica MD-simulate of vitreous silica preparation includes the real manufacturing process, as a rule: silica melting and cooling to room temperature (a long time dynamics). The correct representation of the geometry of the interstitial void space (correct statistics of solubility sites and doorways between them) requires quite a large volume of the system, at least a thousand or more molecules of SiO2. Quantum chemistry methods are too expensive for such systems. The usual practice is the using of molecular dynamics simulation for this. Interatomic interaction is described by the empirical potentials. There are many variants of using such technique, which differ by the forms of empirical potentials, methods of their parameterization, and the quenching modes [10–20]. The most part of these models does not deal with the problem of dissolution and diffusion of gases in glass, unfortunately. Problems of correct simulation of the interstitial void space in glasses have received a little attention up to now. In particular, many well-established empirical potentials give significantly overestimated (sometimes, by several tens of percent) densities of vitreous silica compared to the experimental values [21–25]. Due to high sensitivity of the gas diffusion of the interstitial void space's geometry, we cannot talk about a quantitative diffusion simulation based on such models. The effect of the empirical potentials on the density of vitreous silica prepared by molecular dynamics simulation has been studied in [26,27]. Samples of vitreous silica were prepared by the simulation of the silica melt quenching at the constant pressure. Covalent component of Si\\O bond was specified by Morse potential multiplied by a weighting function, which allows cutting the tail of original Morse potential smoothly. It has been shown [27] that the long-range cutting of the potential (at distances ~3 Å) does not effect on the structural elements of glass, but significantly effects on their packaging. A good correspondence between the density of the resulting vitreous silica and the cutting radius in the range of 2.5–5 Å has been found. The value of the glass density is closed to the experimental value 2.2 g/cm3 when the cutting radius is neared to 3 Å. Thus, this model can be used to molecular dynamics simulation of gas diffusion therein. In this paper the model samples of vitreous silica was prepared according to the [26]. Interatomic potential is given as: vðrÞ ¼

h i2 Z Si Z O 1 þ D 1−e−aðr−r0 Þ −1 0 r eβðr−r1 Þ þ 1

ð1Þ

where the ionic component of the ionic–covalent bond is represented by the Coulomb term, and the covalent component is described by the Morse potential multiplied by the weighting function g(r), which determines the “soft” cut-off of the Morse potential [13,26,27], ZSi = +2.4e and ZO = −1.2e are the charges of the silicon and oxygen ions, respectively. The values of the parameters D, α, and r0 were taken from [28], and parameters β and r1 – from [26]. For r1 ~ 2.95 Å, the density of the model vitreous silica is close to the experimental value. All parameters of potential (2) are presented in Table 1. The samples of vitreous silica, dedicated to simulate the diffusion of helium therein, were prepared as follows. 1000 silicon atoms and 2000

Table 1 Parameters of the potentials. Parameter

D, eV α, 1/Å r0, Å β, 1/Å r1, Å

Bond Si\ \O

O\ \O

Si\ \Si

0.340554 2.006700 2.100000 10.0 2.95

0.042395 1.379316 3.618701 – –

0 0 0 – –

oxygen atoms in the cubic unit cell with periodic boundary conditions was nested. The initial coordinates of the atoms were corresponded to the structure of cubic cristobalite. Velocities of particles was set so as to correspond to Maxwell distribution with T0 = 300 K. The initial parameters of the thermostat (T0) and the barostat (atmospheric pressure P0) were set, and molecular dynamics simulation was started. After a short relaxation period ~ 0.1 ns at T0 = 300 K, the system was heated up to the T1 = 6000 K using the linear increase of the thermostat temperature. The duration of the heating was 0.2 ns. After reaching the temperature T1 and a short period of relaxation (0.1 ns) at this temperature, the system was linearly cooled to room temperature T0 for 10 ns. Normal pressure (0.1 MPa) was maintained in the system during the entire process. The time step was equal to 2 fs. The dynamics of the atoms was calculated with the HOOMD-blue software package [29–32], which effectively implements parallel computations on graphics processors using the CUDA technology [33,34]. The electrostatic component was calculated by the particle–particle/ particle–mesh (PPPM) method implemented in HOOMD-blue. The Morse potential was converted to the tabulated form (for details, see the documentation for HOOMD-blue [29]). The simulation speed (for 3000 atoms, NPT ensemble) with the use of the ASUS GeForce GTX 780Ti DirectCU II graphics cards was approximately 3–4 ns/h. 2.2. Helium-matrix interaction Potential of interaction of helium with matrix atoms is important part of quantitative molecular dynamics simulation of diffusion process. There are some variants based on Lennard-Jones and Born-Mayer potential in literature [35–38], but their parameterization was done in configurations too far from diffusion problems ones. The acting force on an atom of helium placed within interstitial void consists of two parts: a short-range repulsive force owing to the electron shells of the gas atoms and matrix overlapping, and attractive dispersion and polarization ones. The limiting stage of diffusion of helium atoms in the glass is to overcome the potential barrier in the doorway connecting adjacent solubility sites. These doorways are Si\\O-rings with the ring size distribution maximum at about 6-member ring [39] where repulsion prevails [40]. Therefore only repulsive pair empirical Born-Mayer potential was used in this work: U ðr ÞHe−o ¼ A e−Br

ð2Þ

where r is the distance between the helium atom and the oxygen ion, A and B – constants were calculated as follows. Flat siloxane rings can be considered as an idealized model of doorway connecting adjacent solubility sites in silica glass. The potential energy of helium atom placed in ring's center has been calculated using density functional theory for several flat siloxane rings (from 4 to 8 Si\\O members) [40]. These data can be used to parameterization of Eq. (2) as follows. Electric field in center of ring was equal 0. Therefore, as a first approximation, one can assume that the contribution of polarization forces to helium energy is low and repulsive forces between helium and oxygen atoms and silicon ions are dominated. Moreover, if one take into account the fast decay of potential (2) with distance and the fact that oxygen ions are located more closely to the center of ring then silicon ion, one can assume that the oxygen ions give the main contribution to potential energy of helium atom centered in the ring. With these assumptions, one can plot energy of helium versus helium-oxygen distance and parametrizate potential (Eq. (2)) using fitting. 2.3. Preparation of matrix samples with embedded helium atoms For the helium diffusion calculation, it was needed to place a certain amount of helium atoms into solubility sites. This procedure simulates


the helium atoms dissolution into vitreous silica. The ratio of solubility sites occupied by the helium atoms, with respect to all available solubility sites can be evaluated as follows. The concentration of all available solubility sites for helium in vitreous silica is Cs ~ 2.3 ⋅ 1021 cm−3 [41]. The solubility of helium in vitreous silica at atmospheric pressure and room temperature is γ = Ci/Cg = 0.025, where Ci is the concentration of helium atoms inside the vitreous silica, and Cg is the helium concentration in the gas phase [42]. The helium concentration in the gas phase is Cg ~ 2.7 ⋅ 1019 cm−3 at standard temperature and pressure. Hence the concentration of helium atoms in the matrix is Ci ~ 6.8⋅ 1017 cm−3. Let us assume that any solubility site being occupied contains not more than one helium atom. This assumption will give an upper estimate of the number of occupied sites. As a result the number of occupied solubility sites in relation to the total number is KHe = Ci/Cs ~ 3 ⋅ 10−4. This ratio is very small. Hence the dynamics of each dissolved helium atom can be calculated independently of the other helium atoms. Moreover, the volume of the unit cell containing 1000 atoms under normal conditions (density 2.2 g/cm3 ) is approximately

49

equal to V ~ 4.5 ⋅ 10− 20 cm3. Therefore the average number of helium atoms in the cell is 0.03 only. In principle the above estimates make it possible to reduce the diffusion part of the problem to simulation of only single helium atom migrated through unit cell. However, aiming for quantitative modeling of diffusion one should use a large ensemble of helium atoms noninteracting among themselves for better statistics. Fig. 1 shows an example of the trajectory of the helium atom in a matrix with its migration through the 6-membered ring. The example was obtained in this paper. Visualization of the silicon-oxygen network of matrix and helium path was made with VMD package [43,44]. Fig. 1 demonstrates the character of the motion of a helium atom in vitreous silica, which consists in a sufficiently long duration localized motion (vibration) inside the interstitial voids with fairly rare transitions to adjacent solubility sites through the doorways formed by rings of the silicon–oxygen tetrahedra. Because of this nature of motion, the mean square displacement required for the diffusion coefficient evaluating, calculated from the trajectory of a single atom will strongly fluctuate. Therefore, for correct

Fig. 1. Stereo pair of the trajectory of a helium atom during the migration through a six-membered ring of vitreous silica; the change in the color of the trajectory corresponds to the change in the time: from red at the initial time through white to blue at the final time. It is a cross stereo pair; i.e., the right eye looks at the left image, and vice versa (top). Dependence of the coordinates of a helium atom on the migration time (bottom). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

50


calculation the mean square displacement should to introduce a sufficiently large ensemble of helium atoms (usually, 1000 atoms) with the assumption that helium atoms cannot interact among themselves but cat interact with matrix atoms. The preparation of the matrix and the embedding of helium atoms into it were performed as follows. First of all, the matrix was “heated” up to a given temperature. After a short period of relaxation, the positions of the matrix atoms were fixed. After that, helium atoms were embedded in the matrix one after another according to the following procedure. Random velocity (according to the Maxwell distribution at a given temperature) was set for each helium atom. Then random position of the atom inside unit cell was set. Potential energy of helium atom placed at the position was calculated using the Born–Mayer potential (Eq. (2)) and the periodic boundary conditions. If the kinetic energy of the helium atom exceeded the potential energy, the atom was placed at the position. Otherwise, the next attempt was made with new random values of the velocity and coordinates.

2.4. Simulation of helium diffusion After the preparation of the system according to the procedure described in the previous section, the calculation of the helium and matrix atoms moving was started. The molecular dynamics package HOOMDblue [29–32] was used. To avoid a non-physical influence of the pressure of the helium atoms ensemble on the size of the unit cell, the dynamics of the entire system was calculated as NVT ensemble. After an initial period of relaxation, the current coordinates of the helium atoms were stored. The further simulation was accompanied by periodic calculations of the mean square displacement of the ensemble of helium atoms with respect to these coordinates. Some examples of time dependences of the mean square displacement of helium atoms for several temperatures of the system are shown in Fig. 2. It is seen that, owing to the use of a sufficiently large ensemble, fluctuations of the mean square displacement of helium atoms are relatively small. After reaching the linear increase in the mean square displacement Δr2 as a function of the time t, the diffusion coefficients D was calculated according to Einstein formula:

D¼

1 Δr2 : 6 t

ð3Þ

Table 2 Density of vitreous silica samples (g/cm3). Sample no. 1

2

3

4

5

6

2.2028

2.2041

2.1972

2.2021

2.2012

2.2035

Average

Stdev

[45]

2.2018

0.0025

2.2

3. Results 3.1. Characteristics of the vitreous silica samples Six samples of vitreous silica using the technique described in Subsection 2.1 were prepared for molecular dynamics simulation of helium diffusion. The main characteristics of these samples are presented in Tables 2–5. The right side of each table contains the values averaged over six samples, the standard deviations, and the corresponding experimental values taken from the literature. The density calculated for each sample at atmospheric pressure and room temperature are given in Table 2. The positions of the maxima of the first three peaks in each pair radial distribution function for the same samples at room temperature and atmospheric pressure are presented in Table 3. The distributions were calculated taking into account the periodic boundary conditions. The positions of the maxima of the distribution functions of the angles between the Si\\O\\Si and O\\Si\\O bonds are presented in Table 4. The values were calculated by the nonlinear regression (fitting) method with a sinusoidal function for the distributions of the Si\\O\\Si bonds and a Gaussian function for the O\\Si\\O bonds. The ring statistics is one of the most important factors characterizing the spatial network structure of silicon–oxygen tetrahedra forming vitreous silica. The rings are shortest closed chains of tetrahedra (or, equivalently, O\\Si units) in the bulk of the matrix without chords. The distribution of rings over the sizes (numbers of units) is given in Table 5. For convenience, Fig. 3 shows the histogram of the distribution averaged over six samples and distributions taken from the literature.

3.2. Parameterization of the potential describing the interaction of a helium atom with the matrix According to the data presented in [40] the contribution of one oxygen ion to a potential energy of helium atom placed in the ring center for each siloxane ring was calculated using the procedure described in Subsection 2.2. The rings have different diameters. This makes it possible to plot the dependence of the potential energy of the helium atom versus the distance to the oxygen ion. The obtained dependence on a logarithmic scale (energy) is shown in Fig. 4. It is seen that the points lie near a straight line; i.e., this dependence is well described by the exponential potential of the Born–Mayer type (Eq. (2)). The result of the fitting gives the following coefficients of the Born–Mayer potential: A = 2050.84 eV and B = 4.35448 Å−1.

3.3. Diffusion coefficients for helium in vitreous silica

Fig. 2. Examples of time dependences of the mean square displacement of helium atoms. Numbers at the curves are temperatures in units of Kelvin.

The diffusion coefficients for helium in the range from 500 to 1300 K at atmospheric pressure were calculated for each of the six samples of vitreous silica. The coefficients for the each temperature point were averaged over the ensemble for six samples. The results of the calculations are shown in Fig. 5 (asterisks). For comparison, this figure shows the experimentally measured diffusion coefficients (filled circles) obtained by Swets and colleagues [8]. The corresponding activation energies EA and the pre-exponential factors D0 are presented in Table 6.


51

Table 3 Positions of the peaks of the radial pair distribution functions (rdf, Å) for vitreous silica samples. Peak

Si\ \O (I) Si\ \O (II) Si\ \O (III) O\ \O (I) O\ \O (II) O\ \O (III) Si\ \Si (I) Si\ \Si (II) Si\ \Si (III)

Sample no. 1

2

3

4

5

6

1.6182 4.0469 6.3045 2.6281 5.0198 7.3483 3.1372 5.1042 7.5639

1.6189 4.0502 6.2392 2.6258 5.0159 7.3176 3.1367 5.0875 7.4292

1.6195 4.0356 6.2662 2.6254 5.0054 7.3873 3.1396 5.1012 7.6045

1.6186 4.0577 6.2281 2.6280 5.0054 7.3475 3.1399 5.1098 7.5679

1.6180 4.0524 6.2830 2.6282 5.0081 7.3657 3.1345 5.0919 7.6167

1.6182 4.0234 6.2549 2.6282 5.0057 7.3509 3.1418 5.0667 7.5816

4. Discussion

Average

Stdev

[46]

[47]

1.6186 4.0444 6.2627 2.6273 5.0101 7.3529 3.1383 5.0936 7.5606

0.0005 0.0127 0.0282 0.0013 0.0062 0.0230 0.0026 0.0155 0.0676

1.62 4.15 6.4 2.65 5.1

1.62

3.12 5.1

2.625

3.14

The fraction of 8- to 10-membered rings reaches up to 20% and more. These distributions are typical for rapidly quenched or deposited silica glass. Samples obtained and used in this work have properties as the first group of glasses and the second. On the one hand, the ring size distribution function is rather narrow (the right wing of the distribution rapidly decays, and the fraction of 8- to 10-membered rings does not exceed 10%). But the number of 6-membered rings and 7-membered ones are comparable on the other hand. The maximum of the distribution function is placed between 6 and 7. Thus, molecular dynamics simulation with empirical potentials (Eq. (1)) makes it possible to obtain within the NPT-process the samples of vitreous silica with the correct volume of interstitial voids and with the main structural parameters correlated well with the experimental data available in the literature (Tables 2–4). All samples used in this study were prepared by the same cooling rate −5.7 · 1011 K/s. Density is weakly dependent on the cooling rate when the quenching rate becomes less than 1014 K/s. For example, when the cooling rate decreases by the three orders of magnitude, glass density decreases by only 3–4% [26]. However, the quenching rate remains extremely high in MD simulation. It is not critical for the short-range order parameters, but it significantly affects the middlerange order parameters. The truncating of tail of empirical potentials can to improve the situation and to get the correct values of a number of parameters that characterize the middle order, even at such an extremely high cooling rate. These parameters include the density, position of 2 and 3 peaks of RDF and the Si\\O\\Si bond-angle distribution. Unfortunately, the literature contains no direct experimental data for other middle-range parameters, characterizing the geometry of the interstitial space in greater detail (for example, the rings size distribution). This makes it impossible to clearly prove or disprove the validity of the models prepared in this work for gas diffusion problems. Nevertheless, a good agreement of the calculated diffusion coefficients with experimental ones may be considered as some indirect argument for the adequacy of the model used. The correct parameterization of empirical potential approximated the interaction between gas atoms migrated and matrix is the second significant factor determined the quality of the simulation of the gas diffusion through non-porous media. The configuration containing a helium atom in the center of a flat siloxane ring is a good model of the transition state of a helium atom jumping through the doorway connecting the adjacent solubility sites. Therefore, the parameterization based on quantum mechanical calculations of these configurations gives

Molecular dynamics simulation allows to create molecular models for the new functional materials and to test their functionality without using of any a priori empirical data relating to the structure of materials. This process becomes more complicated for the non-porous membrane materials designed for gas separation, because of the high sensitivity of the diffusion process to modeling quality of matrix and gas-matrix interaction. Therefore in this paper special attention was paid to the characteristics of the matrix models and the parameterization of the interaction potentials between the gas atoms and the matrix. The preparation of vitreous silica samples using molecular dynamics simulation models in essence real melting glass manufacturing process at atmospheric pressure, i.e. within the NPT-process. The value of the density and volume of interstitial voids of glass obtained is the result of vitrification rather than an external conditions (for example, as in NTV process). Morse potential with truncated tail in the form (1) makes it possible to obtain models of silica glasses, the density of which is in good agreement with the experimental value. This fact is illustrated in Table 2. In addition, as one can see from Tables 3 and 4, the short-range order parameters (positions of the first peaks of the pair radial distribution functions and the maxima of the angular distributions) are in very good agreement with the experimental values obtained from the diffraction data [46] and, especially, with the refined data reported in [47]. The positions of the second and third peaks of the pair radial distribution functions depend on the structural elements of vitreous silica (silicon–oxygen tetrahedra) packing manner. As can be seen from Table 3, these parameters are close to the experimental values (for all the second peaks and the third peak Si\\O (III)). Because of absence of the direct experimental data on the ring statistics in vitreous silica, the ring size distributions calculated in this paper are compared with the similar data obtained from the numerical simulation available in literature. The data can be divided into two groups. An example of the first group is presented in [48,49]. For this group of models, the ring size distribution function has a maximum corresponding to six-membered rings, and the distribution function itself is relatively narrow (i.e. the fraction of 8- to 10-membered rings is small). The structure of this type models resembles the structure of a highly distorted cristobalite. Models of the second group, on the contrary, have relatively broad size distribution functions of the rings.

Table 4 Maxima of the angular distributions (deg) for vitreous silica samples. Angle

Si\ \O\ \Si O\ \Si\ \O

Sample no. 1

2

3

4

5

6

153.23 109.06

153.21 109.01

153.05 108.92

153.83 108.95

153.51 108.93

153.86 108.98

Average

Stdev

[46]

[47]

153.45 108.97

0.3407 0.0524

144 109.5

153 109.4

52


Table 5 Size distribution of rings (%) in vitreous silica samples. Number of units

2 3 4 5 6 7 8 9 10

Sample no. 1

2

3

4

5

6

0.28 1.18 7.56 20.50 30.14 29.64 10.03 0.67 0.00

0.06 1.23 7.70 22.15 28.85 27.01 11.89 1.06 0.06

0.00 1.12 7.98 19.75 32.37 28.40 9.65 0.73 0.00

0.06 0.72 8.45 19.63 32.31 28.20 9.84 0.72 0.06

0.00 1.14 7.19 21.39 32.34 28.75 8.39 0.74 0.06

0.00 0.84 5.73 20.32 33.69 29.90 8.69 0.84 0.00

Average

Stdev

[48]

[49]

[50]

[51]

0.07 1.04 7.44 20.63 31.62 28.65 9.75 0.79 0.03

0.11 0.21 0.94 0.98 1.77 1.05 1.24 0.14 0.03

1.05 9.06 24.39 36.76 20.56 6.97 1.05 0.17

1.43 8.67 25.47 36.37 17.69 6.26 2.86 1.25

0.39 12.92 16.53 23.54 22.17 16.50 4.87 3.08

0.54 5.81 18.12 26.13 26.40 14.89 5.65 2.31

better quality empirical potentials for the simulation of the gas diffusion in a non-porous matrix. A combination of the correct vitreous silica model in terms of the interstitial void space with the good-quality parameterization of the helium-matrix interaction potentials made it possible to directly simulate the diffusion process by the molecular dynamics method and to obtain numerical values of the diffusion coefficients, which, as can be seen from Fig. 4 and Table 6, correlate with the experimental data over a wide temperature range as well. This result shows that the molecular dynamics simulation can be effectively used for end-to-end computer-aided design of the new membrane materials: from the structure's development and synthesis of model samples, up to their diffusion properties testing. It cannot be too highly stressed that this method no requires any a priory empirical data on the structure of the material tested. 5. Conclusion The ability of the molecular dynamics simulation for quantitative simulation of gas diffusion in non-porous materials was tested by the using the well-studied process of helium diffusion in vitreous silica. The correct

Fig. 3. Ring size distribution functions.

Fig. 5. Calculated (asterisks) and experimental (points) temperature dependences of the diffusion coefficient of helium in vitreous silica.

Table 6 Activation energy EA of helium diffusion and pre-exponential factor D0.

Fig. 4. Dependence of the potential energy of the interaction of a helium atom with an oxygen ion according to the calculation from the data reported in [51] (points) and the results of fitting the Born–Mayer potential (line).

Activation energy EA Pre-exponential factor D0

Calculation (this work)

Experiment [8]

25.019 kJ/mol 4.145 · 10−4 cm2/s

25.368 kJ/mol 5.458 · 10−4 cm2/s


interstitial voids volume value of the glass were prepared by using molecular dynamics simulation of quenching of molten quartz at a constant pressure, using Truncated Morse potentials. Parameterization of the gas-matrix atoms interaction potential made for configurations similar to ones corresponding to the migrating gas atom jump from one solubility site to an adjacent one. It has been shown that taking into account these two aspects allows for quantitative modeling of diffusion of helium in vitreous silica, calculate the diffusion coefficients being close to the experimental data. The method not requires any a priory data on the structure of designed material to test its diffusion properties. Thus, the molecular dynamics simulation can be effectively used for end-to-end computer-aided design of new membrane materials.

[18] [19] [20] [21]

Acknowledgments

[32]

This study was supported by the Russian Science Foundation (project no. 14-13-00289). References [1] S.A. Stern, T.F. Sinclair, P.J. Gareis, N.P. Vahldieck, P.H. Mohr, Ind. Eng. Chem. 57 (1965) 49. [2] E.V. Fomenko, E.S. Rogovenko, L.A. Solovyov, A.G. Anshits, RSC Adv. 4 (2014) 9997. [3] W.H. Zachariasen, J. Am. Chem. Soc. 54 (1932) 3841. [4] A.C. Wright, M.F. Thorpe, Phys. Status Solidi B 250 (2013) 931. [5] A.C. Wright, J. Non-Cryst. Solids 401 (2014) 4. [6] P.L. Studt, J.F. Shackelford, R.M. Fulrath, J. Appl. Phys. 41 (1970) 2777. [7] J.S. Masaryk, R.M. Fulrath, J. Chem. Phys. 59 (1973) 1198. [8] D.E. Swets, R.W. Lee, R.C. Frank, J. Chem. Phys. 34 (1961) 17. [9] G.G. Boiko, G.V. Berezhnoi, Glas. Phys. Chem. 29 (2003) 42. [10] L.V. Woodcock, C.A. Angell, P. Cheeseman, J. Chem. Phys. 65 (1976) 1565. [11] T.F. Soules, J. Chem. Phys. 71 (1979) 4570. [12] T.F. Soules, J. Non-Cryst. Solids 49 (1982) 29. [13] J.D. Kubicki, A.C. Lasaga, Am. Mineral. 73 (1988) 941. [14] C.A. Angell, J.H.R. Clarke, L.V. Woodcock, Adv. Chem. Phys. 48 (1981) 397. [15] S. Tsuneyuki, M. Tsukada, H. Aoki, Y. Matsui, Phys. Rev. Lett. 61 (1988) 869. [16] W. Kob, J. Phys. Condens. Matter 11 (1999) 85. [17] B. Vessal, M. Leslie, C.R.A. Catlow, Mol. Simul. 3 (1989) 123.

[22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

[33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51]

53

B.W.H. van Beest, G.J. Kramer, R.A. van Santen, Phys. Rev. Lett. 64 (1990) 1955. S. Tsuneyuki, Mol. Eng. 6 (1996) 157. A. Takada, P. Richet, C.R.A. Catlow, G.D. Price, J. Non-Cryst. Solids 354 (2008) 181. T.F. Soules, G.H. Gilmer, M.J. Matthews, J.S. Stolken, M.D. Feit, J. Non-Cryst. Solids 357 (2011) 1564. K. Vollmayr, W. Kob, K. Binder, Phys. Rev. B 54 (1996) 15808. K. Yamahara, K. Okazaki, K. Kawamura, J. Non-Cryst. Solids 291 (2001) 32. A. Takada, P. Richet, C.R.A. Catlow, G.D. Price, J. Non-Cryst. Solids 345 (2004) 224. N. Kuzuu, K. Nagai, M. Tanaka, Y. Tamai, Jpn. J. Appl. Phys. 44 (2005) 7550. S.V. Kukhtetskiy, Modern Problems of Science and Education 1, 2015. S.V. Kukhtetskiy, Modern Problems of Science and Education 2, 2015. A. Pedone, G. Malavasi, M.C. Menziani, A.N. Cormack, U. Segre, J. Phys. Chem. B 110 (2006) 11780. HOOMD-blue web page: https://codeblue.umich.edu/hoomd-blue/ J.A. Anderson, C.D. Lorenz, A. Travesset, J. Comput. Phys. 227 (2008) 5342. J. Glaser, T.D. Nguyen, J.A. Anderson, P. Lui, F. Spiga, J.A. Millan, D.C. Morse, S.C. Glotzer, Comput. Phys. Commun. 192 (2015) 97. D.N. LeBard, B.G. Levine, P. Mertmann, S.A. Barr, A. Jusufi, S. Sanders, M.L. Klein, A.Z. Panagiotopoulos, Soft Matter 8 (2012) 2385. CUDA web page: http://www.nvidia.com/object/cuda_home_new.html N. Wilt, The CUDA handbook, A Comprehensive Guide to GPU Programming, 2013 ISBN-13: 978-0321809469. C. Chakravarty, K.V. Thiruvengadaravi, Proceedings of the Conference on Frontiers in Materials Modelling and Design, Kalpakkam, August 1996 20–23. O. Talu, A.L. Myers, Colloids Surf. A Physicochem. Eng. Asp. 187 (2001) 83. M. Reich, R.C. Ewing, T.A. Ehlers, U. Becker, Geochim. Cosmochim. Acta 71 (2007) 3119. B. Guillot, N. Sator, Geochim. Cosmochim. Acta 80 (2012) 51. B. Vessal, M. Amini, C.R.A. Catlow, Mol. Simul. 15 (1995) 123. P. Hacarlioglu, D. Lee, G.V. Gibbs, S.T. Oyama, J. Membr. Sci. 313 (2008) 277. J.F. Shackelford, Int. J. Appl. Glas. Sci. 2 (2011) 85. R.H. Doremus, J. Amer. Ceram. Soc. 49 (1966) 461. VMD web page: http://www.ks.uiuc.edu/Research/vmd/ W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graphics Modell. 14 (1996) 33. N.P. Bansal, R.H. Doremus, Handbook of Glass Properties. NY, 1986. R.L. Mozzi, B.E. Warren, J. Appl. Crystallogr. 2 (1969) 164. J.R.G. Da Silva, D.G. Pinatti, C.E. Anderson, M.L. Rudee, Philos. Mag. 31 (1975) 713. J.P. Rino, I.I. Ebbsjo, R.K. Kalia, A. Nakano, P. Vashishta, Phys. Rev. B 47 (1993) 3053. F. Barmes, L. Soulard, M. Mareschal, Phys. Rev. B 73 (2006) 224108. S. Kohara, K. Suzuya, J. Phys. Condens. Matter 17 (2005) S77. B. Guillot, Y. Guissani, J. Chem. Phys. 105 (1996) 255.