- Email: [email protected]
Monte Carlo simulations of amorphous hydroxylated silica (SiO2) nanoparticles
Computational Materials Science 135 (2017) 90–98
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Monte Carlo simulations of amorphous hydroxylated silica (SiO2) nanoparticles Naveen Kumar Kaliannan Research Group of Prof. Dr. R. Hentschke, Department of Mathematics and Natural Sciences, Bergische Universität Wuppertal, Wuppertal 42119, Germany
a r t i c l e
i n f o
Article history: Received 11 December 2016 Received in revised form 25 February 2017 Accepted 6 April 2017 Available online 20 April 2017 Keywords: Amorphous hydroxylated silica nanoparticles Bulk amorphous silica Atomic structure Spherical nanoparticles Monte Carlo simulations Silanol group concentration Bulk and surface properties
a b s t r a c t Monte Carlo simulations were carried out on bulk amorphous silica (SiO2) and amorphous hydroxylated silica nanoparticles with different particle sizes (diameter 1–10 nm). The potential developed by D. M. Tether, and modified and tested by Cormack, Du et al. was used to model the interatomic interactions of SiO2 in both cases (bulk and nanoparticles). The bulk system was realized using periodic boundary conditions. The SiO2 nanoparticles were generated by cutting out a sphere from the bulk silica melt to the required particle size. Free valences on the nanoparticle surfaces were saturated via additional hydroxyl groups and then quenched to T = 300 K under free boundary conditions. The effect of particle size on bulk and surface properties of the nanoparticles were calculated at T = 300 K and P = 0 GPa and studied via radial distribution functions, bond angle distributions, bond distances, coordination numbers, different types of OH group concentrations (isolated, geminol and vicinal groups) and radial density profiles. In addition, a comparison study was done between amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica properties in order to study differences between them. The study shows that the bulk properties of amorphous hydroxylated SiO2 nanoparticles are sizedependent and different from those of the bulk silica. However, increasing the particle size leads to an approach of the particle’s bulk properties to the bulk properties of the (quasi) infinite system. In addition, the study also shows that decreasing the particle size (i.e., increase in surface area-to-volume and surface-to-core atom ratios) results in increasing the surface effects and surface OH group concentration. Accordingly, small-sized SiO2 nanoparticles have higher surface OH group concentration and larger surface effects than large-sized SiO2 nanoparticles. The simulation results also show that decreasing the particle size not only results in increasing the surface effects but also results in shifting the bond distances and angles towards higher values. Details of the modeling, simulations results and the study are presented in this paper. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Silica (SiO2) nanoparticles have been an active topic of research in material sciences [1–12] in the recent years due to their unique structural and surface properties, which are different from their bulk. In general, the structural and surface properties of SiO2 nanoparticles are size-dependent in nature. Such size-dependent properties are desirable for many material applications. For instance, silica nanoparticles of size less than 10 nm have higher surface atom concentration and OH group concentration, which exhibit excellent adsorption ability to the toxic pollutants [13]. Furthermore, SiO2 nanoparticles are used as reinforcement fillers in rubber and leather applications due to their high mechanical
E-mail address: [email protected] http://dx.doi.org/10.1016/j.commatsci.2017.04.004 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.
strength, high stability and easily modifiable surface [13]. In biomedical applications, SiO2 nanoparticles are used as a carrier for drug and gene delivery systems [1,12,13]. Because of their unique and excellent properties, silica nanoparticles become an important material in a variety of applications ranging from catalysis [13] to microelectronics [1,7,13]. Thus a better understanding of the size-dependent structural and surface properties of the SiO2 nanoparticles is very important in order to use them in an enhanced way for their applications. The size-dependent structural and surface properties of silica nanoparticles have been studied extensively via various experimental (e.g. NMR, Raman and IR spectroscopy, X-ray and neutron diffraction.) [9,11,14,15], theoretical [9,11] and atomistic computer simulation methods (e.g. Monte Carlo and Molecular Dynamics methods) [1–3,5–8,10]. Experimental studies on silica nanoparti-
91
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
cles using diffraction and spectroscopy have concentrated either on structural [14,15] or on size-dependent surface properties [9]. Computer simulation studies, on the other hand, have focused and provided more detailed and valuable insights into both the properties [1–3,6–8,10]. Also, most of the detailed studies on silica have been conducted using computer simulations. One such detailed simulation study was done by Vo Van Hoang [1,2] in 2007. He studied the structure of spherical amorphous SiO2 nanoparticles of 3 different diameters (2 nm, 4 nm and 6 nm), and found that the structural properties of the nanoparticles are size dependent and different from the bulk. He also found that the core of the nanoparticles is similar to the amorphous bulk structure, while the surface structure is different. Similarly, Roder, Kob and Binder [8] studied the structure of the SiO2 spherical nanoparticles of 3 different diameters (2.65 nm, 3.8 nm and 5.35 nm), and found that the shape of nanoparticles doesn’t depend on the temperature. Their study of the core and the surface structure produces similar results as the one by Vo Van Hoang [1,2]. Du and Cormack [7] simulated the structure and hydroxylation of planar silica glass surfaces. They estimated the surface OH group concentration on planar silica to be 4.54/nm2, which agrees with the experimental value of 4.6/nm2 [9]. The major silica studies have focussed mainly on the non-hydroxylated SiO2 nanoparticles, planar silica surfaces, bulk silica and not on the amorphous hydroxylated SiO2 nanoparticles and spherical silica surfaces. Much less research attention has been given to the amorphous hydroxylated SiO2 nanoparticles, their modeling, their surfaces and their size-dependent properties. In this present work, the focus is on amorphous hydroxylated SiO2 nanoparticles. First, the modeling of amorphous hydroxylated SiO2 nanoparticles is described. Then the effect of particle size on the bulk and surface properties of amorphous hydroxylated SiO2 nanoparticles investigated via computer simulations is discussed. In addition, the difference between the amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica properties simulated at the same temperature and pressure is analyzed and presented in detail.
hand it is easy to implement and overall simpler than the MD technique. The most important part of MC simulation is the potential energy calculations. During the simulations, the potential energy is calculated via interatomic potentials, which are parametrized either on the basis of experimental data or the results of ab initio calculations of smaller systems. The quality of the interatomic potentials decides the validity and quality of the MC simulations. So before starting MC simulations, selecting a reliable interatomic potential that describes the system and its properties is very important. Because of the technological importance of silica, several interatomic potentials have been developed to describe the various SiO2 structures and to study their properties. In particular, a potential proposed by Beest, Kramer and Van Saten (BKS) [25] has been widely used for studying the bulk properties and structure of silica glass. Recently, D. M. Teter, Cormack and Du et. al. (TCD) [7,26] have proposed a set of potentials for common oxides that include most of the silica glass interactions, especially the surface hydroxylation. In this work, we use the TCD potential because it has interatomic potentials for both bulk silica and amorphous hydroxylated SiO2 nanoparticles. Like the BKS potential, the TCD potential has a short-range Buckingham and a long-range Coulomb potential. The equation for the pair potential energy is given by the following expression:
U Pair1 ðr ij Þ ¼
1
4p0
qi qj e2 C ij þ Aij eBij rij 6 r ij rij
ð1Þ
where, i, j 2 {Si, O, H}. U Pair1 ðr ij Þ is the potential energy between atoms i and j, rij is the distance between these atoms, qi and qj are their partial charges, and Aij, Bij and Cij are constants, 0 is the permittivity of free space and e is the elementary charge. In order to simulate the OH group on the silica surfaces, i.e. the silanol group (Si-O-H), with a correct bond angle, TCD have introduced a Coulomb-subtracted Morse potential for hydroxyl interactions with a screened harmonic three-body angular potential for silanol (Si-O-H) bond angles [7,26]. The Coulomb-subtracted Morse potential equation for the hydroxyl groups is given by the following expression:
1 qH qO e 2 4p0 r OH
2. Computational details
U Pair2 ðr OH Þ ¼ Dð1 exp½bðr r 0 ÞÞ2 D
2.1. Simulation method
where D, b, r0 are constants, qH and qO are the atomic partial charges of hydrogen and oxygen in the hydroxyl group and rOH is their distance. The three-body angular potential for the silanol group (Si-OH) has he screened harmonic form,
In the past the Molecular Dynamics (MD) and Monte Carlo (MC) simulation techniques have been successfully applied to the simulation of equilibrium properties of both bulk silica and silica nanoparticles. In MD simulations [16,17] the equilibrium and transport properties are calculated from the trajectory of all atoms, i.e., from the time evolution of positions and momentum. The time evolution is calculated by solving Newton’s equation of motion, using effective interatomic potentials. The MD simulations provide structural, thermodynamic and dynamic information. There have been numerous MD simulation studies of silica in the past [1,2,6–8,18–23]. MC simulations [24], on the other hand, generate a large sequence of configurations using random numbers. Each configuration is accepted or rejected according to certain probability conditions (like the so-called Metropolis criterion) that depend mainly on the potential energy differences between configurations and the temperature. The expectation values of the properties of interest are determined by averaging over the accepted configurations. The MC simulations provide only the structural and thermodynamic information. Both these techniques have found widespread use in the field of molecular simulations. In this work, we focus on properties which are not explicitly dependent on time. Thus we employ the MC technique, because for the purpose at
U Angular ðhjik ; rik ; r ij Þ ¼
2 E0 r ij rik þ : hjik h0 exp 2 q1 q2
ð2Þ
ð3Þ
where rij and rik are the distances between the center atom and its neighbors. hijk is the angle between j, i and k. Here i represent the center atom and j and k are the neighbor atoms. In Si-O-H the oxygen is the center atom and silicon and hydrogen are the corresponding neighbor atoms. E0 ; q1 and q2 are the constants. The total potential energy U total of a SiO2 particle is given by the following expression:
U total ¼
X X X U Pair1 ðr ij Þ þ U Pair2 ðr ij Þ þ U Angular ðhjik ; r ik ; r ij Þ i
i
ð4Þ
i
TCD have fitted their potential constants to correctly describe the bulk properties and structure of silica over a wide temperature and pressure range [7,26]. But the potential constants for OH group and silica-OH interactions are fitted in a crystalline system based on the structure and mechanical properties under ambient temperature and pressure [7,26]. Du and Cormack [7] have tested these
92
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
fitted potential constants and proved they are reliable to simulate the hydroxylation on silica surfaces. The fitted potential constants and charges are listed in the Tables 1–4. The bulk silica has only two different ions or partially charged atoms (Si, O) and three different pair interactions (Si-Si, O-O, SiO) between them (see Table 6). This makes the bulk silica easy to model. But unlike the bulk silica, the hydroxylated silica nanoparticles have four different ions (Si, O (in silica), O (in OH groups), H), one angular (Si-O-H) and ten different pair interactions, which make them more complex to model. Because the hydroxylated silica nanoparticles have so many interactions, the potential functions for their corresponding pair and angular interaction must be clearly described and understood before modeling. Here we do this by tabulating their potential functions in a table (see Table 5). In addition, we also show the pair potential energy for the corresponding pair interactions in Fig. 1 using the potential constants and charges in Tables 1–4.
Table 2 Buckingham potential constants of the TCD model [26,7].
2.2. Modeling of hydroxylated silica nanoparticles
Table 4 Three body angular potential constants of the TCD model [26,7].
First, the bulk silica melt was prepared by heating a random structure to T = 4000 K and P = 0 GPa via NPT-ensemble based Monte Carlo simulation using periodic boundary conditions. Here the simulations use TCD potential as described in the previous Section 2.1. Coulomb interactions are calculated using the Ewald summation [28–31]. After generating the bulk silica melt, the SiO2 nanoparticles were generated by cutting out a sphere from the melt structure. Free valences on the nanoparticle surfaces were saturated via additional hydroxyl groups as shown in Fig. 3. Then the saturated SiO2 nanoparticles were quenched to T = 300 K under free boundary conditions. Since the simulation uses free boundary conditions (non-periodic), Ewald summation is no more applicable for calculating the Coulomb interactions. Instead a double loop over all atom pairs for the Coulomb interactions was implemented, as shown in the following Eq. (5).
U Coulomb ¼
e2 4p0
N1 X N X
qi qj r i¼1 j¼iþ1 ij
ð5Þ
Bulk simulation: The bulk simulation was performed with 100SiO2 (100 silicon and 200 oxygen atoms) groups at T = 300 K and P = 0 GPa. The simulation consisted of an initial equilibration of 2.5 million moves followed by additional 2.5 million trail moves (see Fig. 2(a)). The average values of all the properties of interests were calculated by averaging over the final 2.5 million trail moves. From the bulk simulations, the averages of the partial and total radial distribution functions, bond angle distributions, bond distances and coordination number were calculated. Particle simulations: The particle simulations were performed with ten different sized nanoparticles (1–10 nm) at T = 300 K and P = 0 GPa. Each simulation consisted of an initial equilibration of 10 million moves followed by additional 10 million trail moves (see Fig. 2(b)). The average values of all the properties were calculated by averaging over the final 10 million trail moves. From the particle simulations, the averages of the density profile, radial dis-
Table 1 Atomic partial charges of the TCD model [26,7]. Atom
Charges
qSi qO (in silica glass) qO (in hydroxyl group) qH
+2.4 1.2 0.856 +0.256
i-j
Aij (eV)
Bij (Å1)
Cij (eV Å6)
Si2.4-O1.2 O1.2-O1.2 Si2.4-O0.856 O0.856-O0.856 O1.2-O0.856 H0.256-O0.856 H0.256-O1.2
13702.905 1844.7458 12443.824 1844.7458 1844.7458 100.0 100.0
5.1595 2.9099 5.1595 2.9099 2.9099 4.0 4.0
54.681 192.58 54.681 192.58 192.58 0.0 0.0
Table 3 Coulomb-subtracted Morse potential constants of the TCD model [26,7]. i-j
D (eV)
b (Å1)
r0 (Å)
H0.256-O0.856 H0.256-O1.2
7.0525 7.0525
5.26 5.26
0.9485 0.9485
i-j-k
E0 (eV)
h0 (deg)
q1 (Å)
q2 (Å)
H0.256-O0.856-Si2.4
12.0
118.5
2.0
3.2
Table 5 The potential functions corresponding to each pair and angular interactions in the SiO2 nanoparticles. The angular interaction is given in the following table. Pair (i-j)
Potential function
Si2.4-Si2.4 Si2.4-O1.2 Si2.4-O0.856 Si2.4-H0.256 O1.2-O1.2 O1.2-O0.856 O1.2-H0.256 O0.856-O0.856 O0.856-H0.256 H0.256-H0.256
UCoulomb UCoulomb + UBuckingham UCoulomb + UBuckingham UCoulomb UCoulomb + UBuckingham UCoulomb + UBuckingham UCoulomb + UBuckingham + UCoulombsubrt. UCoulomb + UBuckingham UCoulomb + UBuckingham + UCoulombsubrt. UCoulomb
Morse
Morse
Angular (j-i-k)
Potential function
H0.256-O0.856-Si2.4
UAngular
Table 6 The potential functions corresponding to each pair interactions in bulk silica model. Pair (i-j)
Potential function
Si2.4-Si2.4 Si2.4-O1.2 O1.2-O1.2
UCoulomb UCoulomb + UBuckingham UCoulomb + UBuckingham
tribution functions, types of OH groups, bond angle distributions, bond distances and coordination number were calculated.
3. Simulation results and discussion The type of properties investigated here from the simulations are structural and surface properties. First, we study the structural properties for different sizes of amorphous hydroxylated SiO2 nanoparticles at the temperature of 300 K and zero pressure:
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
93
Fig. 1. The pair potential energy of Si-Si, O-O, Si-O, H-H, Si-H and O-H as the function of the interatomic distance r. At large r the energy is very small. When r decreases, the attractive Si-O, O-H interactions and the repulsive Si-Si, H-H, O-O interactions increase. The repulsion region between Si-O is very small. It prevents the Si-O atoms from coming too close during the simulations.
(a) Potential energy per SiO2 group in bulk SiO2 , bulk density vs. MC-steps.
(b) Total potential energy (T. P. E) of the SiO2 nanoparticles vs. MC-steps
Fig. 2. The density and potential energy profiles vs. Monte Carlo steps. Both the profiles have the initial equilibration region (blue line) and the equilibrium data region (red line) used for study and analysis. The expectation values of the properties of interest are determined by averaging over the equilibrium data region. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1. Structural properties The radial density profile for different sizes of hydroxylated SiO2 spherical nanoparticles as a function of r from the center of the particles is shown in Fig. 4. The radial density profile shown in Fig. 4 is determined as follows: we find the average number of atoms in the spherical shell of radius r with the thickness of 2 Å from the center of particles. Then, we compute the radial density profile by dividing the average number of atoms by the volume of each spherical shell. Mathematically the radial atom number density profile qðrÞ is expressed as:
qðrÞ ¼
Nðr; drÞ 4pr2 dr
ð6Þ
where Nðr; drÞ is the number of atoms in a shell of radius r from the center of particle and thickness dr. 4pr 2 dr is the shell volume. All the radial density profiles show similar atom number density fluctuations around 0.066/Å3, i.e. 0.044/Å3 for oxygen and 0.022/Å3 for silicon, respectively. At smaller radius r the density
fluctuation qðrÞ is very random and noisy because of the less number of atoms. As the r becomes larger, the fluctuation qðrÞ becomes smoother. So we consider only the larger spherical shells for our analysis. The radial density profile of our amorphous nanoparticles is compared to the one of the crystalline systems with the same size to study the difference between them. Clearly the density fluctuation qðrÞ is very stronger for crystalline nanoparticles, indicates that the crystalline structure does not have uniform density in all the spherical shells. But the atomic disorder makes the density uniform throughout the system, that is why the weaker fluctuation, smaller peak heights, widths and depths in amorphous SiO2 nanoparticles (see the region between 20 and 30 Å in 7 nm particle). In the density profiles, clearly the density fluctuation of silicon reaches zero first, then the oxygen reaches, followed by hydrogen. Also, the density fluctuation of hydrogen is observed only in the surface regions. This means our simulations successfully simulated the OH group on the silica surfaces. We also observe a sudden drop of the radial density near the surface. The surface region starts
94
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98 Table 7 The surface thickness for different sizes of the amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa.
Fig. 3. Different possibilities of bond breaking and OH group saturation on the spherical silica surfaces. Trisilanol groups and isolated oxygens are removed from the surfaces because they both were not part of silica surfaces. The trisilanol group is removed by removing the corresponding silicon and 3 oxygen atoms [6]. The red line is the surface cut. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) 4.5nm
Particle diameter (nm)
Surface thickness (nm)
2 3 4 5 7 9 10 Other simulation studies
0.63 0.62 0.64 0.63 0.62 0.66 0.61 0.5–1 [7]
from the point where the overall number density suddenly falls, and reaches to zero, which is the end point of surface region. With these two points, starting and ending point of the surface region, the surface thickness of the nanoparticles can be easily estimated. The estimated surface thickness values are shown in Table 7. The surface thickness values almost remain constant, around 0.6 nm, and are unaffected by a change in particle size. Other simulation studies [7] have suggested values between 0.5–1 nm. In order to study amorphous hydroxylated SiO2 nanoparticles structures and their variation with particle size, we also show the partial and total RDFs of hydroxylated SiO2 spherical nanoparticles for 3 different sizes (2 nm, 3 nm and 7 nm) in Fig. 6. Clearly the first (see magnified plot) and second peaks in the RDFs of the bulk and nanoparticles are broader, indicates that the bulk and nanoparticle structures are disordered (i.e., amorphous). This suggests that the melted structure quenched to 300 K via MC simulations leads to an amorphous state. It is also found that the peak heights in the RDFs of the nanoparticles change with the particle size and are different from the bulk. The peak heights lower with decrease in the particle size because of an increase in the surface contribution and surface effects. When the particle size increases, the peak heights approach towards the bulk, indicates that the particle structure approaches the bulk structure as the particle size increases. This RDFs behavior with particle size is similar to the behavior observed in Vo Van Hoang study [1,2]. The bond distances calculated from the RDFs are shown in Table 8. The particle’s bond distances show excellent agreement with the bulk values, except for 2 nm particle, possibly resulting from the larger surface effects. The Si-Si, O-O and Si-O coordination numbers calculated from the bulk silica and amorphous hydroxylated SiO2 nanoparticles
(b) 7nm
Fig. 4. The radial number density profile of hydroxylated SiO2 nanoparticles, for both crystalline (top) and amorphous (bottom). Note: The crystalline hydroxylated SiO2 nanoparticles were prepared by using the crystal structure of SiO2 (beta-cristoabalite).
95
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98 Table 8 The first and second peak positions of the partial and total RDFs for amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica at T = 300 K and P = 0 GPa. Pair
Si-Si O-O Si-O H-O Si-H
Peak
Particle size (nm)
(Å)
2
3
7
First Second First Second First Second First First
3.16 5.17 2.61 5.14 1.60 4.00 0.98 2.2
3.14 5.10 2.61 5.08 1.60 4.05 0.98 2.2
3.13 5.10 2.61 5.08 1.60 4.06 0.98 2.2
Bulk
EXP. [15,14]
3.13 5.15 2.61 5.06 1.60 4.13 – –
3.13 5.18 2.65 4.95 1.62 4.15 0.95–1 –
Table 9 The coordination number (CN) of Si-Si, O-O and Si-O for amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica at T = 300 K and P = 0 GPa. Pair i-j
CN
Si2.4-Si2.4 Si2.4-O1.2 = 0.856 Si2.4-O1.2 O1.2 = 0.856-O1.2 = 0.856 O1.2-O1.2
(a) Si-O-Si bond angle
First First First First First
Particle size (nm)
Bulk
2
3
7
3.08 4.00 3.08 4.90 4.6
3.40 4.00 3.40 5.28 5.15
3.73 4.00 3.73 5.67 5.65
(b) O-Si-O bond angle
4.00 4.00 4.00 6.00 6.00
(c) Si-O-H bond angle
Fig. 5. The bond angle distributions (BAD) for bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa.
structures are shown in Table 9. The Si-O, O-O and Si-Si coordination numbers of bulk silica are found to be 4, 6 and 4. The Si-O coordination number equal to four means that each silicon atom is surrounded by four oxygen atoms. Therefore the basic building block of bulk silica is SiO4 tetrahedron. The coordination numbers of the nanoparticles, on the other hand, change with the particle size, and approach the bulk values as the particle size increases. However this change is not for Si2.4 - O1.2 =0.875. Bond angle distributions have been used to analyze the local structure in the silica. Si-O-Si and O-Si-O are the two major bond angles normally investigated and studied in silica structure. Here the O-Si-O, Si-O-Si and Si-O-H bond angles are investigated to analyze the local structure in the core and surface of the nanoparticles. Fig. 5 shows the investigated bond angle distributions for bulk silica and amorphous hydroxylated SiO2 nanoparticles structures at T = 300 K and P = 0 GPa. Although there is a clear short range ordering characterized by the peaks in the RDFs of the bulk and nanoparticles, the Si-O-Si and O-Si-O bond angle distributions are significantly very broad confirming the disordered structure in both cases (bulk and nanoparticles). Like the RDFs, the Si-O-Si and O-Si-O bond angle distributions of the nanoparticles change with particle size, and are slightly different from the bulk. Also,
the Si-O-Si and O-Si-O BADs of the nanoparticles approach the bulk distribution as the particle size increases. In the Si-O-Si BAD case, there is even a clear size effect evolution with a shift toward higher angles (>160 ) with decreasing particle size (see Table 10). This means the connectivity between SiO4 tetrahedral structures changes significantly as the particle size decreases. On the other hand, the peak maximum in the O-Si-O BAD of the bulk and nanoparticles is located at 109.0 , which agrees with the experimental values (109.47 , 109.50 ). It means that the surface effect present in the nanoparticles doesn’t affects the local structure inside the SiO4 tetrahedral structures. The Si-O-H bond angle is narrowly distributed at 118.5 with a deviation of 0 , which means all the Si-O-H bond angles are reproduced exactly at 118.5 . This value agrees with the experimental and theoretical value of 118.5 , found by McCarthy, Tamassia, Woon, and Thaddeus [11]. The comparison between simulated and experimental bond angles is shown in Table 10. 3.2. Surface properties In most cases, the surface properties of silica depends mainly on the surface hydroxyl group concentration and distribution [9]. The
96
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98 Table 10 Comparison of simulated and experimental bond angles for amorphous hydroxylated SiO2 nanoparticles and bulk silica at T = 300 K and P = 0 GPa. Angle
Si-O-Si O-Si-O Si-O-H
Fig. 6. The total and partial radial distribution functions for bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa. The peak heights in the Si-O (increase), O-O (increase), and O-H (decrease) RDFs of the nanoparticles approach the bulk peak heights as the particle size increases, while the first peak height in the Si-Si RDF of the nanoparticles moves away from the bulk. This RDFs behavior is similar to the one observed in Vo Van Hoang study [1,2].
surface OH group concentration and distribution on the silica determines whether the silica surface is hydrophobic or hydrophilic [9]. For instance, the silica surface with higher surface OH group concentration produces more hydrophilic effect. In contrast the surface with less surface OH group concentration produces
Particle size (nm) 2
3
7
161 109 118.5
160 109 118.5
156 109 118.5
Bulk
EXP. [14,11,15]
155 109 –
144 , 154 109.47–.50 118.5
Fig. 7. (a) Proportion of atoms. (b) Surface area-to-volume ratio, surface-to-core atom ratio and oxygen-to-silicon atom ratio profiles. (c) Surface and core atom concentration of the nanoparticles.
hydrophobic effect. For a better understanding of the surface properties of silica nanoparticles, some investigation is needed into OH group concentration. The OH group concentration investigated from the simulated hydroxylated SiO2 spherical nanoparticles is
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
Fig. 8. The surface OH group concentration on SiO2 spherical nanoparticles for different particle sizes.
shown in Fig. 8. Clearly the OH group concentration decreases as the particle size increases. Therefore, small-sized nanoparticles have higher surface OH group concentration than large-sized nanoparticles. The average surface OH group concentration on SiO2 spherical nanoparticles is estimated to be around 7.4/nm2, which is very close to the theoretical value (7.8/nm2), found by IIler [27]. Other experimental and simulation studies have suggested different OH group concentrations between 4.5 and 9/nm2 [6,7,9]. These difference are mainly due to the distinct geometry, size, initial conditions and potentials used in the simulations. The surface area-to-volume ratio, oxygen-to-silicon atom ratio, surface-to-core atom concentration and surface-to-core atom ratio profiles of the nanoparticles are shown in Fig. 7. Clearly the surface area-to-volume and surface-to-core atom ratio’s increase as the particle size decreases. Accordingly, small-sized nanoparticles have higher surface area-to-volume ratio, higher surface atom concentration and higher surface-to-core atom ratio than large-sized nanoparticles.
4. Summary In this paper, the structural and surface properties of the amorphous hydroxylated SiO2 nanoparticles with different sizes (1– 10 nm) have been studied at T = 300 K and P = 0 GPa via MC simulations using the TCD model potential. In addition, the comparison between the bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles properties simulated at the same temperature and pressure has been analyzed. The most important conclusions drawn from the study and analysis are the followings: 1. The study shows that the TCD potential is able to successfully reproduce the structure of bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa. The size-dependent structural and surface properties were determined and were shown to be in good agreement with other experimental and simulation studies. From the agreement, it can be suggested that the TCD model potential and the methodology described in this paper are suitable and reliable for the modeling and simulation of amorphous hydroxylated SiO2 nanoparticles. 2. In this work, the SiO2 nanoparticles were generated by cutting out a sphere from the bulk silica melt to the required particle size. Free valences on the nanoparticle surfaces were saturated via additional hydroxyl groups and then quenched to T = 300 K under free boundary conditions via MC simulations. The MC simulation results clearly show that the bulk and nanoparticle structures are disordered i.e., amorphous. This means that the melted structure
97
quenched to 300 K via MC simulations leads to an amorphous state. Therefore, MC simulation using the TCD model potential does not succeed to crystallise the bulk SiO2 and SiO2 nanoparticles at 300 K. 3. The calculations of the structural properties show that the structural properties of the nanoparticles are size-dependent and different from those of bulk silica. As expected, increasing the particle size leads to an approach of the particle’s bulk properties to the bulk properties of the (quasi) infinite system. The calculations of the density profile of the nanoparticles show that the core structure and surface thickness of the nanoparticles are unaffected by a change in particle size. 4. The simulation results suggest that there is a clear size effects on the bond angle distributions and radial distribution functions. In the Si-O-Si BAD and Si-Si RDF cases, there is even a clear size effect evolution with a shift toward higher angles and bond distances with decrease in particle size. This means the decreasing in particle size not only results in increasing the surface effects but also results in shifting the bond distances and angles towards higher values. 5. The calculations of the surface properties of the nanoparticles show that decreasing the particle size results in increasing the surface OH group concentration. Accordingly, small-sized SiO2 nanoparticles have higher surface OH group concentration than large-sized SiO2 nanoparticles. The average surface OH group concentration on the SiO2 spherical nanoparticles is estimated to be around 7.4/nm2, which is very close to the theoretical value of 7.8/nm2 found by Iler [27]. Other experimental and simulation studies have suggested different OH group concentrations between 4.5 and 9/nm2 [6,7,9]. These difference are mainly due to the distinct geometry, size, initial conditions and potentials used in the simulations.
Acknowledgements I would like to thank Prof. Dr. Reinhard Hentschke (Email: [email protected]) for giving me a chance to write this paper in his research group (http://constanze.materials.uni-wuppertal.de) and for his encouragement and support. During this work, I have learned a lot about silica modeling and simulations from his great experience and knowledge. Again I would like to thank him for proofreading this paper and for his comments and corrections. Finally I would like to thank my family and friends for their support and encouragement.
References [1] V.V. Hoang, Molecular dynamics simulation of amorphous SiO2 nanoparticles, J. Phys. Chem. 111 (2007) 44. [2] N.T. Huynh, V.V. Hoang, H. Zung, Structural properties of amorphous SiO2 nanoparticles, 2016.. [3] Naveen Kumar Kaliannan, A study on the structure and properties of silica glass and silica nanoparticles via Monte Carlo simulations, Master thesis, University of Wuppertal, Germany, 2016, submitted for publication, https:// www.researchgate.net/publication/311138439_A_study_on_the_structure_ and_properties_of_silica_glass_and_silica_nanoparticles_via_Monte_ Carlosimulations, doi: http://dx.doi.org/10.13140/RG.2.2.34922.11202. [4] C. Silvina, A.S. Gustov, O. Jon, W. Stephen, Dielectric study of hydration water in silica nanoparticles, J. Phys. Chem. 116 (45) (2012) 24340–24349. [5] D. Makimura, C. Metin, T. Kabashima, T. Matsuoka, Q.P. Nguyen, Caetano, R. Miranda, Combined modeling and experimental studies of hydroxylated silica nanoparticles, J. Mater. Sci. 45 (2010) 5084–5088. [6] J.M. Stallons, E. Iglesia, Simulations of the structure and properties of amorphous silica surfaces, Chem. Eng. Sci. 56 (2001) 4205–4216. [7] J. Du, A.N. Cormack, Molecular dynamics simulation of the structure and hydroxylation of silica glass surfaces, J. Am. Ceram. Soc. 88 (9) (2005) 2532– 2539. [8] A. Roder, W. Rob, K. Binder, Structure and dynamics of amorphous silica surfaces, J. Chem. Phys. 114 (2001) 7602.
98
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
[9] L.T. Zhuravlev, The surface chemistry of amorphous silica, Zhuravlev Model Colloids Surfaces A: Physiochem. Eng. Aspects 173 (1–3) (2000) 1–38. [10] I.V. Schweigert, K.E.J. Lehtinen, M.J. Carrier, M.R. Zachariah, Structure and properties of silica nanoclusters at high temperatures, Phys. Rev. B 65 (2002) 235410. [11] M.C. McCarthy, F. Tamassia, D.E. Woon, P. Thaddeus, A laboratory and theoretical study of silicon hydroxide Si-O-H, J. Chem. Phys. 129 (18) (2008) 184301. [12] D. Douroumis, I. Onyesom, M. Maniruzzaman, J. Mitchell, Mesoporous silica nanoparticles in nanotechnology, Crit. Rev. Biotechnol. 33 (3) (2012). [13] R. Nandanwar, P. Singh, F. zia Haque, Synthesis and properties of silica nanoparticles by sol-gel method for the application of green chemistry, Mater. Sci. India 10 (1) (2013) 85–92. [14] R.L. Mozzi, B.E. Warren, Structure of vitreous silica, J. Appl. Crystallogr. 2 (4) (1969) 164–172. [15] J.R.G. DaSilva, D.G. Pinatti, C.E. Anderson, M.L. Rudee, A refinement of the structure of vitreous silica, Philos. Mag. Part B 31 (3) (1974) 713–717. [16] R. Khanna, V. Sahajwalla, Treatise on process metallurgy, Process Fundamentals, first ed., vol. 1, Elsevier, 2013. [17] D. Beljonne, J. Cornil, Multiscale Modelling of Organic and Hybrid Photovoltaics, Springer, Berlin, 2014. [18] J. Horbach, Molecular dynamics computer simulation of amorphous silica under high pressure, J. Phys.: Condensed Matter 20 (2008) 24. [19] A. Carre, Development of empirical potentials for amorphous silica, Doctoral thesis, Deparment of Physics, Johannes Gutenberg University mainz, September 2007.
[20] K. Glaser, Computational studies of silica, Doctoral thesis, Department of Chemistry, University College London, May 2011. [21] D. Herzbach, Comparison of model potentials for molecular dynamics simulation of crystalline silica, Doctoral thesis, Johannes Gutenberg, Universitat Mainz, July 2004. [22] J.D. Kubicki, A.C. Lasaga, Molecular dynamics simulations of SiO2 melt and glass: Ionic and covalent models, Am. Mineral. 73 (1988) 941–955. [23] P. Tangney, S. Scandolo, An ab initio parametrized interatomic force field for silica, J. Chem. Phys. 117 (2002) 19. [24] F. Rouquerol, J. Rouquerol, K. Sing, Adsorption by Powders and Porous Solids Principles, Methodology and Applications, Academic Press, San Diego, 1998. [25] V. Beest, G. Kramer, V. Santen, Force fields for silicas and aluminophosphates based on ab initio calculations, Phys. Rev. Lett. 64 (1990) 1955. [26] C. Massobrio, J. Du, M. Bernasconi, P.S. Salmon, Molecular Dynamics Simulations of Disordered Materials, Springer Series in Material Sciences (2015) 215. [27] Ralph K. Iler, The Chemistry of Silica: Solubility, Polymerization, Colloid and Surface Properties and Biochemistry of Silica, John Wiley and Sons, New York, 1979. [28] Leeuw, Perram, Smith, Simulation of electrostatic systems in PBC. I. Lattice sums and dielectric constants, Roy. Soc. Lond. Ser. A 373 (1980) 27–56. [29] B.A. Wells, A.L. Chaffee Ewald, Summation for molecular simulations, J. Chem. Theory Comput. 11 (2015) 3684–3695. [30] U. Essmann, L. Perera, Berkowitz, H. Lee, L. Pedersen, A smooth particle Mesh Ewald potential, J. Chem. Phys. 103 (1995) 8577–8592. [31] T. Darden, D. York, L. Pedersen, Particle mesh Ewald: a N.log(N) method for Ewald sums in large systems, J. Chem. Phys. 98 (12) (1993) 10089.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Monte Carlo simulations of amorphous hydroxylated silica (SiO2) nanoparticles Naveen Kumar Kaliannan Research Group of Prof. Dr. R. Hentschke, Department of Mathematics and Natural Sciences, Bergische Universität Wuppertal, Wuppertal 42119, Germany
a r t i c l e
i n f o
Article history: Received 11 December 2016 Received in revised form 25 February 2017 Accepted 6 April 2017 Available online 20 April 2017 Keywords: Amorphous hydroxylated silica nanoparticles Bulk amorphous silica Atomic structure Spherical nanoparticles Monte Carlo simulations Silanol group concentration Bulk and surface properties
a b s t r a c t Monte Carlo simulations were carried out on bulk amorphous silica (SiO2) and amorphous hydroxylated silica nanoparticles with different particle sizes (diameter 1–10 nm). The potential developed by D. M. Tether, and modified and tested by Cormack, Du et al. was used to model the interatomic interactions of SiO2 in both cases (bulk and nanoparticles). The bulk system was realized using periodic boundary conditions. The SiO2 nanoparticles were generated by cutting out a sphere from the bulk silica melt to the required particle size. Free valences on the nanoparticle surfaces were saturated via additional hydroxyl groups and then quenched to T = 300 K under free boundary conditions. The effect of particle size on bulk and surface properties of the nanoparticles were calculated at T = 300 K and P = 0 GPa and studied via radial distribution functions, bond angle distributions, bond distances, coordination numbers, different types of OH group concentrations (isolated, geminol and vicinal groups) and radial density profiles. In addition, a comparison study was done between amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica properties in order to study differences between them. The study shows that the bulk properties of amorphous hydroxylated SiO2 nanoparticles are sizedependent and different from those of the bulk silica. However, increasing the particle size leads to an approach of the particle’s bulk properties to the bulk properties of the (quasi) infinite system. In addition, the study also shows that decreasing the particle size (i.e., increase in surface area-to-volume and surface-to-core atom ratios) results in increasing the surface effects and surface OH group concentration. Accordingly, small-sized SiO2 nanoparticles have higher surface OH group concentration and larger surface effects than large-sized SiO2 nanoparticles. The simulation results also show that decreasing the particle size not only results in increasing the surface effects but also results in shifting the bond distances and angles towards higher values. Details of the modeling, simulations results and the study are presented in this paper. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Silica (SiO2) nanoparticles have been an active topic of research in material sciences [1–12] in the recent years due to their unique structural and surface properties, which are different from their bulk. In general, the structural and surface properties of SiO2 nanoparticles are size-dependent in nature. Such size-dependent properties are desirable for many material applications. For instance, silica nanoparticles of size less than 10 nm have higher surface atom concentration and OH group concentration, which exhibit excellent adsorption ability to the toxic pollutants [13]. Furthermore, SiO2 nanoparticles are used as reinforcement fillers in rubber and leather applications due to their high mechanical
E-mail address: [email protected] http://dx.doi.org/10.1016/j.commatsci.2017.04.004 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.
strength, high stability and easily modifiable surface [13]. In biomedical applications, SiO2 nanoparticles are used as a carrier for drug and gene delivery systems [1,12,13]. Because of their unique and excellent properties, silica nanoparticles become an important material in a variety of applications ranging from catalysis [13] to microelectronics [1,7,13]. Thus a better understanding of the size-dependent structural and surface properties of the SiO2 nanoparticles is very important in order to use them in an enhanced way for their applications. The size-dependent structural and surface properties of silica nanoparticles have been studied extensively via various experimental (e.g. NMR, Raman and IR spectroscopy, X-ray and neutron diffraction.) [9,11,14,15], theoretical [9,11] and atomistic computer simulation methods (e.g. Monte Carlo and Molecular Dynamics methods) [1–3,5–8,10]. Experimental studies on silica nanoparti-
91
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
cles using diffraction and spectroscopy have concentrated either on structural [14,15] or on size-dependent surface properties [9]. Computer simulation studies, on the other hand, have focused and provided more detailed and valuable insights into both the properties [1–3,6–8,10]. Also, most of the detailed studies on silica have been conducted using computer simulations. One such detailed simulation study was done by Vo Van Hoang [1,2] in 2007. He studied the structure of spherical amorphous SiO2 nanoparticles of 3 different diameters (2 nm, 4 nm and 6 nm), and found that the structural properties of the nanoparticles are size dependent and different from the bulk. He also found that the core of the nanoparticles is similar to the amorphous bulk structure, while the surface structure is different. Similarly, Roder, Kob and Binder [8] studied the structure of the SiO2 spherical nanoparticles of 3 different diameters (2.65 nm, 3.8 nm and 5.35 nm), and found that the shape of nanoparticles doesn’t depend on the temperature. Their study of the core and the surface structure produces similar results as the one by Vo Van Hoang [1,2]. Du and Cormack [7] simulated the structure and hydroxylation of planar silica glass surfaces. They estimated the surface OH group concentration on planar silica to be 4.54/nm2, which agrees with the experimental value of 4.6/nm2 [9]. The major silica studies have focussed mainly on the non-hydroxylated SiO2 nanoparticles, planar silica surfaces, bulk silica and not on the amorphous hydroxylated SiO2 nanoparticles and spherical silica surfaces. Much less research attention has been given to the amorphous hydroxylated SiO2 nanoparticles, their modeling, their surfaces and their size-dependent properties. In this present work, the focus is on amorphous hydroxylated SiO2 nanoparticles. First, the modeling of amorphous hydroxylated SiO2 nanoparticles is described. Then the effect of particle size on the bulk and surface properties of amorphous hydroxylated SiO2 nanoparticles investigated via computer simulations is discussed. In addition, the difference between the amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica properties simulated at the same temperature and pressure is analyzed and presented in detail.
hand it is easy to implement and overall simpler than the MD technique. The most important part of MC simulation is the potential energy calculations. During the simulations, the potential energy is calculated via interatomic potentials, which are parametrized either on the basis of experimental data or the results of ab initio calculations of smaller systems. The quality of the interatomic potentials decides the validity and quality of the MC simulations. So before starting MC simulations, selecting a reliable interatomic potential that describes the system and its properties is very important. Because of the technological importance of silica, several interatomic potentials have been developed to describe the various SiO2 structures and to study their properties. In particular, a potential proposed by Beest, Kramer and Van Saten (BKS) [25] has been widely used for studying the bulk properties and structure of silica glass. Recently, D. M. Teter, Cormack and Du et. al. (TCD) [7,26] have proposed a set of potentials for common oxides that include most of the silica glass interactions, especially the surface hydroxylation. In this work, we use the TCD potential because it has interatomic potentials for both bulk silica and amorphous hydroxylated SiO2 nanoparticles. Like the BKS potential, the TCD potential has a short-range Buckingham and a long-range Coulomb potential. The equation for the pair potential energy is given by the following expression:
U Pair1 ðr ij Þ ¼
1
4p0
qi qj e2 C ij þ Aij eBij rij 6 r ij rij
ð1Þ
where, i, j 2 {Si, O, H}. U Pair1 ðr ij Þ is the potential energy between atoms i and j, rij is the distance between these atoms, qi and qj are their partial charges, and Aij, Bij and Cij are constants, 0 is the permittivity of free space and e is the elementary charge. In order to simulate the OH group on the silica surfaces, i.e. the silanol group (Si-O-H), with a correct bond angle, TCD have introduced a Coulomb-subtracted Morse potential for hydroxyl interactions with a screened harmonic three-body angular potential for silanol (Si-O-H) bond angles [7,26]. The Coulomb-subtracted Morse potential equation for the hydroxyl groups is given by the following expression:
1 qH qO e 2 4p0 r OH
2. Computational details
U Pair2 ðr OH Þ ¼ Dð1 exp½bðr r 0 ÞÞ2 D
2.1. Simulation method
where D, b, r0 are constants, qH and qO are the atomic partial charges of hydrogen and oxygen in the hydroxyl group and rOH is their distance. The three-body angular potential for the silanol group (Si-OH) has he screened harmonic form,
In the past the Molecular Dynamics (MD) and Monte Carlo (MC) simulation techniques have been successfully applied to the simulation of equilibrium properties of both bulk silica and silica nanoparticles. In MD simulations [16,17] the equilibrium and transport properties are calculated from the trajectory of all atoms, i.e., from the time evolution of positions and momentum. The time evolution is calculated by solving Newton’s equation of motion, using effective interatomic potentials. The MD simulations provide structural, thermodynamic and dynamic information. There have been numerous MD simulation studies of silica in the past [1,2,6–8,18–23]. MC simulations [24], on the other hand, generate a large sequence of configurations using random numbers. Each configuration is accepted or rejected according to certain probability conditions (like the so-called Metropolis criterion) that depend mainly on the potential energy differences between configurations and the temperature. The expectation values of the properties of interest are determined by averaging over the accepted configurations. The MC simulations provide only the structural and thermodynamic information. Both these techniques have found widespread use in the field of molecular simulations. In this work, we focus on properties which are not explicitly dependent on time. Thus we employ the MC technique, because for the purpose at
U Angular ðhjik ; rik ; r ij Þ ¼
2 E0 r ij rik þ : hjik h0 exp 2 q1 q2
ð2Þ
ð3Þ
where rij and rik are the distances between the center atom and its neighbors. hijk is the angle between j, i and k. Here i represent the center atom and j and k are the neighbor atoms. In Si-O-H the oxygen is the center atom and silicon and hydrogen are the corresponding neighbor atoms. E0 ; q1 and q2 are the constants. The total potential energy U total of a SiO2 particle is given by the following expression:
U total ¼
X X X U Pair1 ðr ij Þ þ U Pair2 ðr ij Þ þ U Angular ðhjik ; r ik ; r ij Þ i
i
ð4Þ
i
TCD have fitted their potential constants to correctly describe the bulk properties and structure of silica over a wide temperature and pressure range [7,26]. But the potential constants for OH group and silica-OH interactions are fitted in a crystalline system based on the structure and mechanical properties under ambient temperature and pressure [7,26]. Du and Cormack [7] have tested these
92
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
fitted potential constants and proved they are reliable to simulate the hydroxylation on silica surfaces. The fitted potential constants and charges are listed in the Tables 1–4. The bulk silica has only two different ions or partially charged atoms (Si, O) and three different pair interactions (Si-Si, O-O, SiO) between them (see Table 6). This makes the bulk silica easy to model. But unlike the bulk silica, the hydroxylated silica nanoparticles have four different ions (Si, O (in silica), O (in OH groups), H), one angular (Si-O-H) and ten different pair interactions, which make them more complex to model. Because the hydroxylated silica nanoparticles have so many interactions, the potential functions for their corresponding pair and angular interaction must be clearly described and understood before modeling. Here we do this by tabulating their potential functions in a table (see Table 5). In addition, we also show the pair potential energy for the corresponding pair interactions in Fig. 1 using the potential constants and charges in Tables 1–4.
Table 2 Buckingham potential constants of the TCD model [26,7].
2.2. Modeling of hydroxylated silica nanoparticles
Table 4 Three body angular potential constants of the TCD model [26,7].
First, the bulk silica melt was prepared by heating a random structure to T = 4000 K and P = 0 GPa via NPT-ensemble based Monte Carlo simulation using periodic boundary conditions. Here the simulations use TCD potential as described in the previous Section 2.1. Coulomb interactions are calculated using the Ewald summation [28–31]. After generating the bulk silica melt, the SiO2 nanoparticles were generated by cutting out a sphere from the melt structure. Free valences on the nanoparticle surfaces were saturated via additional hydroxyl groups as shown in Fig. 3. Then the saturated SiO2 nanoparticles were quenched to T = 300 K under free boundary conditions. Since the simulation uses free boundary conditions (non-periodic), Ewald summation is no more applicable for calculating the Coulomb interactions. Instead a double loop over all atom pairs for the Coulomb interactions was implemented, as shown in the following Eq. (5).
U Coulomb ¼
e2 4p0
N1 X N X
qi qj r i¼1 j¼iþ1 ij
ð5Þ
Bulk simulation: The bulk simulation was performed with 100SiO2 (100 silicon and 200 oxygen atoms) groups at T = 300 K and P = 0 GPa. The simulation consisted of an initial equilibration of 2.5 million moves followed by additional 2.5 million trail moves (see Fig. 2(a)). The average values of all the properties of interests were calculated by averaging over the final 2.5 million trail moves. From the bulk simulations, the averages of the partial and total radial distribution functions, bond angle distributions, bond distances and coordination number were calculated. Particle simulations: The particle simulations were performed with ten different sized nanoparticles (1–10 nm) at T = 300 K and P = 0 GPa. Each simulation consisted of an initial equilibration of 10 million moves followed by additional 10 million trail moves (see Fig. 2(b)). The average values of all the properties were calculated by averaging over the final 10 million trail moves. From the particle simulations, the averages of the density profile, radial dis-
Table 1 Atomic partial charges of the TCD model [26,7]. Atom
Charges
qSi qO (in silica glass) qO (in hydroxyl group) qH
+2.4 1.2 0.856 +0.256
i-j
Aij (eV)
Bij (Å1)
Cij (eV Å6)
Si2.4-O1.2 O1.2-O1.2 Si2.4-O0.856 O0.856-O0.856 O1.2-O0.856 H0.256-O0.856 H0.256-O1.2
13702.905 1844.7458 12443.824 1844.7458 1844.7458 100.0 100.0
5.1595 2.9099 5.1595 2.9099 2.9099 4.0 4.0
54.681 192.58 54.681 192.58 192.58 0.0 0.0
Table 3 Coulomb-subtracted Morse potential constants of the TCD model [26,7]. i-j
D (eV)
b (Å1)
r0 (Å)
H0.256-O0.856 H0.256-O1.2
7.0525 7.0525
5.26 5.26
0.9485 0.9485
i-j-k
E0 (eV)
h0 (deg)
q1 (Å)
q2 (Å)
H0.256-O0.856-Si2.4
12.0
118.5
2.0
3.2
Table 5 The potential functions corresponding to each pair and angular interactions in the SiO2 nanoparticles. The angular interaction is given in the following table. Pair (i-j)
Potential function
Si2.4-Si2.4 Si2.4-O1.2 Si2.4-O0.856 Si2.4-H0.256 O1.2-O1.2 O1.2-O0.856 O1.2-H0.256 O0.856-O0.856 O0.856-H0.256 H0.256-H0.256
UCoulomb UCoulomb + UBuckingham UCoulomb + UBuckingham UCoulomb UCoulomb + UBuckingham UCoulomb + UBuckingham UCoulomb + UBuckingham + UCoulombsubrt. UCoulomb + UBuckingham UCoulomb + UBuckingham + UCoulombsubrt. UCoulomb
Morse
Morse
Angular (j-i-k)
Potential function
H0.256-O0.856-Si2.4
UAngular
Table 6 The potential functions corresponding to each pair interactions in bulk silica model. Pair (i-j)
Potential function
Si2.4-Si2.4 Si2.4-O1.2 O1.2-O1.2
UCoulomb UCoulomb + UBuckingham UCoulomb + UBuckingham
tribution functions, types of OH groups, bond angle distributions, bond distances and coordination number were calculated.
3. Simulation results and discussion The type of properties investigated here from the simulations are structural and surface properties. First, we study the structural properties for different sizes of amorphous hydroxylated SiO2 nanoparticles at the temperature of 300 K and zero pressure:
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
93
Fig. 1. The pair potential energy of Si-Si, O-O, Si-O, H-H, Si-H and O-H as the function of the interatomic distance r. At large r the energy is very small. When r decreases, the attractive Si-O, O-H interactions and the repulsive Si-Si, H-H, O-O interactions increase. The repulsion region between Si-O is very small. It prevents the Si-O atoms from coming too close during the simulations.
(a) Potential energy per SiO2 group in bulk SiO2 , bulk density vs. MC-steps.
(b) Total potential energy (T. P. E) of the SiO2 nanoparticles vs. MC-steps
Fig. 2. The density and potential energy profiles vs. Monte Carlo steps. Both the profiles have the initial equilibration region (blue line) and the equilibrium data region (red line) used for study and analysis. The expectation values of the properties of interest are determined by averaging over the equilibrium data region. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3.1. Structural properties The radial density profile for different sizes of hydroxylated SiO2 spherical nanoparticles as a function of r from the center of the particles is shown in Fig. 4. The radial density profile shown in Fig. 4 is determined as follows: we find the average number of atoms in the spherical shell of radius r with the thickness of 2 Å from the center of particles. Then, we compute the radial density profile by dividing the average number of atoms by the volume of each spherical shell. Mathematically the radial atom number density profile qðrÞ is expressed as:
qðrÞ ¼
Nðr; drÞ 4pr2 dr
ð6Þ
where Nðr; drÞ is the number of atoms in a shell of radius r from the center of particle and thickness dr. 4pr 2 dr is the shell volume. All the radial density profiles show similar atom number density fluctuations around 0.066/Å3, i.e. 0.044/Å3 for oxygen and 0.022/Å3 for silicon, respectively. At smaller radius r the density
fluctuation qðrÞ is very random and noisy because of the less number of atoms. As the r becomes larger, the fluctuation qðrÞ becomes smoother. So we consider only the larger spherical shells for our analysis. The radial density profile of our amorphous nanoparticles is compared to the one of the crystalline systems with the same size to study the difference between them. Clearly the density fluctuation qðrÞ is very stronger for crystalline nanoparticles, indicates that the crystalline structure does not have uniform density in all the spherical shells. But the atomic disorder makes the density uniform throughout the system, that is why the weaker fluctuation, smaller peak heights, widths and depths in amorphous SiO2 nanoparticles (see the region between 20 and 30 Å in 7 nm particle). In the density profiles, clearly the density fluctuation of silicon reaches zero first, then the oxygen reaches, followed by hydrogen. Also, the density fluctuation of hydrogen is observed only in the surface regions. This means our simulations successfully simulated the OH group on the silica surfaces. We also observe a sudden drop of the radial density near the surface. The surface region starts
94
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98 Table 7 The surface thickness for different sizes of the amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa.
Fig. 3. Different possibilities of bond breaking and OH group saturation on the spherical silica surfaces. Trisilanol groups and isolated oxygens are removed from the surfaces because they both were not part of silica surfaces. The trisilanol group is removed by removing the corresponding silicon and 3 oxygen atoms [6]. The red line is the surface cut. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(a) 4.5nm
Particle diameter (nm)
Surface thickness (nm)
2 3 4 5 7 9 10 Other simulation studies
0.63 0.62 0.64 0.63 0.62 0.66 0.61 0.5–1 [7]
from the point where the overall number density suddenly falls, and reaches to zero, which is the end point of surface region. With these two points, starting and ending point of the surface region, the surface thickness of the nanoparticles can be easily estimated. The estimated surface thickness values are shown in Table 7. The surface thickness values almost remain constant, around 0.6 nm, and are unaffected by a change in particle size. Other simulation studies [7] have suggested values between 0.5–1 nm. In order to study amorphous hydroxylated SiO2 nanoparticles structures and their variation with particle size, we also show the partial and total RDFs of hydroxylated SiO2 spherical nanoparticles for 3 different sizes (2 nm, 3 nm and 7 nm) in Fig. 6. Clearly the first (see magnified plot) and second peaks in the RDFs of the bulk and nanoparticles are broader, indicates that the bulk and nanoparticle structures are disordered (i.e., amorphous). This suggests that the melted structure quenched to 300 K via MC simulations leads to an amorphous state. It is also found that the peak heights in the RDFs of the nanoparticles change with the particle size and are different from the bulk. The peak heights lower with decrease in the particle size because of an increase in the surface contribution and surface effects. When the particle size increases, the peak heights approach towards the bulk, indicates that the particle structure approaches the bulk structure as the particle size increases. This RDFs behavior with particle size is similar to the behavior observed in Vo Van Hoang study [1,2]. The bond distances calculated from the RDFs are shown in Table 8. The particle’s bond distances show excellent agreement with the bulk values, except for 2 nm particle, possibly resulting from the larger surface effects. The Si-Si, O-O and Si-O coordination numbers calculated from the bulk silica and amorphous hydroxylated SiO2 nanoparticles
(b) 7nm
Fig. 4. The radial number density profile of hydroxylated SiO2 nanoparticles, for both crystalline (top) and amorphous (bottom). Note: The crystalline hydroxylated SiO2 nanoparticles were prepared by using the crystal structure of SiO2 (beta-cristoabalite).
95
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98 Table 8 The first and second peak positions of the partial and total RDFs for amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica at T = 300 K and P = 0 GPa. Pair
Si-Si O-O Si-O H-O Si-H
Peak
Particle size (nm)
(Å)
2
3
7
First Second First Second First Second First First
3.16 5.17 2.61 5.14 1.60 4.00 0.98 2.2
3.14 5.10 2.61 5.08 1.60 4.05 0.98 2.2
3.13 5.10 2.61 5.08 1.60 4.06 0.98 2.2
Bulk
EXP. [15,14]
3.13 5.15 2.61 5.06 1.60 4.13 – –
3.13 5.18 2.65 4.95 1.62 4.15 0.95–1 –
Table 9 The coordination number (CN) of Si-Si, O-O and Si-O for amorphous hydroxylated SiO2 nanoparticles and bulk amorphous silica at T = 300 K and P = 0 GPa. Pair i-j
CN
Si2.4-Si2.4 Si2.4-O1.2 = 0.856 Si2.4-O1.2 O1.2 = 0.856-O1.2 = 0.856 O1.2-O1.2
(a) Si-O-Si bond angle
First First First First First
Particle size (nm)
Bulk
2
3
7
3.08 4.00 3.08 4.90 4.6
3.40 4.00 3.40 5.28 5.15
3.73 4.00 3.73 5.67 5.65
(b) O-Si-O bond angle
4.00 4.00 4.00 6.00 6.00
(c) Si-O-H bond angle
Fig. 5. The bond angle distributions (BAD) for bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa.
structures are shown in Table 9. The Si-O, O-O and Si-Si coordination numbers of bulk silica are found to be 4, 6 and 4. The Si-O coordination number equal to four means that each silicon atom is surrounded by four oxygen atoms. Therefore the basic building block of bulk silica is SiO4 tetrahedron. The coordination numbers of the nanoparticles, on the other hand, change with the particle size, and approach the bulk values as the particle size increases. However this change is not for Si2.4 - O1.2 =0.875. Bond angle distributions have been used to analyze the local structure in the silica. Si-O-Si and O-Si-O are the two major bond angles normally investigated and studied in silica structure. Here the O-Si-O, Si-O-Si and Si-O-H bond angles are investigated to analyze the local structure in the core and surface of the nanoparticles. Fig. 5 shows the investigated bond angle distributions for bulk silica and amorphous hydroxylated SiO2 nanoparticles structures at T = 300 K and P = 0 GPa. Although there is a clear short range ordering characterized by the peaks in the RDFs of the bulk and nanoparticles, the Si-O-Si and O-Si-O bond angle distributions are significantly very broad confirming the disordered structure in both cases (bulk and nanoparticles). Like the RDFs, the Si-O-Si and O-Si-O bond angle distributions of the nanoparticles change with particle size, and are slightly different from the bulk. Also,
the Si-O-Si and O-Si-O BADs of the nanoparticles approach the bulk distribution as the particle size increases. In the Si-O-Si BAD case, there is even a clear size effect evolution with a shift toward higher angles (>160 ) with decreasing particle size (see Table 10). This means the connectivity between SiO4 tetrahedral structures changes significantly as the particle size decreases. On the other hand, the peak maximum in the O-Si-O BAD of the bulk and nanoparticles is located at 109.0 , which agrees with the experimental values (109.47 , 109.50 ). It means that the surface effect present in the nanoparticles doesn’t affects the local structure inside the SiO4 tetrahedral structures. The Si-O-H bond angle is narrowly distributed at 118.5 with a deviation of 0 , which means all the Si-O-H bond angles are reproduced exactly at 118.5 . This value agrees with the experimental and theoretical value of 118.5 , found by McCarthy, Tamassia, Woon, and Thaddeus [11]. The comparison between simulated and experimental bond angles is shown in Table 10. 3.2. Surface properties In most cases, the surface properties of silica depends mainly on the surface hydroxyl group concentration and distribution [9]. The
96
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98 Table 10 Comparison of simulated and experimental bond angles for amorphous hydroxylated SiO2 nanoparticles and bulk silica at T = 300 K and P = 0 GPa. Angle
Si-O-Si O-Si-O Si-O-H
Fig. 6. The total and partial radial distribution functions for bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa. The peak heights in the Si-O (increase), O-O (increase), and O-H (decrease) RDFs of the nanoparticles approach the bulk peak heights as the particle size increases, while the first peak height in the Si-Si RDF of the nanoparticles moves away from the bulk. This RDFs behavior is similar to the one observed in Vo Van Hoang study [1,2].
surface OH group concentration and distribution on the silica determines whether the silica surface is hydrophobic or hydrophilic [9]. For instance, the silica surface with higher surface OH group concentration produces more hydrophilic effect. In contrast the surface with less surface OH group concentration produces
Particle size (nm) 2
3
7
161 109 118.5
160 109 118.5
156 109 118.5
Bulk
EXP. [14,11,15]
155 109 –
144 , 154 109.47–.50 118.5
Fig. 7. (a) Proportion of atoms. (b) Surface area-to-volume ratio, surface-to-core atom ratio and oxygen-to-silicon atom ratio profiles. (c) Surface and core atom concentration of the nanoparticles.
hydrophobic effect. For a better understanding of the surface properties of silica nanoparticles, some investigation is needed into OH group concentration. The OH group concentration investigated from the simulated hydroxylated SiO2 spherical nanoparticles is
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
Fig. 8. The surface OH group concentration on SiO2 spherical nanoparticles for different particle sizes.
shown in Fig. 8. Clearly the OH group concentration decreases as the particle size increases. Therefore, small-sized nanoparticles have higher surface OH group concentration than large-sized nanoparticles. The average surface OH group concentration on SiO2 spherical nanoparticles is estimated to be around 7.4/nm2, which is very close to the theoretical value (7.8/nm2), found by IIler [27]. Other experimental and simulation studies have suggested different OH group concentrations between 4.5 and 9/nm2 [6,7,9]. These difference are mainly due to the distinct geometry, size, initial conditions and potentials used in the simulations. The surface area-to-volume ratio, oxygen-to-silicon atom ratio, surface-to-core atom concentration and surface-to-core atom ratio profiles of the nanoparticles are shown in Fig. 7. Clearly the surface area-to-volume and surface-to-core atom ratio’s increase as the particle size decreases. Accordingly, small-sized nanoparticles have higher surface area-to-volume ratio, higher surface atom concentration and higher surface-to-core atom ratio than large-sized nanoparticles.
4. Summary In this paper, the structural and surface properties of the amorphous hydroxylated SiO2 nanoparticles with different sizes (1– 10 nm) have been studied at T = 300 K and P = 0 GPa via MC simulations using the TCD model potential. In addition, the comparison between the bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles properties simulated at the same temperature and pressure has been analyzed. The most important conclusions drawn from the study and analysis are the followings: 1. The study shows that the TCD potential is able to successfully reproduce the structure of bulk amorphous silica and amorphous hydroxylated SiO2 nanoparticles at T = 300 K and P = 0 GPa. The size-dependent structural and surface properties were determined and were shown to be in good agreement with other experimental and simulation studies. From the agreement, it can be suggested that the TCD model potential and the methodology described in this paper are suitable and reliable for the modeling and simulation of amorphous hydroxylated SiO2 nanoparticles. 2. In this work, the SiO2 nanoparticles were generated by cutting out a sphere from the bulk silica melt to the required particle size. Free valences on the nanoparticle surfaces were saturated via additional hydroxyl groups and then quenched to T = 300 K under free boundary conditions via MC simulations. The MC simulation results clearly show that the bulk and nanoparticle structures are disordered i.e., amorphous. This means that the melted structure
97
quenched to 300 K via MC simulations leads to an amorphous state. Therefore, MC simulation using the TCD model potential does not succeed to crystallise the bulk SiO2 and SiO2 nanoparticles at 300 K. 3. The calculations of the structural properties show that the structural properties of the nanoparticles are size-dependent and different from those of bulk silica. As expected, increasing the particle size leads to an approach of the particle’s bulk properties to the bulk properties of the (quasi) infinite system. The calculations of the density profile of the nanoparticles show that the core structure and surface thickness of the nanoparticles are unaffected by a change in particle size. 4. The simulation results suggest that there is a clear size effects on the bond angle distributions and radial distribution functions. In the Si-O-Si BAD and Si-Si RDF cases, there is even a clear size effect evolution with a shift toward higher angles and bond distances with decrease in particle size. This means the decreasing in particle size not only results in increasing the surface effects but also results in shifting the bond distances and angles towards higher values. 5. The calculations of the surface properties of the nanoparticles show that decreasing the particle size results in increasing the surface OH group concentration. Accordingly, small-sized SiO2 nanoparticles have higher surface OH group concentration than large-sized SiO2 nanoparticles. The average surface OH group concentration on the SiO2 spherical nanoparticles is estimated to be around 7.4/nm2, which is very close to the theoretical value of 7.8/nm2 found by Iler [27]. Other experimental and simulation studies have suggested different OH group concentrations between 4.5 and 9/nm2 [6,7,9]. These difference are mainly due to the distinct geometry, size, initial conditions and potentials used in the simulations.
Acknowledgements I would like to thank Prof. Dr. Reinhard Hentschke (Email: [email protected]) for giving me a chance to write this paper in his research group (http://constanze.materials.uni-wuppertal.de) and for his encouragement and support. During this work, I have learned a lot about silica modeling and simulations from his great experience and knowledge. Again I would like to thank him for proofreading this paper and for his comments and corrections. Finally I would like to thank my family and friends for their support and encouragement.
References [1] V.V. Hoang, Molecular dynamics simulation of amorphous SiO2 nanoparticles, J. Phys. Chem. 111 (2007) 44. [2] N.T. Huynh, V.V. Hoang, H. Zung, Structural properties of amorphous SiO2 nanoparticles, 2016.
98
N.K. Kaliannan / Computational Materials Science 135 (2017) 90–98
[9] L.T. Zhuravlev, The surface chemistry of amorphous silica, Zhuravlev Model Colloids Surfaces A: Physiochem. Eng. Aspects 173 (1–3) (2000) 1–38. [10] I.V. Schweigert, K.E.J. Lehtinen, M.J. Carrier, M.R. Zachariah, Structure and properties of silica nanoclusters at high temperatures, Phys. Rev. B 65 (2002) 235410. [11] M.C. McCarthy, F. Tamassia, D.E. Woon, P. Thaddeus, A laboratory and theoretical study of silicon hydroxide Si-O-H, J. Chem. Phys. 129 (18) (2008) 184301. [12] D. Douroumis, I. Onyesom, M. Maniruzzaman, J. Mitchell, Mesoporous silica nanoparticles in nanotechnology, Crit. Rev. Biotechnol. 33 (3) (2012). [13] R. Nandanwar, P. Singh, F. zia Haque, Synthesis and properties of silica nanoparticles by sol-gel method for the application of green chemistry, Mater. Sci. India 10 (1) (2013) 85–92. [14] R.L. Mozzi, B.E. Warren, Structure of vitreous silica, J. Appl. Crystallogr. 2 (4) (1969) 164–172. [15] J.R.G. DaSilva, D.G. Pinatti, C.E. Anderson, M.L. Rudee, A refinement of the structure of vitreous silica, Philos. Mag. Part B 31 (3) (1974) 713–717. [16] R. Khanna, V. Sahajwalla, Treatise on process metallurgy, Process Fundamentals, first ed., vol. 1, Elsevier, 2013. [17] D. Beljonne, J. Cornil, Multiscale Modelling of Organic and Hybrid Photovoltaics, Springer, Berlin, 2014. [18] J. Horbach, Molecular dynamics computer simulation of amorphous silica under high pressure, J. Phys.: Condensed Matter 20 (2008) 24. [19] A. Carre, Development of empirical potentials for amorphous silica, Doctoral thesis, Deparment of Physics, Johannes Gutenberg University mainz, September 2007.
[20] K. Glaser, Computational studies of silica, Doctoral thesis, Department of Chemistry, University College London, May 2011. [21] D. Herzbach, Comparison of model potentials for molecular dynamics simulation of crystalline silica, Doctoral thesis, Johannes Gutenberg, Universitat Mainz, July 2004. [22] J.D. Kubicki, A.C. Lasaga, Molecular dynamics simulations of SiO2 melt and glass: Ionic and covalent models, Am. Mineral. 73 (1988) 941–955. [23] P. Tangney, S. Scandolo, An ab initio parametrized interatomic force field for silica, J. Chem. Phys. 117 (2002) 19. [24] F. Rouquerol, J. Rouquerol, K. Sing, Adsorption by Powders and Porous Solids Principles, Methodology and Applications, Academic Press, San Diego, 1998. [25] V. Beest, G. Kramer, V. Santen, Force fields for silicas and aluminophosphates based on ab initio calculations, Phys. Rev. Lett. 64 (1990) 1955. [26] C. Massobrio, J. Du, M. Bernasconi, P.S. Salmon, Molecular Dynamics Simulations of Disordered Materials, Springer Series in Material Sciences (2015) 215. [27] Ralph K. Iler, The Chemistry of Silica: Solubility, Polymerization, Colloid and Surface Properties and Biochemistry of Silica, John Wiley and Sons, New York, 1979. [28] Leeuw, Perram, Smith, Simulation of electrostatic systems in PBC. I. Lattice sums and dielectric constants, Roy. Soc. Lond. Ser. A 373 (1980) 27–56. [29] B.A. Wells, A.L. Chaffee Ewald, Summation for molecular simulations, J. Chem. Theory Comput. 11 (2015) 3684–3695. [30] U. Essmann, L. Perera, Berkowitz, H. Lee, L. Pedersen, A smooth particle Mesh Ewald potential, J. Chem. Phys. 103 (1995) 8577–8592. [31] T. Darden, D. York, L. Pedersen, Particle mesh Ewald: a N.log(N) method for Ewald sums in large systems, J. Chem. Phys. 98 (12) (1993) 10089.