Simulation of borosilicate glasses with non-constant force field molecular dynamics

In this study the simulation of microscopical behavior of borosilicate glasses was conducted with non-constant force field molecular dynamics. The suggested model consists of classical pair potentials in the Buckingham form, long range Coulomb interaction, intramolecular bonded interactions and possibility of bond breaking and formation. The latter effects are accompanied by changes in the types of the bond-forming particles. The simulated system corresponds to the structure of borosilicate glasses with predominantly four-coordinated boron atoms. Different structure groups are formed due to the dissociation / formation of intramolecular bonds, and the processes of the glass network rearrangement intensifies with temperature increasing.


Introduction
Borosilicate glasses are materials with a wide range of applications.They can be used for nuclear waste utilization [1], as implants for tissue repair [2], in electrochemical devices as ionic conductors [3][4][5], as sealants for solid oxide fuel cells [6] and in many other areas of human activity.From this point of view, microscopical behavior of these glasses comes to researchers' attention.Glasses consisting of several network former oxides (B 2 O 3 , SiO 2 , P 2 O 5 and some others) can contain a colossal amount of possible structure groups and their combinations.For this reason numerical simulations of mixed glasses are very promoting because it can "view" glass structure at molecular level.To date, a molecular dynamics (MD) simula-tion of silicate, borate and mixed borosilicate glasses was conducted in series of papers [1][2][3][4][5][7][8][9][10][11][12][13][14][15][16][17][18][19].Most of these works deal with classical MD, except researches [2,17] performed with Car-Parinello (MD) [20] and reax force field [21] (ReaxFF), correspondingly.The Car-Parinello method brings quantum mechanical calculation into MD simulation.The ReaxFF is constucted to provide chemical reactivity during MD simulation.The last two methods are used to include effects of chemical bonds dissociation and recombination to the simulation.These effects are needed to be considered for correct simulation of glass behavior.To perform such calculations with classical MD is an ambiguous and complicated, if not impossible, task.
In the present paper I demonstrate how one can perform MD simulation of borosilicate glass with non-constant force field methodology.This approximation consists of possibility of bond breaking and formation (including changing in particle types of bond forming atoms) runtime.

Simulation details
MD simulations were performed with an original program written in C in canonical (NVT) ensemble.Newton's equations of motion were integrated by Velocity Verlet algorithm [22] with timestep of 0.5 fs during 100'000 steps.Equilibration time was 0.5 ps (1000 timesteps).Electrostatic interactions were accounted using the Ewald summation.Nosé-Hoover thermostat [23] with relaxation time of 0.02 ps was used for maintaining the temperature around desired temperature.There was tested three temperatures: 298, 500 and 1000 K.The considered system consists of 94 Si atoms, 80 B atoms, 429 O atoms and 242 alkaline metal (Me) atoms (the total number of particles is 845).The oxygen atoms were divided into bridging (O c ) and non-bridging (O t ) species.Bridging oxygen connects with two atoms of glass-forming element (Si / B) and non-bridging -with only one.Initial configuration of the system was generated with own special code.The box was cubic with the edge length of 20.4 Å to match an approximate density of glassy silicates.Van der Waals interactions in the system were given by Buckingham pair potential: where U ij is the potential energy between the i-th and the j-th atoms, r ij is the distance between them, and A ij , ρ ij and C ij are empirical parameters.Besides the Van der Waals short-range interactions, covalent bonds between oxygen and silicon / boron were set in the form of an intramolecular harmonic potential: where k ij is the spring constant and r 0ij is the equilibrium distance between particles.Unlike interactions provided by Eq. ( where k tb is the spring constant of the threebody potential, 100.0 eV, θ ijk is the angle between ij and ik bonds, θ 0 is the equilibrium angle, 148.3° (taken from [15] for Si -O -Si angle); i, j and k are the indexes of O c atom and its covalent-bonded neighbors, correspondingly.
The values of A ij , ρ ij and C ij were taken from papers [1,11,12,15,19].Simulations were performed with different combination of these parameters.The more suitable for the glass structure description set of the parameters is summarized in Table 1 with corresponding references.The parameters of the valent bond potential (2) are given in Table 2.The atomic charges were suggested to be found as 4δ, 3δ, a, -2δ and -δ-a for Si, B, Me, O c and O t , correspondingly.The values of δ and a are set to 0.4e and 1.0e.Our MD program allows deleting (adding) valent bonds from (to) the corresponding list to mimic processes of bonds dissociation and glass forming: and If the bonds automatically break (form) then the interatomic distance is higher (lower) than a maximal bond distance parameter (r m ).The values of this parameter are given in Table 2.The further dissocia-tion with participation of Si+ or B+ species is not allowed.The charges of Si+ and B+ species are chosen to keep electroneutrality during reactions (4) and ( 5).

Results and Discussion
The obtained structure of the system is presented in Fig. 1.One can see that the simulated glass consists of branchy chains of Si -O / B -O bonds which can form various combinations and loops.In principle, there are two possible boron coordinations in borate glasses, triangular and tetrahedral [24].In our case almost all boron in the glass is four-coordinated, Fig. 1.
Radial distribution functions (RDFs) are presented in Fig. 2. All RDF-curves consist of one sharp and several broad maximums.This is typical of liquids or amorphous systems and indicates the presence of the short-range and the absence of the long-range orders.There are no series of well-resolved maxima, specific for crystalline solids.The first maximum of RDF-curve corresponds to the most probable distance in the near-neighbor coordination.According to Fig 2(a), the distance between Me and O t is less than that between Me and O c .It can be explained by more negative charge of O t species (Simulation Details section) and, therefore, stronger Coulomb interaction.RDFs almost do not depend on temperature, the differences are observed only for the RDFs with the sharpest maxima, for example, Si -O.For these RDFs the value of the first maximum decreases with temperature increasing and the peak becomes broader, Fig. 2 As the system initially does not have crystalline lattice and a corresponding look of RDF-curves, it is necessary to find another criterion of transition between liquid and solid states.Such a criterion can be, for example, mean square displacement (MSD), which characterizes a mobility of ions.In the simulated borosilicate glass, MSDs of most ions tend to reach a plateau with some oscillation at temperature of 298 K, except one for Me ions increased linearly, Fig. 3(a).The time derivative of MSD grows with temperature increasing, Fig. 3(b).This means that all species start to diffuse and the system behavior becomes liquid-like.Analogous effects can be seen on trajectories of the ion motion.At temperature of 298 K most ions vibrate around one point (although with big enough amplitude), Fig. 4(a).An exception is the Me ion, its trajectory consists of a series of hops between positions with long enough oscillation time.With temperature increasing, Me ion spends less time in one position, more often hops and covers greater distance, Fig. 4(b).Its movement character becomes more and more liquid-like.Fig. 5 demonstrates the number of Si and B species as a function of time.The changes in these quantities are provided by dissociation of B -O c and Si -O c bonds according to equations ( 4) and (5).Since dissociation in the suggested model occurs upon reaching the determined distance (Table 2), the observed fluctuations of the particles numbers are related with changing in the B / Si -O c distances.According to Fig. 5, the suggested pair potentials provide stronger B -O c bonding, then Si -O c one, as the amplitude of the Si number oscillation is bigger.Also, Fig. 5 demonstrates that the amplitudes of these oscillations increase with temperature.This is a result of thermal motion which promotes overcoming of forces of interatomic attraction, as noted above in the discussion of the radial distribution functions.The bonds breaking and formation realized by methodology of non-constant force field allow glass structure to be "dynamical, " i.e. to destroy ones atomic groups and to cre-ate others.A variety of formed molecular groups can see in Fig. 1.

Conclusions
In the paper reference classic pair potentials have been tested for molecular dynamic simulation of borosilicate glasses.The suggested combination of the potentials and ionic charges describes interatomic distances in borosilicate glasses reasonably well.Oxygen atoms were divided into bridging and non-bridging by the value of electrical charge.The transition between these states is possible due to bond formation / dissociation.These phenomena are implemented by non-constant force field molecular dynamics.It has been shown that the degree of the bond dissociation increases with temperature growth.
), this potential is applied only to covalent bonded pairs (S -O t , Si -O c , B -O t and B -O c ) listed separately.In addition, a three-body potential energy term was applied to maintain a feasible valent angle for Si -O c -Si, B -O c -B and Si -O c -B bonds: (b).This process occurs, most likely, due to the thermal motion.The first maximum for B -O t / O c and Si -O t / O c pairs locates at distance of 1.375-1.425and 1.575-1.675Å, correspondingly, which is close to B -O и Si -O interatomic distances in borosilicate glasses [15].

Fig. 3 .Fig. 4 .
Fig. 3.The mean square displacement (MSD) of the simulated borosilicate glass: (a) for different ions at 298 K; (b) for Me and Si ions at different temperatures

Fig. 5 .
Fig. 5.The number of Si and B species as a function of simulation time