Simulation of Polyhedral Crystal Growth Based on the Estimated Surface Energy of Crystallographic Planes

Polyhedral shapes can be found in crystalline materials ranging from macroscopic natural mineral solids to microscopic or nanoscopic particles. These shapes originate from the crystallographic properties of the constituting material, and the outer shape depends on several unique habit planes. In this study, polyhedral crystal growth was simulated considering the surface energy and crystallographic characteristics. A series of polyhedrons, including cube, truncated hexahedron, cuboctahedron, truncated octahedron, and regular octahedron, was targeted. First, the polyhedron’s static surface energy and dynamic energy variation during crystal growth were computed. Then, the crystal-growth process was simulated based on the energy minimization policy. Interestingly, when the simulation began with truncated hexahedral nucleus, the shape changed to a cuboctahedron; however, a certain type of truncated octahedron was obtained when starting with different types of truncated octahedrons. In addition, once converged cuboctahedron abruptly changed the shape to a truncated octahedron as the crystal became larger. These results were supported by the static and dynamic energy curves. Furthermore, the method was applied to different materials by assuming virtual parameters, yielding various morphologies.


Introduction
Crystalline materials have distinct crystallographic properties, and various natural minerals exhibit specific shapes [1], such as the cubic shape in NaCl and DOI: 10.4236/msa.2021.1211034 520 Materials Sciences and Applications hexagonal prism in quartz. Habitat planes, which are intrinsically based on the crystal structure, construct such polyhedral shapes, and atomistic characteristics determine the shape of macroscopic objects. It is interesting that the polyhedral shapes are also observed at the nanoscale. Nanoparticles are small particles with a diameter of tens or hundreds of nanometers and are used in various fields, including medicine, pharmacy, and functional materials [2] [3]. It is crucial to control the shape of the particles because various properties strongly depend on the outer shape. However, artificial shape-control methods, such as mechanical machining or molding of nanoparticles are difficult to use; thus, natural formation of shapes by crystal growth is the only possible procedure [4] [5]. In contrast to macro-crystals, the shapes of nanoparticles directly depend on crystallographic characteristics. These characteristics possibly change the shape during particle growth because the balance of the volume and surface area depends on the particle size. To predict the stable shape of the particles and the change in shape during crystal growth, the surface energy must be evaluated. However, experimentally measuring the surface energy is challenging; thus, computer simulation is available alternative, and especially molecular dynamics simulation is effective for capturing the atomistic and crystallographic properties [6] [7] [8] [9] [10], although the focused scale is limited in nanoscale. A phase-field model is suitable for simulating the crystal-growth process on a larger scale [11]. The model is often applied for crystal growth in solidification process, and complicated shapes such as dendrites can be successfully regenerated. This model is also applied for faceted surfaces, and polyhedral shapes are obtained [12] [13].
However, the surface energy should be expressed in a continuous function with respect to the crystal orientation, which makes it difficult to introduce the atomistic discrete properties. Polyhedral shapes can also be seen in a foam structure [14]. Interfacial energy plays an important role in determining the shape of the foam cells, where energy minimization leads to optimum morphology. The Kelvin cell is a well-known shape that is an equal-edged truncated octahedron, whereas the Weaire-Phelan cell is the alternate best shape [15] [16]. The authors applied phase-field model to the foam pattern optimization and obtained Kelvin cell as a stable structure [17] [18]. Although these results are not directly related to the morphology of a single particle, the interfacial energy is frequently used in conjunction with the surface energy and is applicable for this study. Nevertheless, these estimates mostly focus on macroscopic structures, ignoring the vast range from nano to macroscale. Therefore, we are motivated to demonstrate crystal growth simulation from a small nucleus to a large size and morphological change in the process using a simplified equation. To evaluate the surface energy of polyhedral crystal from atomistic model, we have demonstrated molecular dynamics simulation of polyhedral nanoparticles and analyzed the surface energy of typical crystallographic planes and edges [19]. In that study, we proposed a polynomial equation and demonstrated that the equation accurately reflected the surface energy of the polyhedral nanoparticle. In this study, the equation is applied to wide range of

Target Polyhedrons
This paper aims to study the morphological shifts occurring during the crystal growth. The volumetric dilatation process is simulated with an initial nucleus created with a specific polyhedral shape. The change in shape is accepted if the energy decreases because of the operation, and this process is repeated continuously. Therefore, the target polyhedrons should be expressed using a continuous parameter. In [19], a series of shapes illustrated in Figure 1(a) are considered. In present study, these shapes are studied, and the definition of the parameter follows. A cube is considered as a reference of the series (Figure 1(a)(i)). A truncated hexahedron is generated by equally cutting the eight vertices of the cube, as shown in Figure 1(a)(ii). When the length is cut equal to the half of the edge, the shape is called a cuboctahedron, as shown in Figure 1(a)(iii). If the vertices are cut more deeply, a truncated octahedron is generated, as shown in Figure 1(a)(iv). Finally, when the full length of the original cube is cut, the remaining body becomes a regular octahedron, as shown in Figure 1(a)(v). Hereafter, the cutting length from the cube's vertex is used as the shape parameter c, i.e., c = 0 represents a cube, 0.5 a cuboctahedron, and c = 1 a regular octahedron.
These shapes are suitable for the face-centered cubic (fcc) crystal because the faces appearing in this series are (100) and (111) planes only. Figure 1

Surface Energy Equation of Polyhedral Crystals
A regular polyhedron is defined as a polyhedron in which all the faces, edges, and vertices have equivalent shape, area, length, and geometrical symmetry.
Further, it is assumed that the surface energy U of such polyhedron is different from the contribution of the face, edge, and vertex, i.e., U F , U E and U V . In the case of U F , equal contributions are expected from all faces, and it is expressed as U F = FSu F , where F is the number of faces in the polyhedron, S is the area of a single face, and u F is the facial energy per unit area. Similarly, U E is expressed as U E = ELu E and U V = Vu V , where E and V are the number of edges and vertices in the polyhedron, respectively, L is the length of a single edge, and u E and u V are the edge energy per unit length and vertex energy per point, respectively. Then the surface energy of the polyhedron is expressed as follows: However, Equation (1) is invalid when we consider a crystalline material because not all faces are equivalent and have different characteristics depending on the crystallographic feature. Therefore, Equation (1) is modified as follows: Here, i, j, and k are indexes identifying a face, edge, and vertex, respectively, S i is the area of face i, L j is the length of edge j, and F i u , Following [19], the face, edge and vertex energies were calculated using molecular dynamics simulation, and obtained values for the specific crystallographic elements, i.e., face, edge, or vertex, are listed in Table 1. These energies, e F , e E , and e V were calculated by taking an average of the atoms on the target element, and hence the values are those per atom. In this study, the relative difference corresponds to determinant, and we assume these values as those per area, per length and per point of vertex.
Surface energy is defined as the energy difference between the surface and bulk region. Therefore, the bulk energy was also calculated using the same atomic model and estimated as e B = −8.224. The surface energy U is, conse-  (2). In addition, the general surface energy is expressed as the value per unit area, i.e., g = U/S A , where g is energy per area and S A is the total surface area of the solid, but we will denote U as the surface energy to show the size dependency explicitly.

Energy Calculation
The surface energy of the polyhedral crystals is calculated using Equation (2) and parameters in Table 1. The results are shown in Figure 2(a). In addition, the total energy G of the polyhedron, which is calculated as the sum of the surface energy U and bulk energy G B = e B V (Figure 2

Surface Energy and Total Energy
The cube has higher surface energy than tetrahedron, as shown in Figure 1(a) because the cube's (100) face has a higher energy than that of the tetrahedron's (111) face. As the shape parameter c becomes larger from 0 to 0.5, the energy decreases steadily because the high-energy area of the (100) decreases and the low-energy area of the (111) increases. This tendency changes at c = 0.5 for cuboctahedron. The decreasing gradient becomes steep instantaneously above c = 0.5, and the energy reaches a minimum between c = 0.65 and 0.70. Then the energy increases gradually until c reaches 1 (octahedron). Geometrically, a sphere has the minimum surface area for the same volume body. Among the polyhedrons in this study, the truncated octahedron has the minimum surface area. Nevertheless, the total edge length is the shortest for the cube, and the longest for the cuboctahedron. Because of the contradictory effects, the minimum point appears around c = 0.65 -0.70, and the minimum point slightly varies depending on the size. The size dependency is more apparent in the total energy, as shown in Figure  2(b). As the volume becomes larger, the effect of the bulk energy becomes dominant, and the energy difference between different shapes reduces. Nevertheless, it should be noted that this is a relative effect of the shape on the total energy, and that the absolute difference in energies between different shapes is still significant.

Surface Energy of Different Crystal
The surface energy depends on the crystallographic characteristics. Here, the reference values of the surface energy are virtually changed. Assuming a crystal for which the (100) face is more stable than the (111) face, the energies for (100) and (111) in Table 1 were transposed. Likewise, the energies of the 111-111 and 100-100 edges were also exchanged, and this model is referred as (100)-base model. In addition, another model that has an identical u F for all faces, identical u E for all edges, and identical u V for all vertices was considered, and this model is referred as constant-value model. The calculated results for the (100)-base model and constant-value model are shown in Figure 3(a) and Figure 3(b), respectively. In the (100)-base model, low energy is assigned to the (100) face and the cube exhibits the lowest energy for L10 model, as shown in Figure 3(a). However, when the size increases, the energy in the intermediate range reduces and the minimum point appeared around c = 0.5. This is because the total surface area is smaller than the cube, and the total surface energy decreases even when it has a high-energy face. A similar tendency is also observed when a constant value is  assigned, and the result is shown in Figure 3(b). This model demonstrates the direct effect of the surface area, and the minimum point appears around c = 0.6. In the case of L10 size, local minima are also observed at c = 0 and c = 0.45, whereas it disappears as the size increases. This indicates that the morphological change occurs during crystal growth from a small nucleus, which is discussed in the later section.

Variation in Energy during Crystal Growth
According to the well-known classical nucleation theory, the crystal-growth process consists of nucleation and growth. As the crystal grows larger, the surface area S and volume V increase. The increasing surface area increases the energy, whereas increasing volume decreases the bulk energy. Because the relative effect of the surface is more evident than volume when the size is small, the solid shrinks to vanish to minimize energy. However, if the solid exceeds the critical size, it expands because the total energy decreases as it grows. When the crystal is assumed to have a spherical shape, the total energy of the crystal is expressed as follows: where g is the surface energy per area, Δg is the bulk energy per volume, and R is the radius of the crystal. Then, the critical radius is obtained as This theory is used to investigate the polyhedral crystal targeted in this study.
At first, the effect of the crystallographic character is neglected, and constant surface energy is assumed for all faces. The result is shown in Figure 4(a), where the surface and bulk energies are assumed as γ= 2 and Δg = 3, respectively.
Overall, the initial increase below the peak at the critical volume and the continuous decrease afterward are commonly observed for all shapes. The energy for a certain volume follows the sequence: sphere < cuboctahedron < truncated octahedron (denoted as Kelvin) < octahedron < cube.  Then, the crystallographic character was considered, and the values in Table 1 were applied. The calculated energies are shown in Figure 4(b). Here, the sphere is out of target and excluded from the plot. Then the order from the lowest energy is replaced from Figure 4(a), and the octahedron becomes the smallest.
Interestingly, the curves for other shapes are crossed, and the truncated octahedron (Kelvin) is higher than the others in the small range, whereas it gets lower in the large range.
Assuming the different crystallographic properties, the reference energy values of the (100)-base and constant-value models were applied as in the previous section. Figure 4(c) and Figure 4(d) show the energy curves for these models, respectively. The cube has the lowest energy in the (100)-base model, followed by the cuboctahedron and truncated octahedron. The octahedron curve crosses the cuboctahedron and truncated octahedron as the volume increases. The crossover in the energy curves during crystal growth is not observed in the constant-value model, though the difference in shape decreases. Therefore, the crossing in the energy curve is considered the typical result brought by the crystallographic characteristics.

Simulation Procedure
The crystal-growth process was simulated using the surface energy calculated in the preceding section. A nucleus is set initially with the size represented by a referential length of a cube L 0 and the shape expressed by shape parameter c 0 , and the corresponding initial energy G 0 is calculated. First, the volumetric change is tried. The reference length is tentatively increased by +ΔL with keeping the shape parameter, and the total energy G' is calculated. If the energy decreases and G' < G 0 , the volumetric change is accepted and moved to the shape-change trial. If the energy increases by the volumetric change, the reference length is reduced by −ΔL, and moves to shape change. In the shape-change trial, the shape parameter is varied randomly by +Δc or −Δc. The preliminary total energy is calculated, and if the energy decreases, the shape modification is accepted.
This procedure is repeated, and if ten successive unaccepted trials occur, the shape is deemed to be stable for that size. The calculation then continues to the next step of the volumetric change, and these processes are repeated until the crystal becomes large enough or diminishes.

Simulation Results
When the initial size of a nucleus was set below L 0 = 2.  Figure 6 shows the variation of the shape parameter. When the calculation is initiated with 0 < c 0 < 0.5, the parameter immediately changes to c = 0.5, and the shape is maintained until the 95th calculation step. When the simulation is started with 0.5 < c 0 < 1.0, the parameter changes to 0.7 c ≅ for each scenario, but the merged value gradually decreases, and finally converges to c = 0.67. As a result, at the early stages of formation, the crystal of this material has the shape of a cuboctahedron or truncated octahedron. Interestingly, at the 95th calculation step, the cuboctahedron suddenly changes shape and merges with the   The surface energy calculated in Sections 4 and 5 is used to explain these changes in shape. As shown in Figure 2, the surface energy distributes continuously, but there is a singular minimum point at c = 0.5. When a nucleus is set in the range 0 < c 0 < 0.5, the shape changes toward the energy-minimizing direction.
Because the gradient of the energy curve is negative in this range, the shape changes toward increasing c and falls into the cusp at c = 0.5. The cusp is deep when the size is small but eventually becomes shallower, and finally, the stable state pops out. The energy curve in the range 0.5 < c 0 < 1.0, there is a local minimum between c = 0.65 and 0.7. Therefore, a nucleus in this range changes its shape to the truncated octahedron represented by this minimum. Because the minimum point slightly moves as the size becomes larger, the crystal shape also changes according to the transfer. Figure 4(b) shows the transition from a cuboctahedron to a truncated octahedron in the later stages. When the volume is small, the energy curve for the cuboctahedron is lower than the truncated octahedron (Kelvin). Meanwhile, the two curves cross together, and the truncated octahedron shows lower energy. The shape shift occurred at this point, according to the energy advantage.

Different Crystal Model
Values of the reference energy of the polyhedrons were virtually varied, and simulations were performed. The (100)-base and constant-value models were investigated after Section 4.3. Figure 8(a) and Figure 8  When the initial shape parameter is 0 < c 0 < 0.5 in the (100)-base model, the shape changes to a cube. The shape changes toward cuboctahedron when the initial shape is 0.5 < c 0 < 1.0, though the converging value is not exactly 0.5 and slightly larger. Figure 9(a) shows the visual illustration of the shape. As shown in the snapshot of the 8th calculation step, the shape is almost cuboctahedron, but there is a short edge between two yellow triangular faces. The energy curve in Figure 3(a) supports these variations. Because the gradient is positive for the curve between c = 0 and 0.5, the shape shifts toward c = 0. Meanwhile, in the range between c = 0.5 and 1.0, a local minimum exists at a point slightly larger than 0.5 for a small crystal. Therefore, the stable shape is not a complete cuboctahedron but a truncated octahedron near c = 0.5.
When using the constant-value model, the shape settles into a cubic shape when the simulation is started by c 0 = 0.2, whereas the stable shape shifts to a cuboctahedron when started by c 0 = 0.4. In the energy curve for this model presented in Figure 3(d), a local maximum exists around c = 0.2. Therefore, when the initial shape is smaller than this point, the shape moves to c = 0, and if larger, it goes to c = 0.5. When the initial nucleus is a truncated octahedron between c 0 = 0.5 and 1.0, the shape parameter varies into around c = 0.6, where the local minimum exits as shown in Figure 3(d). The corresponding shape-change process is illustrated in Figure 9(b).
In these two cases, a drastic shape change did not occur in the later stage. This result is also indicated in Figure 5(c) and Figure 5

Conclusions
Polyhedral crystal growth was simulated based on the surface energy considering the crystallographic characteristics. A series of polyhedrons were targeted, including cube, truncated hexahedron, cuboctahedron, truncated octahedron, and regular octahedron. The static surface energy of the polyhedron was estimated first. Then, based on the classical nucleation and growth theory, the variation in energy as a function of the growing volume was determined. Finally, a dynamic crystal-growth process was simulated. As a result, morphological change was observed in the early stages; when starting with truncated hexahedron, the shape converged to a cuboctahedron, while it converged to a certain truncated octahedron when starting with any truncated octahedron. Furthermore, once converged octahedron drastically changed its shape to a truncated octahedron as the crystal became larger. These changes were reasonably explained by the static and dynamic surface-energy curves. Additionally, the method was applied to different materials by assuming virtual parameters, and different shapes were obtained.
Comparison with real material was not performed yet, but similar polyhedral shapes are observed in real materials [20]. Therefore, it can be concluded that the applicability of the present method to various crystalline materials was indicated.
The reference energy was taken from the atomistic model in this work, although no quantitative verification was performed. A more complete and precise evaluation of the reference energy is required for comparison with the real material system, which can be done using a suitable molecular dynamics model. In this

Conflicts of Interest
The author declares no conflicts of interest regarding the publication of this paper.