Prediction of crystal properties by numerical simulation

The prediction of crystal properties by numerical simulation has become commonplace in the last 20 years as computers have grown more powerful and theoretical techniques more sophisticated. High accuracy prediction of elastic, electronic, transport and phase properties is possible with modern methods.

Ab Initio Calculations

Ab initio or first principles calculations are any of a number of software packages making use of density functional theory to solve for the quantum mechanical state of a system. Perfect crystals are an ideal subject for such calculations because of their high periodicity. Since every simulation package will vary in the details of its algorithms and implementations, this page will focus on a methodological overview.

Basic theory

Density functional theory seeks to solve for an approximate form of the electronic density of a system. In general, atoms are split into ionic cores and valence electrons. The ionic cores (nuclei plus non-bonding electrons) are assumed to be stable and are treated as a single object. Each valence electron is treated separately. Thus, for example, a Lithium atom is treated as two bodies – Li+ and e- – while oxygen is treated as three bodies, namely O²⁺ and 2e⁻.

The “true” ground state of a crystal system is generally unsolvable. However, the variational theorem assures us that any guess as to the electronic state function of a system will overestimate the ground state energy. Thus, by beginning with a suitably parametrized guess and minimizing the energy with respect to each of those parameters, an extremely accurate prediction may be made. The question as to what one's initial guess should be is a topic of active research.

In the large majority of crystal systems, electronic relaxation times are orders of magnitude shorter than ionic relaxation times. Thus, an iterative scheme is adopted. First, the ions are considered fixed and the electronic state is relaxed by considering the ionic and electron-electron pair potentials. Next, the electronic states are considered fixed and the ions are allowed to move under the influence of the electronic and ion-ion pair potentials. When the decrease in energy between two iterative steps is sufficiently small, the structure of the crystal is considered solved.

Boundary conditions

A key choice that must be made is how many atoms to explicitly include in one's calculation. In Big-O notation, calculations general scale as O(N3) where N is the number of combined ions and valence electrons. For structure calculations, it is generally desirable to choose the smallest number of ions that can represent the structure. For example, NaCl is a bcc cubic structure. At a first guess, one might construct a cell of two interlocked cubes – 8 Na and 8 Cl – as one's unit cell. This will give the correct answer but is computationally wasteful. By choosing appropriate coordinates, one might simulate it with just two atoms: 1 Na and 1 Cl.

Crystal structure calculations rely on periodic boundary conditions. That is, the assumption is that the cell you have chosen is in the midst of an infinite lattice of identical cells. By taking our 1 Na 1 Cl cell and copying it many times along each of the crystal axes, we will have simulated the same superstructure as our 8 Na 8 Cl cell but with much reduced computational cost.

Raw output

Only a few lists of information will be output from a calculation, in general. For the ions, the position, velocity and net force on each ion are recorded at each step. For electrons, the guess as to the electronic state function may be recorded as well. Finally, the total energy of the system is recorded. From these three types of information, we may deduce a number of properties.

Calculable properties

Unit cell parameters

Unit cell parameters (a,b,c,α,β,γ) can be computed from the final relaxed positions of the ions. In a NaCl calculation, the final position of the Na ion might be (0,0,0) in picometer Cartesian coordinates and the final position of the Cl ion might be (282,282,282). From this, we see that the lattice constant would be 584 pm. For non-orthorhombic systems, the determination of cell parameters might be more complicated, but many ab-initio numerical packages have utilities to make this calculation simpler.

Once the lattice cell parameters are known, patterns for single crystal or powder diffraction can be readily predicted via Bragg's Law.

Temperature and pressure

The temperature of the system can be estimated by use of the Equipartition Theorem, with three degrees of freedom for each ion. Since ionic velocities are generally recorded at each step in the numerical simulation, the average kinetic energy of each ion is easy to calculate. There exist schemes which attempt to control the temperature of the simulation by, e.g. enforcing each ion to have exactly the kinetic energy predicted by the Equipartition Theorem (Berendsen thermostat) or by allowing the system to exchange energy and momentum with a (more massive) fictitious enclosing system (Nose-Hoover thermostat).

The net force on each ion is generally calculated explicitly at each numerical step. From this, the stress tensor of the system can be calculated and usually is calculated by the numerical package. By varying the convergence criteria, one can either seek a lowest energy structure or a structure that produces a desired stress tensor. Thus, high pressures can be simulated as easily as ambient pressures.

Elastic properties

The Young's modulus of a mineral can be predicted by varying one cell parameter at a time and observing the evolution of the stress tensor. Because the raw output of a simulation includes energy and volume, the integrated version of the Birch-Murnaghan equation of state is often used to determine bulk modulus.

Electronic density of states

The electronic density functional is explicitly used in the calculation of the electronic ground state. Packages such as VASP have an option to calculate the electronic density of states per eV to facilitate the prediction of conduction bands and band gaps.

Thermal transport properties

The Green-Kubo relations can be used to calculate the thermal transport properties of a mineral. Since the velocities of the ions are stored at each numerical step, one can calculate the time correlation of later velocities with earlier velocities. The integral of these correlations is related to the Fourier thermal coefficient.

Diffusion

By recording the ionic positions at each time step, one can observe how far, on average, each ion has moved from its original position. The mean squared displacement of each ion type is related to the diffusion coefficient for a particle undergoing Brownian motion.