Simulation and visualization of crystal shapes and interfaces

Nature often arranges atoms in the shape of perfect crystals, but sometimes she creates defects and multiple domains. The optimal crystal shape at zero kelvin can be found via the Wulff construction, which can be only be carried out analytically for those trivial cases where next nearest neighbour interactions can be neglected. For our system of interest - videlicit the HCP non-Bravais case, numerical simulation is needed. This system is of relevance because we are modeling helium crystals. We have modeled two adjacent crystallites with different orientations in contact creating twist or tilt grain boundaries, and calculated the surface energy of the interface. Experience gained from several aspects of this project have a wider applications, including the condensed matter simulation application to sample construction for multi-domain crystals, and a visualization one for representation in the presence of grain boundaries. The optimization of sample shapes into their groundstates is also related to wavefront optimization in multimirror telescopes.


Introduction
There is a shift in materials simulations from studying ideal single crystals to crystals made of several grains with distinct orientations, and boundaries between them. This change in paradigm has resulted in interest in automated approaches to sample preparation. One general approach is that of the EU FP7 SimPhoNy project, where different codes are connected by wrappers and a formalism enabling easy interoperability, [1]. Another is the preparation of specific codes such as that of P. Lazic's code (CELLMATCH) to combine two unit cells into a supercell for a graphene system, and its substrate for Density Functional Theory (DFT) calculations with VASP, [2]. Of course such sample preparation can be carried out within Computer Aided Design (CAD) codes such as nCAD by Sgenia (part of SimPhoNy), but unlike CELLMATCH and the code we will introduce below, these do not include the simulation process.
In order to understand experiments on helium crystals, we needed to prepare crystals with several grains in order to study grain boundaries in helium HCP crystals. We are interested in boundaries between grains whose orientations are tilted or twisted. Thus we need to prepare adjacent crystal grains with different types of boundaries between them. The present study calculations HCP crystal shapes, cuts the crystals and then recombines them with tilt ed or twisted boundaries. A full account is provided on the project website [3] and additional material is in a supplementary pdf [4].

Experiments
Images of the evolution of a helium crystal on heating from the laboratory of S.G. Lipson at the Technion were collated by JA, published in [5] and presented in [4]. The crystal edges become round and facets become rough as the temperature increases. The change from flat facets with sharp edges to rough and rounded is known as the roughening transition. All crystals exhibit such behaviour, but due to equilibration times it is harder to observe in metals.

Wulff construction at T=0
Crystals at T=0 have an equilibrium shape with sharp edges and flat facets. This shape can be calculated by the Wulff construction which gives the minimal shape formed by the planes that have the minimal distance from the origin at each angle, [6]. The procedure is as follows: Consider a circular surface element, ∆s, at angle θ, φ with a normaln = (sin(θ), cos(φ), sin(θ), cos(θ)). A missing bondǫ would be contained in a unit circle if, the bond energy, b i (θ, φ) = 1. We note that b i (θ, φ) = 1 ifn ·ǫ > 0 and 0 ifn ·ǫ ≤ 0. The sum of all bonds in ∆s for the same b i is: The surface tension for this shape is: where σ b is the bond energy. The Wulff construction is analytic for simple Bravais nearest neighbour cases, and can be done numerically otherwise, with details given in [5].
Research into roughening transitions -Theoretical, Experimental and computational has been ongoing at the Technion for 40 years. Several faculty and many graduate students have been involved -Avron, Landau, Lipson, Adler, Balfour, Carmi, Baum, Geminterm, Hashibon, Polturak and most recently Saltoun.

The HCP lattice
The Hexagonal Close Packed (HCP) lattice is non-Bravais and like the other common non-Bravais lattice, that of diamond, is common to many interesting elements and compounds. HCP structures include one allotrope of solid helium, and the metal magnesium. Numerical studies on HCP are complicated by the two sublattice nature, and it requires 4 Miller indices. A discussion of this structure was given by two of us in a study of melting in magnesium [7]. The rich structure of HCP solid Helium makes the additional effort worthwhile again here. Fig.  1 shows the lattice structure at left.

Technion simulations of roughening transitions and their visualization
Crystal interfaces and their roughening can be studied numerically with SOS (Solid-On-Solid) models, or as interfaces in Ising models, either by direct simulation or for simple cases by series expansions. The SOS models have also been applied to telescope phasing. For the simulations, visualization is very important especially to check the complicated boundary conditions. Early examples are presented in [4].
The surface tension of the Next Nearest Neighbour (NNN) HCP is shown in the central frame of Fig. 1 (Mathematica Image by G. Baum, published in [5]) where the indentations indicate facets. After the initial observation of the roughening transition, interest turned to the order of roughening temperatures of the different facets. The HCP crystal shape is shown in Fig. 2 at upper left, with the 3 facets, known as a, c, and s, in the NN version and at lower left in the NNN version which has additional facets, cs and a2. In the central panel of Fig. 2 we show the Wulff surface with outlines of the crystal facets inside. The additional facets were researched with a simulation of surfaces of Ising models in the MSc thesis of A. Hashibon. We concluded that: A large NNN term was needed in order to achieve experimental/simulation agreement

Cutting and recombining crystals
The code used for the calculation of the surface free energy was originally published in the appendix of G. Baum's MSc. thesis, Technion, 1994, and restored from backups to be recycled for the present project, and took about 5 mins to recompile in LINUX gfortran after 20 years. Once we could again create crystals with the Wulff shape we began to cut and recombine in order to study the boundaries. The Wulff surface and crystal crystal shape from the present calculation are shown in the rightmost panels of Figs. 1 and 2 providing validation.

Steps in the process
Our process has several steps, whose details are presented in [4]. After the execution of an initial fortran code that calculates the surface tension, the next stages of the calculations were carried out with interactive Graphical User Interfaces (GUIs) in MATLAB. These stages included using the surface tension to calculate the Wulff shape and create an HCP crystal. This crystal was then cut, into crystal1 and crystal2, and the resulting crystals rotated and then recombined with

The MATLAB and AViz codes
There are five distinct GUI-enabled MATLAB codes for these processes, each with its own interface. They are described in detail on the website [3], and summarised in [4]. We show one example, the central cutting rot and conecting crystals.m GUI routine in Fig. 3. Datafiles of crystal boundary coordinates were found and then fed into surface energy calculation routines. We generated .xyz files as another interactive option and a tilted interface is shown on the right of Fig. 3 drawn with AViz [8].

Estimation of surface energy
The main application of the new code is the calculation of the boundary surface energy. Three such graphs as functions of twist and tilt direction are presented in Fig. 4, with the leftmost panel's boundary conditions (BC) being infinite, and the center's finite. These relate to whether or not the outer contour surface is ignored when calculating bond projections to approximate cases of missing bonds.
The parameters 0.5 v. 085 relate to a threshold for "close enough" to decide if a bond is present in the infinite boundary case. In the tilted case only 0.5 was used. In both twisted cases surface energy close to zero when rotation is zero. In the infinite case surface energy close to constant because no. of bonds do not change a lot while rotating. In the tilted case a much stronger angular dependence is seen.

Interface visualization with AViz
A very "transparent" way to categorise the structure of the interface is to use AViz for several layers on both sides, and to draw in bonds between atoms in adjoining crystals. The uncut is