Factors Affecting Molecular Dynamics Simulated Vitreous Silica Structures

Ersan Demiralp, Norman T. Huff*, Tahir Cagin, William A. Goddard III

Materials & Process Simulation Center, 400 South Wilson, Beckman Institute (139-74),

California Institute of Technology, Pasadena, CA 91125

* Owens Corning Science & Technology Center, 2790 Columbus Road, Granville, OH 43023-1200

Abstract

Molecular dynamics (MD) simulations are being used to gain insights into the atomic level structure of glasses. However, the particular MD cooling cycle algorithms used to generate the simulated structures can influence some calculated structural parameters. In this work, the magnitude of the differences in structural characteristics are studied as a function of the cooling cycle algorithms used to generate simulated glass structures. It is shown that run to run variations for MD simulations containing several hundred atoms can be significant. As a result, one should be extremely cautious in drawing conclusions regarding structural characteristics based upon single simulations.

Introduction

Over the past 15 years, the use of MD techniques to simulate glass structures has become an accepted technique for gaining insight into the structure of glasses at an atomistic level. However, the manner in which the MD simulation is carried out can have an impact upon the observed glass structure. Typically, these simulations are carried out by performing a series of MD simulations on a 3-dimensional periodic system using a constant number of particles in the model cell, a constant volume for the model cell, and a constant average temperature of the particles in the model cell (NVT dynamics). The simulation begins by performing an NVT MD simulation for tens of picoseconds at some high temperature (typically 6000 K or higher) in an attempt to remove any "memory" of the starting arrangement of the atoms. The temperature of the system is then reduced in steps of several hundred to perhaps 2000 K until ambient temperatures are reached. The NVT simulation time at each of the temperature steps is typically tens of picoseconds. As one would expect, the structures generated with these MD simulations are sensitive to the particular temperatures and times used in the multi-step simulation. However, there has been no systematic published study of the effects of the cooling cycle, initial structure, and cell size upon the final structural parameters. As a result, different researchers utilize different cooling algorithms. The implicit assumption made is that significant results are independent of the cooling cycle. The aim of this paper is to quantify the impact of the particular cooling cycle and cell size upon the simulated room temperature glass structures.

Procedure

In this paper, all simulations are performed upon silica. The program utilized is Release 3.0 of Cerius^{2} by Molecular Simulations Incorporated. The standard starting cell contains 648 atoms (216 Si and 432 O atoms). This cell is a cube, 21.48 Å on a side. The standard starting cell that contained 1536 atoms is a cube, 28.64 Å on each side. The initial arrangement of atoms for both sizes of cells is that of crystalline cristobalite.

NVT dynamics are utilized for most of the annealing cycle. The volume utilized for the NVT dynamics corresponds to a silica density of 2.17454 g/cc. For many force fields, this results in very high stress in the unit cell at 300 K, typically more than 2 Gpa. As a result, in this work the NVT dynamics is followed by NPT dynamics during which time the structure relaxes to the density which the force field predicts for atmospheric pressure.

Most of the calculations are carried out using the glass force field 2.01 supplied with the Cerius2 software. This force field utilizes a 3 body term similar to that employed by Garofalini et. al. In addition, some calculations are performed using force fields developed by Tsuneyuki and by Demiralp.

Two types of cooling cycles are utilized. These are presented in Tables I and II. In Table I, Cycle 1-I consisted of 20 ps of NVT dynamics at 8000 K followed by 20 ps of NVT and 4000 K, 2000 K and 1000 K. This is followed by 10 ps of NVT dynamics at 300 K and then 30 ps of NPT dynamics at 300 K. In Table II, the between 4000 K and 300 K, dynamics are performed at 100 K intervals from 3900 K to 300 K for the indicated time. For example, in Cycle 2-I, 4 ps of NVT dynamics are run at 3900 K, 4 ps at 3800 K, 4 ps at 3700 K, etc. to 400 K.

TABLE I: Type 1 Cooling Cycles

Temperature |
Cycle 1-I |
Cycle 1-II |
Cycle 1-III |
Cycle 1-IV |
Cycle 1-V |
Cycle 1-VI |

8000 K NVT |
20 ps |
40 ps |
40 ps |
40 ps |
0 ps |
0 ps |

6000 K NVT |
0 ps |
0 ps |
0 ps |
0 ps |
0 ps |
40 ps |

4000 K NVT |
20 ps |
0 ps |
40 ps |
40 ps |
40 ps |
0 ps |

2000 K NVT |
20 ps |
0 ps |
20 ps |
40 ps |
0 ps |
0 ps |

1000 K NVT |
20 ps |
0 ps |
0 ps |
0 ps |
20 ps |
20 ps |

300 K NVT |
10 ps |
10 ps |
10 ps |
10 ps |
10 ps |
10 ps |

300 K NPT |
30 ps |
30 ps |
30 ps |
30 ps |
30 ps |
30 ps |

TABLE II: Type 2 Cooling Cycles

Temperature |
Cycle 2-I |
Cycle 2-II |
Cycle 2-III |
Cycle 2-IV |
Cycle 2-V |

8000 K NVT |
20 ps |
20 ps |
20 ps |
40 ps |
40 ps |

4000 K NVT |
60 ps |
60 ps |
20 ps |
20 ps |
0 ps |

100 K steps |
4 ps |
1 ps |
1 ps |
1 ps |
1 ps |

300 K NVT |
10 ps |
10 ps |
10 ps |
10 ps |
10 ps |

300 K NPT |
30 ps |
30 ps |
30 ps |
30 ps |
30 ps |

Structural Characterization

The structures are characterized by calculating the Si-O-Si and O-Si-O bond angle distributions, the coordination number of Si and O atoms and the density of the glass after NPT dynamics. The bond angles and coordination numbers are calculated based upon an algorithm by Garofalini.

Results

The results are shown in Tables III and IV. The first column shows which of the cooling cycles shown in Tables I and II are used. Each cooling cycle appears twice showing the results for each of the duplicate runs. The columns labeled O-Si-O and Si-O-Si indicate the location of the peak and the width of the distribution in the O-Si-O and Si-O-Si bond angles. The location of the peak and width of the peak at one-half the maximum peak height are calculated using a one degree increment for the angle "bins". The next two columns show the percentage of Si and O atoms which are 4 and 2 coordinated respectively. Two atoms are considered coordinated if they are within 1.1 times the sum of the van der Waals radii of each other. The next column, labeled Density Range, shows the range in the densities averaged over 200 fs increments over the last 15 ps of the NPT simulations. The column labeled average density is the average density of the structure over the last 15 ps of the NPT simulations.

Table III: Effect of Cooling Cycle on Structural Characteristics for MSI Glass FF 2.01

Cooling Cycle |
O-Si-O peak°/fwhm° |
Si-O-Si peak°/fwhm° |
% Four Coord Si |
% Two Coord O |
Density Range |
Average Density |

1-I |
108/15 |
154/35 |
97.5 |
96.0 |
2.142-2.163 |
2.15005 |

1-I |
107/15 |
161/35 |
98.4 |
97.0 |
2.141-2.167 |
2.15293 |

1-II |
108/15 |
150/35 |
96.6 |
93.8 |
2.210-2.232 |
2.22012 |

1-II |
108/17 |
156/38 |
96.8 |
94.6 |
2.210-2.238 |
2.22253 |

1-III |
108/15 |
155/33 |
99.2 |
97.6 |
2.161-2.187 |
2.17403 |

1-III |
107/14 |
156/33 |
98.9 |
97.8 |
2.104-2.128 |
2.11527 |

1-IV |
108/17 |
154/37 |
98.8 |
95.8 |
2.167-2.192 |
2.17688 |

1-IV |
108/15 |
153/36 |
98.5 |
94.5 |
2.131-2.163 |
2.14837 |

1-V |
108/5 |
175/7 |
100 |
100 |
1.898-1.940 |
1.92004 |

1-V |
108/5 |
175/7 |
100 |
100 |
1.897-1.943 |
1.92017 |

2-I |
107/14 |
150/36 |
99.3 |
98.5 |
2.129-2.154 |
2.14216 |

2-I |
107/15 |
153/34 |
99.3 |
98.6 |
2.210-2.238 |
2.22253 |

2-II |
107/15 |
150/34 |
99.6 |
98.3 |
2.153-2.176 |
2.16453 |

2-II |
107/14 |
152/33 |
99.8 |
98.2 |
2.126-2.148 |
2.13646 |

2-III |
107/14 |
155/32 |
98.8 |
97.9 |
2.122-2.150 |
2.13599 |

2-III |
107/16 |
155/29 |
99.0 |
98.6 |
2.155-2.177 |
2.16600 |

2-IV |
108/18 |
157/35 |
98.8 |
96.8 |
2.188-2.212 |
2.20108 |

2-IV |
109/14 |
158/33 |
99.1 |
98.1 |
2.149-2.171 |
2.16089 |

2-V |
107/16 |
155/34 |
98.4 |
96.1 |
2.167-2.189 |
2.17725 |

2-V |
108/15 |
157/34 |
98.0 |
97.2 |
2.166-2.190 |
2.17742 |

Table IV: Effect of Force Field and Cell Size Upon Structural Characteristics Using Cooling Cycle 2-I

Cooling Cycle |
O-Si-O peak/fwhm |
Si-O-Si peak/fwhm |
Four Coord Si |
Two Coord O |
Density Range |
Average Density |

MSI gff 2.01 (648) |
107/14 |
150/36 |
99.3 |
98.5 |
2.129-2.154 |
2.14216 |

MSI gff 2.01 (648) |
107/15 |
153/34 |
99.3 |
98.6 |
2.210-2.238 |
2.22253 |

MSI gff 2.01 (1536) |
107/14 |
153/36 |
99.3 |
98.6 |
2.157-2.169 |
2.16278 |

MSI gff 2.01 (1536) |
107/14 |
155/30 |
99.1 |
97.9 |
2.148-2.160 |
2.15476 |

Tsuneyuki (648) |
106/16 |
143/38 |
99.9 |
99.6 |
2.276-2.324 |
2.29877 |

Tsuneyuki (648) |
107/15 |
146/40 |
99.6 |
99.7 |
2.291-2.336 |
2.31401 |

Demiralp (648) |
107/18 |
144/35 |
99.2 |
96.7 |
2.236-2.284 |
2.26411 |

Demiralp (648) |

**Discussion**

The most prominent item in Table III is the very narrow peak distribution and coordination number distribution for cooling cycle 1-V. In this particular cooling cycle, 40ps of dynamics at 4000° K was not sufficient to greatly disturb the atoms from their crystalline lattice positions. As a result, the subsequent lower temperature dynamics simulations allowed the atoms to reachieve a very nearly crystalline lattice arrangement. Note, however, that the density is quite low for the cristobalite like structure which was generated. This shows that the MSI glass force field does not well represent this crystalline form of silica.

The data indicate that simulations generated with longer times at high temperatures (8000 K) tend to have the Si-O-Si peak closer to the tetrahedral angle of 109° 27’. For example, the average angle for runs 1-I, 2-I, 2-II, and 2-III which have only 20ps of NVT dynamics at 8000 K was 107.1°. For runs 1-II, 1-III, 1-IV and 2-IV, (40 ps of NVT dynamics at 8000 K) the average was 108°. This difference is statistically significant at the 95% confidence level.

There appears to be a small increase in the Si-O-Si angle with increasing time of 8000 K NVT dynamics but the difference is not statistically significant. The peak breadth also appears to increase slightly with increasing high temperature dynamics but again the difference is not statistically significant.

Looking at the coordination of Si and O, it appears that longer times at temperatures of 2000 K to 4000K are favorable for higher percentages of 4 coordinated Si and 2 coordinated O. The lowest percentages of 4 coordinated Si and 2 coordinated O are found in 1-II where there was a quench directly from 8000 K to 300 K. The highest percentages are found in 2-I and 2-II where 40 ps were spent at 4000 K in addition to significant time spent gradually cooling from 4000 K to 300 K.

The density range over the final 15 ps of the 30 ps NPT runs presented in Table III is 0.02 to 0.03 g/cm^{3} for all of the runs. This means that depending upon the instant that a simulation is terminated, the density of that structure can vary by more than 0.02 to 0.03 g/cm^{3}. However, the difference in average density from one run to the next is even more variable. For example, for cooling cycle 2-V, the difference in average density for the two runs is just under 0.002 g/cm^{3} but for cooling cycle 2-I, the difference is over 0.08 g/cm^{3}. If we look at the lowest and highest density in all of these runs (except 1-V) the density varies from 2.104 g/cm^{3} (1-III) to 2.238 g/cm^{3} (1-II and 2-I), about a 6.5% variation. It is interesting to observe that the density variation is not correlated with either the calculated bond angle peaks or the coordination numbers.

Table IV shows the differences in angles, coordination and density resulting from different force fields and cell size. The data indicate that all three force fields predict the O-Si-O bond angle to be about 107°. However, the MSI force field predicts that the Si-O-Si bond angle is a little over 150° whereas the Tsuneyuki and Demiralp force fields both predict that the Si-O-Si bond angle is in the mid 140° range. It is interesting to note that the larger cell size does not appear to have a significant impact upon the bond angle distributions or the coordination numbers. However, the variation in density over the final 15 ps of NPT dynamics appears to be reduced.

Looking at the coordination numbers, it appears that the Tsuneyuki force field tends to have a higher percentage of 4 coordinated Si and 2 coordinated O than the other two force fields. The number of 2 coordinated O atoms is particularly low for the Demiralp force field.

The average density of the MSI force field gives the lowest density and the density closest to the experimental value of 2.19 g/cm^{3 }while the Tsuneyuki force field gives the highest density. Run to run variation appears to be a little larger for the Tsuneyuki force field than for the other two force fields.

Conclusion

It is apparent from these data that calculated structural factors and densities can be influenced by the cooling cycle used to simulate amorphous silica structures. In addition, it is apparent that even when the same cooling cycle is used, successive runs can have significant structural differences. Therefore, multiple runs must be carried out to determine glass structural features from MD simulations.

References