# Computational Materials Engineering: Recent Applications of VASP in the MedeA^{®} Software Environment

## Article information

## Abstract

Electronic structure calculations have become a powerful foundation for computational materials engineering. Four major factors have enabled this unprecedented evolution, namely (i) the development of density functional theory (DFT), (ii) the creation of highly efficient computer programs to solve the Kohn-Sham equations, (iii) the integration of these programs into productivity-oriented computational environments, and (iv) the phenomenal increase of computing power. In this context, we describe recent applications of the Vienna Ab-initio Simulation Package (VASP) within the MedeA^{®} computational environment, which provides interoperability with a comprehensive range of modeling and simulation tools. The focus is on technological applications including microelectronic materials, Li-ion batteries, high-performance ceramics, silicon carbide, and Zr alloys for nuclear power generation. A discussion of current trends including high-throughput calculations concludes this article.

**Keywords:**Materials engineering; Computations; ab initio; Applications

## 1. Introduction

More than ever before, our society depends on a perplexing multitude of materials to meet needs such as housing, heating/cooling, clean water, production of food, energy, infrastructure, communication, transportation, and health care, as well as allowing recreational and artistic activities. It is fair to say that until now the materials necessary for all these purposes have been developed by experimental methods. Although the fundamental physical and chemical laws that govern the properties of materials are known, the goal of designing materials with specific properties by a purely theoretical/computational approach remains elusive. The reasons for this are: (i) practical materials such as stainless steel, ceramic thermal barrier coatings of turbine blades, or carbon-fiber reinforced composites contain dozens of chemical elements; (ii) their functional properties depend on the microstructure and complex interfaces as well as on the properties of individual phases; (iii) the length- and timescales underlying engineering properties span more than fifteen orders of magnitude; and (iv) all of the above properties are a function of the processing history such as heat treatment of alloys, sintering of ceramics, or curing of composite materials. For these reasons, materials development will depend on experimentation for a long time to come, yet as we show here, computational methods will make increasingly important contributions to materials engineering going forward.

The experimental development of novel materials remains an expensive and time-consuming activity and will become even more so as the requirements for high performance, low cost, and environmental compatibility become more stringent. Hence, any technology which can focus and accelerate the improvement of existing and the development of new materials is highly valuable. This is indeed the role and the challenge for computational materials engineering, as has been expressed in the concept of “Integrated Computational Materials Engineering” (ICME).1)

Simulations on the macroscopic scale are a well-established engineering practice, especially in structural analysis, in computational fluid dynamics (CFD), and in electrical circuit design. In these highly successful macroscopic simulations, materials properties data have traditionally been taken from experiment. As the predictive power of atomistic simulations increases, computed materials property data can become the input for these macroscopic simulations. This opens the exciting and unprecedented opportunity to close the loop from system design to materials design as illustrated in Fig. 1. However, as explained above, it is unrealistic to expect that one can completely replace experiment by simulations. Rather, one needs to create the synergy between atomistic simulations and experiments in the positive spirit of yin-yang or taegeuk.

In the following sections, recent applications will be presented, which demonstrate current capabilities of predicting materials properties from *ab initio* computations as part of materials engineering. These applications include published as well as unpublished work. In all cases the Vienna Ab-initio Simulation Package (VASP)2–4) within the MedeA^{®} software environment with its various tools5) has been used.

## 2. Recent Applications

### 2.1 The TiN/HfO_{2} interface in high-k dielectrics

In complementary metal oxide semiconductor (CMOS) technology, steadily diminishing device sizes have mandated the introduction of high-k dielectrics such as hafnium dioxide, which are replacing pure silicon dioxide dielectric layers. As a consequence, to maintain a low threshold voltage for switching, the material for the gate metal has had to be changed. Titanium nitride has emerged as a suitable choice in this role. A key requirement for energy efficient switching of CMOS devices is the alignment of the Fermi level (*i.e.* the energy of the highest occupied states) of the metallic gate with the band edges of the semiconducting channel of the device, as illustrated in Fig. 2.

Empirically it was found that annealing of the as-deposited TiN in an oxygen atmosphere increased the work function, as desired. Secondary ion mass spectroscopy (SIMS) measurements showed that oxygen atoms had penetrated into the layer of TiN. It was thus concluded that the replacement of oxygen in TiN causes the increase of the work function. Using MedeA^{®}-VASP, detailed electronic structure calculations of models of the HfO_{2}/TiN interface revealed, however, that replacement of N by O inside the TiN layer did not change the work function (*cf.* Fig. 2(a) and (c)). One could have concluded that computed results and experiment are in contradiction. Actually, this is not the case. While annealing of the stack in an oxygen-containing atmosphere leads to ingress of oxygen into the TiN layer, the O atoms inside the TiN layer are not the cause for the work function increase. Rather, calculations revealed that, driven by the ingress of O atoms, the diffusion of N atoms, and the filling of O vacancies in the HfO_{2} layer, the replacement of O atoms by N atoms *exactly at the interface* between HfO_{2} and TiN caused a dramatic increase of the work function,6) thus reconciling computations and experiment. The origin of this behavior is the different chemical interaction of oxygen *vs*. nitrogen with the transition metal atoms Hf and Ti. Changes in the distribution of electronic charges between the HfO_{2} and TiN layers at the interface determine the effective work function. This detailed understanding and control of the chemistry at the interface is thus critical to fabrication processes of energy efficient transistors.

### 2.2 Strength of metal/ceramic interfaces and thermal expansion: Al and Si_{3}N_{4}

Silicon nitride is a fascinating ceramic material with a wide range of applications including engine parts and ignition systems in cars, ball bearings, for example in wind turbines, in rocket thrusters due to its resistance to thermal shocks, but also in medical orthopedic devices and in the semiconductor industry as insulator and diffusion barrier. As in many practical applications, interfaces play a critical role. An illustrative example in this context is the strength of an interface between aluminum and silicon nitride.

Figure 3 shows two models of Al/Si_{3}N_{4} interfaces, namely one with Si-terminated silicon nitride and the other with N-termination. The models were constructed using the interface building tools of MedeA^{®} combined with energy minimization performed with VASP, resulting in the structures shown in Fig. 3. As a measure of the strength of the interface, the work of separation is computed from the energy difference between the interface model and the energies of the corresponding free and relaxed surfaces.

The difference between the two models is striking. While the Si-terminated interface results in a rather weak bonding between the two materials, the presence of N atoms between the surface Si atoms and the Al atoms of the metallic layer leads to a very strong cohesion between silicon nitride and aluminum.

The ceramic and the metallic phases have significantly different coefficients of thermal expansion, as shown in Fig. 4. If the interface is formed at ambient temperature, heating of the system will create a high compressive stress within the metallic phase, which is likely to lead to misfit dislocations and partial decohesion. Modeling such a complex non-equilibrium process requires a multi-scale approach, as will be discussed in the last sections of this review.

Here, the coefficient of thermal expansion is computed on the *ab initio* level using the so-called quasi-harmonic approximation. From phonon calculations for a range of different lattice parameters one obtains the vibrational entropy. Combined with the electronic energy this in turn gives expressions for the Gibbs free energy as a function of temperature. As the temperature is increased, the minimum of the Gibbs free energy shifts to larger lattice parameters. Analysis of this temperature dependence gives the coefficient of thermal expansion as used, for example, in the case of Mg_{2}SiO_{4}.7) Developed by K. Parlinski,8) the integration and automation of this capability within MedeA^{®} greatly facilitates this task.

### 2.3 Design of low-strain cathode materials for Li-ion batteries

The volume change of active materials that accompanies charge and discharge of Li-ion batteries is a major source of degradation which limits the overall lifetime of such a battery. While a zero-strain anode material exists, namely Li_{4}Ti_{5}O_{12}, there have not been any suitable zero-or low-strain materials for cathodes. By using systematic DFT calculations, three low-strain materials have been found within the class of LiMn_{x}Cr_{y}Mg_{z}O_{4}. The most promising materials have been synthesized and characterized by X-ray diffraction and electrochemical techniques. The results are consistent with the *ab initio* predictions.9)

This work focused on oxides of the composition LiM^{1}_{x}M^{2}_{y}M^{3}_{z}O^{4} crystallizing in the spinel structure. Low-strain compounds were identified by performing systematic calculations exploiting Vegard’s law as shown in Fig. 5 for selected structures. All calculations were carried out using VASP in MedeA^{®} with the PBEsol exchange-correlation functional.10,11)

The DFT calculations also provide detailed insight into the mechanisms resulting in a near zero-strain behavior. A synergistic compensation mechanism underlies the desired property as illustrated in Fig. 6. With increasing Li concentration, the Mg-O bond lengths tend to decrease, the Mn-O bond lengths remain similar, while the Cr-O bonds tend to increase. As a result, the overall volume of the crystal structure changes little upon charging and discharging with Li ions.

This behavior is in close analogy to observation for the zero-strain mechanism for Li_{4}Ti_{5}O_{12}, where local distortions in the crystal structure likewise allow this material to keep the volume nearly unchanged upon lithium insertion.

### 2.4 The structure and properties of boron carbide

Boron carbide, B_{4}C, is one of the hardest materials known, close to cubic boron nitride and diamond. Due to these mechanical properties, it is used in applications such as armor and bulletproof vests. In nuclear power reactors boron carbide is used to control the neutron flux due to the high neutron absorption of ^{10}B and the radiation hardness and chemical stability of B_{4}C.

According to the boron-carbon phase diagram,12) a boron carbide phase exists between approximately 8 at% and 20 at% C with a melting point reaching 2450°C. The crystal structure of boron carbides consists of icosahedra connected with short linear rods of three atoms. However, the distribution of the C atoms in this structure is far from obvious. Clark and Hoard13) give a structure for B_{4}C where the icosahedra consist only of B atoms connected with C-C-C linear rods. For the boron-rich compound B_{13}C_{2} the structure given by Larson14) shows B_{12} icosahedra connected with linear rods of the composition C-B-C.

The universal cluster expansion (UNCLE) method15) based on *ab initio* calculations with VASP as implemented in MedeA^{®} offers a unique methodology for the investigation of the distribution of C atoms in boron carbides as a function of C concentration. The result for the B–C system is shown in Fig. 7.

At a concentration of 20 at % (x_{B} = 0.8) the most stable structure consists of icosahedra of the composition B_{11}C with the connecting rods of C-B-C arrangements. This is consistent with the earlier theoretical work by Mauri *et al.*16) At the higher boron concentration of the compound B_{13}C_{2}, all sites of the icosahedra are occupied by boron atoms while the rods maintain the C-B-C motifs.

The elastic coefficients, which are readily computed with MedeA^{®} using a fully automated, symmetry general approach17) reveal a stiffening of the material with increasing C concentration between B_{9}C and B_{4}C as shown in the insert in Fig. 7.

Furthermore, the computed phonon dispersions for B_{4}C reveal an isolated high-frequency mode, which originates from bond-stretching vibrations of B atoms in the C-B-C linear rods as illustrated in Fig. 8 consistent with the analysis of Lazzari *et al.* 18)

### 2.5 Optical properties of Y_{2}O_{3}

The optical properties in the spectral range of visible and ultraviolet light are determined by electronic transitions from occupied to unoccupied states. Quantitative predictions of these states require a level of theory beyond standard density functional calculations. So-called hybrid functionals such as HSE0619,20) offer a practical approach to compute excitation energies. Using this approach in VASP and the optical analysis tools in MedeA^{®}, the computed refractive index of Y_{2}O_{3} (yttria) is in good agreement with experimental data,21) as illustrated in Fig. 9.

Important optical properties can be computed with good accuracy for a range of materials including transition metal oxides such as yttria as illustrated here.

### 2.6 Hydride formation in Zircaloy

The formation of zirconium hydrides is of high concern in the operation of nuclear reactors. Corrosion of zirconium alloys used in the core of nuclear power plants produces hydrogen, and a fraction of the hydrogen diffuses into the zirconium material. When the hydrogen concentration exceeds the terminal solid solubility, the excess hydrogen starts to precipitate as hydrides. This process may lead to embrittlement with crack formation due to lower ductility of the hydrides than that of the Zr matrix.

VASP as integrated in the MedeA^{®} computational environment has been employed to study structural, thermodynamic, and elastic properties of the Zr-H system.22) The computational accuracy of this method is needed to quantify and determine the behavior of hydrogen in Zr. This becomes clear considering the small energy difference between the octahedral and tetrahedral sites for hydrogen in the Zr lattice. The electronic total energy difference between the sites is computed to be only 5.9 kJ/mol, with the tetrahedral site being energetically favorable. Vibrational effects can readily be added using MedeA^{®}-Phonon. Inclusion of vibrational effects change the energy difference between the sites to 0.5 kJ/mol at T = 0 K and to 8.6 kJ/mol at 600 K.

Using a thermodynamic model of the solution of H_{2} into the octahedral and tetrahedral sites, solubility isotherms and terminal solubility of H in Zr can be computed in very good agreement with experimental data. For example, the simulations predict that the *γ*-hydride phase forms at H–Zr ratios between 1.1 at high temperatures and 1.4 at low temperatures. The reported existence range of the *γ*–phase is for H-Zr ratios between 1.1 and 1.5.

The calculations also show that the hydrogen solubility increases under tensile strain and decreases under compressive strain. This leads to hydrogen migration and accumulation in expanded regions of the Zr lattice resulting in hydride precipitation. Examples of regions under tensile stress can be at a Zr/ZrO_{2} interface, at the front of a crack tip, or even in regions around Zr self-interstitial atoms. Furthermore, hydrogen is attracted to Zr vacancies and voids. The simulations show that up to six hydrogen atoms are strongly bound inside a single Zr vacancy. Clustering of vacancies into dislocation loops can lead to regions with very high local hydrogen concentration. The simulations show that hydrogen inside the vacancy loops can delay or in some cases even prevent collapse of the loops. Each of these situations lead to regions highly supersaturated with hydrogen and could be potential nucleation sites of zirconium hydrides.

A systematic study of the zirconium hydrides has been performed by successively filling tetrahedral sites in the zirconium lattice by hydrogen, probing a large number of configurations for H-Zr ratios between zero (pure α-Zr) up to complete filling of the sites at a ratio of 2.0 (ɛ-ZrH_{2}). Computation of the elastic properties of the hydrides is conveniently carried out using the automated approach.17) in MedeA^{®}. Some of the hydride structures display elastic instability, such as cubic δ-hydride with full hydrogen occupancy which can be stabilized by introducing vacancies on the hydrogen sites or by a tetragonal distortion into ɛ-ZrH_{2}. The elastic moduli of the most stable hydrides at each stoichiometry are shown in Fig. 10. The bulk modulus increases almost linearly with hydrogen concentration from pure α-Zr to ZrH_{2}. The shear moduli of the hydrides are similar to that of α-Zr while Young’s moduli of the hydrides typically are lower than for α-Zr. The clear exception is ZrH_{1.25} which has high elastic moduli. This is identified as a γ-hydride of *P*4_{2}/mmc symmetry.

### 2.7 H-induced formation of nanotunnels on SiC surfaces

The final example is the formation of nanotunnels on surfaces of silicon carbide.23) The discovery of this type of surface features is the result of a close interaction between precise experimentation and systematic density functional calculations as detailed below. Silicon carbide is a fascinating ceramic material with a range of practical applications. It is a wide band gap semiconductor of interest for power electronics and high-frequency applications. The material remains operational at high temperatures and it is resistant to radiation. In nuclear energy technology, silicon carbide is used as fuel cladding, for example in the TRISO fuel. Furthermore, silicon carbide is a bio-compatible ceramic of interest to medical applications. While the ingredients, silicon and carbon, are readily available, the synthesis of high-purity SiC wafers requires highly sophisticated approaches24) and the hardness of SiC makes processing difficult.

The richness of phenomena on silicon carbide surfaces has stimulated a number of detailed investigations of this material. For example, C-rich surfaces exhibit formation of carbon chains;25) hydrogen, which usually passivates surfaces, can induce metallization;26) Si atomic lines can form on Si-rich surfaces;27) and complex surface reconstructions have been characterized.28)

One of the origins of this richness of silicon carbide surfaces is the fact that in crystalline silicon carbide, both Si and C atoms are tetrahedrally coordinated as in crystalline silicon and in diamond. The second nearest neighbors are arranged either in a cubic arrangement or in a hexagonal lattice. The energy difference between these two arrangements is small and hence SiC exists in a multitude of polymorphs with cubic and hexagonal SiC forming the end members. The Si-C bond length of 1.89 Å in SiC is close to the geometric mean of the Si-Si bond length of 2.35 Å and that of the C-C bond length of 1.55 Å in diamond. In other words, the lattice of bulk SiC is compressively strained in comparison with the lattice of pure silicon while, conversely, it is in a state of tensile strain in comparison with the diamond lattice. These opposing strains are equilibrated in bulk silicon carbide, but the situation at surfaces is different, especially if these surfaces are either rich in carbon or in silicon.

On a Si-rich surface of cubic SiC the release of the surface stress gives rise to a remarkable behavior, when this surface is exposed to atomic hydrogen, namely the formation of nanotunnels. When a clean SiC(001) Si-rich 3 × 2 reconstructed surface is exposed to atomic hydrogen, the most reactive Si atoms in the top layer react strongly exothermally by forming the structure 2H shown in Fig. 11. A plausible binding site for the subsequent H atoms seems to be the bridge position between the Si atoms in the third layer. In fact, this has been assumed in earlier investigations, but it leads to inconsistencies with experimental vibrational data. In particular, the vibrational modes of H-atoms in the bridge position are far too low to explain the presence of a high-frequency Si-H stretching mode which would be expected for Si atoms bonded to C atoms in the form C-Si-H. Using systematic *ab initio* calculations with VASP in MedeA^{®}, another reaction scheme has been identified as illustrated in Fig. 11 Rather than binding in a bridge position in the third layer of Si atoms, H atoms can bind exothermally to Si atoms in the second layer forming a structure denoted 6H. This reaction is remarkable, because it breaks Si-Si bonds between the second and third layer leading to an outward movement or “puckering” of the Si atoms in the second layer.

The reaction from structure 2H to 6H can be interpreted as a H-induced relief of the stress in the Si-rich surface of SiC(001). An illustration of the nanotunnel thus created is shown in Fig. 12.

The computed vibrational frequencies (*cf.* Fig. 13), obtained from *ab initio* phonon calculations using MedeA^{®}-Phonon and VASP for the nanotunnel structure, agree very well with the experimental data obtained from high-resolution electron energy loss spectroscopy for SiC(001) 3 × 2 surface exposed to hydrogen and deuterium.23)

Furthermore, the computed frequencies are also consistent with the notion that Si-H stretch frequencies are shifted to higher values if the Si atoms are bonded to C atoms. Earlier experiments using infrared spectroscopy showed absorption at 2118 and 2140 cm^{−1} (Δν = 22 cm^{−1}).26) Computations using the nanotunnel model result in frequencies for the stretch modes Si1a-H and Si3a-H of 2087 cm^{−1} (not marked explicitly in Fig. 13) and 2120 cm^{−1} (Δν = 33 cm^{−1}) as discussed in Ref. 23. The earlier bridge-bonded model is inconsistent with these experimental data. Thus, *ab initio* calculations have been essential in the clarification of the remarkable nanotunnel structure of a silicon carbide surface.

## 3. Trends and Perspectives

During the past decades we have witnessed steady progress in computational materials engineering of growing industrial value. It is probably fair to say that the ability to compute total energies of ensembles of any types of atoms using density functional theory is a cornerstone for this remarkable development. This fundamental capability has enabled in-depth understanding of rather complex systems and the prediction of a range of materials properties as illustrated in the previous section for selected cases.

As example, for electronic applications we have presented a computational analysis of the interface between HfO_{2} and TiN in the context of enhancing the efficiency of transistors with high-k dielectrics. The strength of metal/ceramic interfaces as a function of the interface chemistry has been illustrated here for the case of Al/Si_{3}N_{4} together with the *ab initio* calculation of thermal expansion coefficients. The accurate prediction of lattice parameters of compounds such as transition metal oxides in the spinel structure has enabled the design of low-strain cathode materials for Li-ion batteries. An application of this methodology to boron carbides using the cluster expansion method has helped to clarify the energetically preferred arrangement of boron and carbon atoms. Furthermore, the computation of the phonon dispersions of boron carbide has allowed the assignment of vibrational frequencies to specific atomic motions, such as the oscillation of boron atoms in linear C-B-C rods in boron carbides. Building on the results of DFT calculations, it is possible to predict optical properties such as the refractive index as a function of energy, which has been shown here for yttria. In the case of metal alloys, this *ab initio* approach provides detailed insight into the interaction of hydrogen atoms with a metal leading to the formation of hydrides and embrittlement, as illustrated here for the case of hydrogen in zirconium. In surface science, *ab initio* calculations have proven to be an invaluable tool to unravel complex surface reconstructions and to investigate the interaction of atoms and molecules with surfaces. In fact, an important industrial application of *ab initio* solid state calculations is related to heterogeneous catalysis. Here, the application to a surface science problem has been demonstrated for the case of H atoms inducing the formation of nanotunnels on Si-rich 3C-SiC(001) 3 × 2 reconstructed surfaces. In combination with precise experimental determinations of the vibrational properties of this system, *ab initio* calculations have revealed energetically favorable structures and have allowed the characterization of these structures by aligning experimental and computed vibrational frequencies, thus leading to the discovery of novel nanotunnel features.

These examples are only a small selection of all the *ab initio* calculations which are currently being performed worldwide for a vast number of different materials. In fact, it is realistic to estimate that several million DFT calculations are currently being performed per year accumulating an unprecedented volume of data. This signifies a paradigm shift. Electronic structure calculations on solid state materials have a long history dating back into the 1930’s, when the first augmented plane wave (APW) calculations were performed in the group of John Slater at the Massachusetts Institute of Technology. Throughout the second half of the past century, electronic structure methods have become increasingly sophisticated, but for the most part remained an academic discipline. Researchers typically spent months and sometimes years to investigate one or a handful of systems. For example, in the early 1980’s the all-electron *ab initio* calculation of the equilibrium geometry of a graphite monolayer (at the time not yet called graphene) took many hours of precious computing resources on supercomputers.29) Today, such a calculation is completed in minutes on a laptop. However, the level of theory has not changed dramatically in the last 30 years. Indeed, we still rely on density functional theory. In the 1980’s the common level of theory was the local density approximation (LDA) and the Kohn-Sham equation in the above case was solved with the then newly developed full-potential linearized-augmented-plane-wave (FLAPW) method.30) Now, more than 30 years later the most common approach is the generalized gradient approximation (GGA) using the projector augmented wave (PAW) method,31) which in essence is not all that different from the original APW method: the wave functions are still a combination of plane waves and localized atomic orbitals with numerical radial functions multiplied by spherical harmonics. The resulting C-C equilibrium distance given in 1982 was 2.450 Å or about 0.4% smaller than the experimental value in graphite. Today’s *ab initio* values are essentially the same.

One can say that during the past four or five decades, density functional calculations have evolved to a mature level. In fact, a recent systematic comparison of the major DFT codes resulted in a remarkable consistency of the computed results despite quite different algorithmic implementations.32) What is novel and perhaps revolutionary at present is the ease and rate new data can be computed for increasingly complex systems. Stimulated by the Materials Genome Initiative in the USA, high-throughput calculations for hundreds of thousands of compounds have become possible and are driving research in a number of academic groups in the world. This is truly exciting and one can expect new algorithms, new computational procedures, and new forms of data analytics to be developed. Large databases of computed results are being created for existing as well as hypothetical structures. Mining this richness of information will undoubtedly reveal new insights and novel materials. While this large volume of data is valuable, it will not resolve other major challenges of computational materials engineering, namely (i) the multi-scale and multi-physics aspects, (ii) the inherently non-equilibrium character of materials, and – last but not least – (iii) the accuracy of the *ab initio* calculations. Let us consider these aspects as they provide a perspective on necessary future developments.

In the overwhelming number of cases, materials engineering has a multi-scale and multi-physics character. For example, if one is interested in designing fracture resistant high-performance ceramics, one needs models of a polycrystalline material which are the result of a sintering process. Such models need to incorporate information about grain size, porosity, composition and properties of the crystalline grains with their defects, the structure of the intergranular interfaces, chemical segregation effects, residual stresses and perhaps charges trapped in defects. The fracture process involves long-range strain fields, elastic and plastic deformations on the mesoscale, grain boundary sliding, crack propagation, bond breaking at crack tips, and diffusion processes. Bond breaking and chemical rearrangements may cause local changes in vibrational energy which may entail thermal transport phenomena. Quantitative and statistically relevant modeling and simulations of such a system are all but routine with current simulation technology and yet the above case is “simple” from an engineering point of view. A more complicated case would be stress corrosion cracking involving an aqueous phase, which adds electrochemical aspects. While a decomposition of such materials problems into coupled discrete and continuum models is tempting, the dynamic fluctuations of atomistic and continuum domains, the vast configurational space, and the many orders of time scales involved make such a direct approach challenging. Such examples point to the need for novel and innovative theoretical and computational approaches combining coarse-graining in length and time scales in moving from electronic structure calculations to a continuum description with “fine-graining” in regions such as crack tips, where atomic-scale phenomena are decisive.

The non-equilibrium character of materials represents a major and fundamental challenge for computational materials engineering. This means that the processing history needs to be included in the modeling. Because of the inherent uncertainties, a statistical approach is needed to establish probabilities and to estimate boundaries. Possibly techniques from manufacturing and quality control may have to be combined with the methods of computational materials science to capture this aspect.

Finally, there remains the question of accuracy of current *ab initio* methods. While quantum chemists have developed approaches, which - at least in principle - can be converged to any desired degree of accuracy, this is not the case for the current form of density functional calculations. No systematic and practical way is known today to converge to the exact density functional. The limitations of DFT-GGA calculations are fairly well known, but despite intense efforts by a number of leading research groups in the world, there is no systematic and practical *ab initio* many-body approach which would allow one to compute, for example, the temperature of solid-liquid phase transitions to within a few degrees even for systems such as pure silicon.

This present situation is by no means unusual in the evolution of science and technology. Rather, it should be taken as stimulus for pushing the frontiers of science. One also has to keep in mind that engineering in all its forms is not and never will be perfect in all respects, but good engineering implies reliable control of the limits. It is not necessary to predict the exact yield stress of a particular sample of a material. What is required is the knowledge of the upper bound where this sample will *not* fail.

This general engineering principle is a good guide in computational materials engineering. Rather than seeking the ultimate accuracy, one needs computational protocols and approaches, which give reliable boundaries. This means a clear understanding of the key physical and chemical mechanisms which determine the properties of a material. The true art of computational materials engineering is knowing which aspects can be neglected while keeping the key characteristics of the problem. Sophisticated computational software environments with all their tools and capabilities facilitate this task, but in final analysis the best multi-scale and multi-physics tool remains the creative human mind abetted by the finest tools developed by the combined effort of the scientific and engineering disciplines.

## Acknowledgments

The authors are grateful to all colleagues of our company Materials Design as well as our many customers, who are providing the support for our work. Most of the calculations reported in this paper have been carried out at the computing facilities of Materials Design, Inc. and special thanks go to Jonathan Bell for his professional support of our computing infrastructure.

## References

*Ab initio*Molecular Dynamics for Liquid Metals. Phys Rev B 47(1):558. 1993;

*ab initio*Total-Energy Calculations Using a Plane-Wave Basis Set. Phys Rev B 54(16):11169. 1996;

^{®}- Materials Exploration and Design Analysis Materials Design, Inc. Angel Fire, NM, USA: 2016.

_{2}Gate Stacks. Appl Phys Lett 96:103502. 2010;

*g*-Mg

_{2}SiO

_{4}from First-Principles Calculations. J Chem Phys 117(7):3340–44. 2002;

_{2}. Phys Rev Lett 78(21):4063. 1997;

*Phys. Rev. Lett*.

**78**, 1396 (1997).

*Phys. Rev. Lett*

**102**039902 (2009).

_{4}C. p. 109–13. In : AIP Conference Proceedings. Albuquergue, NM, USA; 1986.

_{4}C Boron Carbide from a First Principles Analysis of NMR Spectra. Phys Rev Lett 87(8):085506. 2001;

*ab initio*Calculations of Stress. Phys Rev B 65(10):104104. 2002;

_{4}C Boron Carbide. Phys Rev Lett 83:3230. 1999;. and “Erratum”

*Phys. Rev. Lett.*

**85**4194 (2000).

*J. Chem. Phys.*

**124**219906 (2006).

*b*-SiC(100) Surface by Controlled

*sp*→

*sp*

^{3}Transformation. Phys Rev Lett 81(26):5868. 1998;

*C*-SiC(001)

*c*(4x2) Surface Reconstruction: Experiment and Theory. Phys Rev B 75(19):195315. 2007;

_{2}Molecule. Phys Rev B 24(2):864. 1981;