Empirical Potential Molecular Dynamics

[prev] [prev-tail] [tail] [up]

3.4 Empirical Potential Molecular Dynamics

Molecular dynamics has been established as a powerful tool for the generation of amorphous structures [151, 152]. It can simulate the time evolution of a group of atoms at a certain temperature, where the bonding between atoms is mimicked by interatomic empirical potentials. Even though it performs considerably fast, the simulation times are still restricted to a few thousand picoseconds [153, 154, 155, 156]. For this reason, empirical potential molecular dynamics is not capable of simulating the processing of $a$ - $SiO2$ . Nevertheless, a combination of experimental and theoretical investigations have shown that realistic amorphous structures [151, 152] can be produced by cooling down a random configuration of silicon and oxygen atoms from $3000K$ to room temperature within a few tens of a picosecond. It has been found [151, 152, 153, 155] that $a$ - $SiO2$ is composed of slightly deformed tetrahedral $SiO4$ units with one $Si$ atom in their centers. These units are randomly connected to each other so that they form $Si$ - $O$ - $Si$ chains at their corners. In this way, each silicon atom is fourfold coordinated to oxygen atoms and each oxygen atom in turn is bonded to two silicon atoms. The distributions of the $dSi-O$ , $dO-O$ , and $dSi-Si$ bond lengths as well as the $ϕO -Si-O$ and $ϕSi-O-Si$ angles have been used in the following to check the quality of the generated $a$ - $SiO2$ structures.

3.4.1 Fundamentals of Molecular Dynamics

The atomistic dynamics are accurately described by Newton’s law of motion, which is applied for classical molecular dynamics [157].

mi ddtRi(t) = Fi(t) (3.39) F (t) = - ∇ V (R (t),...,R (t)) (3.40) i i 1 N

$Ri(t)$ and $mi$ denotes the position and the mass of atom $i$ . The term on the right-hand side of equation (3.39) represents the force $Fi(t)$ acting on the atom $i$ and is evaluated by the derivative of the interatomic empirical potential $V (R1(t),...,RN(t))$ with respect to $Ri(t)$ . This differential equation is solved numerically using an appropriate time integration algorithm, such as the leap-frog Verlet algorithm.

Ri(t+ Δt) = Ri(t)+ vi(t+ 1 Δt)⋅Δt (3.41) 2 vi(t+ 1Δt ) = vi(t- 1Δt)+ Fi-⋅Δt (3.42) 2 2 mi

In this procedure, the current positions $Ri(t)$ and the accelerations $Fi∕mi$ are stored together with the mid-step velocities $vi(t- 12Δt)$ . The structure generation method used in this thesis is based on a rapid quench of a molten atomic system. Therefore, a thermostat is required to control the temperature of the atomic system during the quenching procedure. For this purpose, the Nose-Hoover thermostat [157] as implemented in the GULP code was employed throughout this thesis. It relies on a sophisticated method to couple the atomic system to a heat bath with the desired temperature. Since this method correctly produces the thermodynamical temperature fluctuations as well as the dynamics of the atomic system, the Nose-Hoover thermostat is usually considered as the working horse for molecular dynamics simulations.

3.4.2 Procedure for Structure Generation

The silicon and oxygen atoms were randomly placed in the periodic simulation cells. In order to avoid any overlapping between the atoms, exclusion radii ( $rSi-O = 1.5˚A$ , $rO-O = 2.5˚A$ , $rSi-Si = 3.1˚A$ ) were used. The edge length of the simulation cells ( $11.79˚A$ ) was chosen to match a mass density of $2.19g∕cm3$ [155]. The resulting random structures were taken as a starting configuration for the subsequent molecular dynamics equilibration step, which was performed at $3000K$ for $30ps$ with a time step of $1fs$ . In this step the atomic structure is evolved from an unnatural random configuration to a liquid that should resemble molten $SiO2$ . It was followed by a quenching step to $0K$ for $30ps$ with a time step of $1fs$ , where the liquid was cooled down to an amorphous solid.

The simulations were performed using the popular Beest-Kramer-van Santen (BKS) potential [158]. This consists of Buckingham potentials, which were extended by a Coulombic term and parametrized to reproduce the interatomic interactions obtained from DFT. These two-body potentials feature artificial singularities at their origins and small separating barriers to the next energy minimum. But since the structures were not heated above $5000K$ , corrections within a certain cut-off radii as applied in [159] could be omitted. The interatomic interactions were only represented by $Si$ - $O$ and $O$ - $O$ pair-potentials that describe the $Si$ - $O$ bonding and ensure the tetrahedral arrangement. Despite these strong simplifications, a series of studies have proven their successful application for $SiO2$ structure generation [158, 159, 152, 151].

In order to prove the correctness of the applied production procedure, the obtained samples were evaluated based on the pair-correlation functions and angle distributions as shown in Fig. 3.5 and 3.6. Due to the fact that edge-sharing tetrahedra are energetically unfavored [152], only samples containing none of these edge-sharing tetrahedra were used for further investigations while the others were simply discarded. The remaining samples exhibited no miscoordination, such as broken $Si$ - $O$ bonds or threefold coordinated $Si$ atoms. As demonstrated in Table 3.1, satisfying agreement has been achieved with previously published results [152, 153, 155]. The selected structures were minimized on a DFT level in order to prepare them for the following defect calculations. During this step, a small structural relaxation was observed indicating that no bonds had been broken.


Ref.	$dSi- O$	$dO -O$	$dSi-Si$	$ϕO-Si-O$	$ϕSi-O -Si$


Present study	1.64	2.66	3.08	109.42	142.62
[152]	1.62	2.64	3.10	109.6	142.0
[155]	1.63	2.67	3.11	109.4	146.8
[153]	1.62	2.68	2.98	109	136

Table 3.1: Comparison of the characteristic properties of the produced $a$ - $SiO2$ structures. $dSi-O$ , $dO-O$ , and $dSi- Si$ denote the first maximum in the corresponding pair-correlation functions. $ϕO-Si-O$ and $ϕSi- O-Si$ are the maxima of the respective angle distributions. The obtained values compare reasonably well with the published values in [152, 155, 153]. It should be mentioned here that the $ϕO -Si-O$ angles are quite sensitive to the details of the used structure generation method. Thus, their values in the literature [152, 153, 155, 159, 156] are subject to a strong variation ranging from $136∘$ and $152∘$ and thus are still debated [156].

Figure 3.5: The pair-correlation functions [151] for the $dSi-O$ (left), the $dO- O$ (middle), and the $dSi-Si$ (right) bonds. The solid black lines represent the data obtained in this thesis while the dashed ochre and the dashed-dotted green line are extracted from the studies of Sarnthein [153] and Rino [152], respectively. In order to improve the accuracy of statistics, the data have been collected from a $1ps$ DFT molecular dynamics run at $300K$ for five different samples. The first sharp peak corresponds to the length of the respective bond type and compare well with the values from literature (see Table 3.1). It is noted that also the other features are found to be in qualitative agreement with the extracted data. The integral over the first peak yields the number of first nearest neighbors for the respective atom type. The values $4.0$ , $6.19$ , and $4.05$ have been obtained for the $dSi- O$ , the $dO-O$ , and the $dSi- Si$ bonds in agreement with [151, 155, 153]. The slight deviations have been attributed to the still small statistics associated with a small resolution with respect to $r$ . Furthermore, the first peaks overlap with the next features at the right integration limit (indicated by the arrow) leading to an inaccurate determination of the numbers of first nearest neighbors.

Figure 3.6: The angle distribution of the $O$ - $Si$ - $O$ (left) and the $Si$ - $O$ - $Si$ (right) chains for the data of this thesis and (ochre) Sarnthein [153] and (green) Rino [152]. The maxima as well as the full width at half maximum of the $O$ - $Si$ - $O$ angle distribution are in reasonable agreement with the extracted data in [152, 153]. The values for the $O$ - $Si$ - $O$ angles are subject to appreciable deviations, which originate from the complications mentioned in Table 3.1. Nevertheless, the maxima of Gaussian fits are found to be in an acceptable agreement.

[prev] [prev-tail] [front] [up]