Symposium CP09—Mathematical Aspects of Materials Science—Modeling, Analysis and Computations

Electric field induced collective reorientation of nematic liquid crystal molecules is of importance for fundamental science and practical applications. This reorientation is either homogeneous over the area of electrodes, as in displays, or periodically modulated, as in electroconvection. The question is whether spatially localized three-dimensional solitary waves of molecular reorientation could be created. Here we demonstrate that the electric field can produce particle-like propagating solitary waves representing self-trapped “bullets” of oscillating molecular director. The oscillations are of the same frequency as the frequency of the driving electric field. These director bullets lack fore-aft symmetry and move with a very high speed perpendicularly to the electric field. They can move either parallel or perpendicularly to the initial alignment direction. The solitons moving perpendicularly to the director have been described in recent publication by our research team, see B.-X. Li et al, Nature Communications, volume 9, 2912 (2018). The bullets are true solitons that preserve spatially-confined shapes and survive collisions. The solitons are topologically equivalent to the uniform state and have no static analogs, thus exhibiting a particle-wave duality. Their shape, speed, and interactions depend strongly on the material parameters, which opens a door for a broad range of future studies. The work is supported by NSF DMREF grant DMS-1729509.

Lyotropic Chromonci Liquid Crystals (LCLCs) consist of disk, plank-like molecules with rigid cores and ionic gourps in the periphery, that undergo phase transitions from isotropic to hexagonal columnar phases as the concentration increaases. The chromonic liquid crystal family embraces a very broad range of materials ranging from nucleotoides, DNA and RNA to dyes and food additives (e.g. Sunset Yellow), proteins and pharmaceutical products. With the addition of condensing agents, such as pEG and spermine, the columns are found to aggregate into faceted, elongated toroidal clusters, with a typical size of 200 microns, with pressures of the order of 0.1 atmospheres. Double-stranded DNA in free solution also forms toroidal aggregates in the presence of condensing agents, but with cluster sizes about 100,000 times smaller than those formed by DSCG. An important class of chromonic aggreages are those formed by bacteriophage viruses as their dsDNA spools into an ordered structure inside the capsid. Typical viral capsids measure about 50 nm and sustein pressures of the order of 30-60 atmospheres. These large pressure values result from the competition of the high bending resistence of the DNA filament with the repulsive electrostatic force between neighboring strands.
The increased awareness of the relevant role that bacteriophage viruses may have in controlling bacterial infection processes motivates our study of such systems. In this presentation, I will discuss two models of chromonic cluster formation, one purely mechanical and a more comprehensive second model that accounts for the permeability of the capsid and changes in its polyelectrolyte environment.

We analyze the Lawrence-Doniach model for highly anisotropic superconductors with layered structure in 3D. For an applied magnetic field that is perpendicular to the layers, there are two critical values for the strength of the magnetic field at which the superconductor has phase transitions. In this talk, we will present some recent work on characterizations of the first critical field at which the magnetic field first penetrates the superconductor to create defects in the material. Such characterizations are achieved by analyzing a mean field model that is the Gamma-limit of the Lawrence-Doniach model in appropriate regimes. This is joint work with Andres Contreras at New Mexico State University.

We describe recent progress on materials and devices for the direct conversion of heat to electricity applicable to the small temperature difference regime, 10-200 C. We are pursuing an idea based on the use of first order phase transformations in epitaxial ferroelectric oxide films. It is a “direct” method in the sense that there is no separate electrical generator. We develop a mathematical theory of this method and
apply it to the design of the materials and devices and the analysis of cycles. We compare theoretical predictions and the behavior of a prototype under cyclic heating/cooling. These devices provide interesting possible ways to recover the vast amounts of energy stored on earth at small temperature difference. They move heat produced by natural and man-made sources from higher to lower temperature and therefore contribute negatively to global warming.

We work within the framework of Landau-de Gennes variational theory of nematics and show how one can obtain from a given nematic material, say A, a desired aposteriori prescribed nematic material, say B, through a suitable homogenisation procedure of suitably chosen (depending on A and B) colloidal inclusions. This is joint work with Giacomo Canevari (Basque Center for Applied Mathematics).

Formulations of problems governing the behavior of the fields in composite materials and heterogeneous media with microstructure often contain internal or hidden variables that model the influence of the processes on the microstructural level. The talk discusses a set of internal variables that are naturally related to the structure of a composite. These hidden variables determine the characteristic internal scales in the composite corresponding to the microlevel processes of different characteristic wavelengths. The parameters of composite’s microgeometry are incorporated into the internal variables through the spectral measure in the Stieltjes analytic representation of the effective properties of the composite; this spectral measure determines the equations for the evolution of the internal variables. The spectral measure contains all information about the geometry of the media on the microscale and links it through the hidden variables to the behavior of the material on the macroscale.

Electromagnetic (EM) heat exchangers are devices which absorb EM radiation converting it into thermal energy to do mechanical work. Applications involving both solar and microwave energy have garnered increasing attention but have yet to utilize the potential benefits of short-wave interactions with the local geometry when the wavelength is comparable to a material’s microstructure. It has been shown for a three-layer laminate that Bragg resonance occurs in this case, and this resonance can be used to control steady state temperature and thermal runaway effects. We investigate these effects in a mathematical model approximating a porous medium for use in a heat exchanger. Classical homogenization methods average over the microscale to obtain a macroscopic description of the material, however, they are incapable of describing short-wave behavior. High-frequency homogenization methods have been developed, but they are restricted to cases incompatible with modeling heat exchangers, such as assuming a lossless medium, spatially uniform dielectric constants, and reducing Maxwell’s equations to the Helmholtz equation. Here we develop a high-frequency homogenization technique that relax these assumptions, considering a laminate geometry composed of alternating layers of lossy dielectric material and lossless fluid channels in the homogenization limit. The presence of flow in the three-layer laminate is shown to induce sharp, but controlled, increases in the spanwise average temperature. The effects of slow flow along with thermal transport throughout the effective medium are also discussed.

Metal-Organic Frameworks(MOFs) are found to be one of the best materials for gas storage because of their large surface area and high flexibility. MOFs are compound with metallic nodes and organic linkers. We have developed an automated system to design the linkers for MOF by computer. This approach formally encompasses a design space of molecules that is astronomically large and beyond the brute force exploration by computer. In order to accelerate the search of this spare, we have tested machine learning methods for forecasting several different forms of linker molecules‘ mechanical behavior from the structure. Two systems representing the structures, high order geometric moments and blurred voxelization, were taken as well-performed features for deep learning artificial neural networks to be trained for dynamic, physics and chemical properties with high accuracy. The capability of predicting the evaluation results rather than experiment or simulating calculation will accelerate the search of new materials in a vast design space. The result shows that with the system size grows, the prediction will be more accurate.

Dislocations are fundamental defects that mediate plastic deformation in metals. Studies of dislocation core structures have relied on a variety of coupling techniques to manage the far-field strain fields; one very successful approach for pure materials has been the flexible boundary condition method based on the lattice Green function. However, the treatment of dislocations in random alloys with atomic-level detail is signicantly more challenging due to the range of environments: both local relaxation of atoms as well as the long-range field must be considered, and there is lack of efficient and accurate approaches to treat both. Here, we describe a method to calculate the dislocation core of non-dilute BCC binary random alloy efficiently by combining the first-principle calculation with an optimized multicomponent Gaussian Approximation Potential (GAP) specically designed for use with flexible boundary conditions. We use density functional theory to compute the energies, forces, elastic constants and force constants of BCC-based Ti-Mo ordered and disordered structures at zero and finite temperatures. This provides the database for our optimized Ti-Mo GAP model that accurately reproduces first and second derivatives of energies for a range of structures relevant to a dislocation. The GAP model relaxes a large dislocation geometry with random chemical species, and then the force constants for the full dislocation are computed and the lattice Green function constructed. This initial step is crucial, as flexible boundary condition methods leave far-field forces unchanged through the relaxation process. Finally, boundary conditions can be used with density-functional theory to compute the core structure of a random alloy.

Electrochemical oxidation of aluminum in electrolytes that dissolve the oxide produces porous anodic oxides (PAO). For specific ranges of anodic voltages and electrolytes, highly ordered hexagonally arranged pores are formed with a well-defined scaling relationship between the spacing of nanoscale pores and the anodizing voltage. Nano PAO has been used for a wide range of applications, including dye-sensitized solar cells, biosensors, and templates for secondary nanomaterials such as carbon, polymers and semiconductor materials. Despite the diverse applications, the mechanism for self-ordering is not yet clearly understood. Recent studies have demonstrated the importance of oxide flow and mechanical stress generation during PAO formation. Here a mathematical model is developed that includes oxide flow as an essential mechanism that drives porous growth, combined with ionic migration and electrochemical reactions in the oxide. The oxide flow is driven by mechanical stresses at oxide film interfaces, volume expansion associated with electrochemical reactions and transport of species, and divergence of ionic migration flux in the oxide. A morphological stability analysis is performed at the inception of instability that leads to the formation of steady-state porous film. Shape evolution at both metal and solution interfaces is included, in contrast to our prior model in which the metal interface was fixed. Instability at the solution interface is determined by the competitive effects of destabilizing surface stress driven oxide flow and stabilizing oxygen ionic migration flux. Metal interface evolution results in a new instability mode due to pressure-driven flow driven by the surface stress at the solution interface. As a result, the present model predicts voltage-dependent pore spacing in good agreement with experimentally observed value. The model determines the sensitivity of pattern formation to critical process parameters and predicts their range of values required for the self-organization of PAO.

As the demand for lithium-ion batteries (LIBs) with high energy density, faster charging times and longer lifetimes has been increasing explosively, it has become very important to understand the huge volume change phenomenon during lithium-ion insertion and extraction in electrodes since those can affect mechanical stability in LIBs. Several finite element (FE) simulation models have been developed to investigate the diffusion of lithium-ion and the corresponding mechanical behavior in anode materials. However, all of the models share a common drawback in that they did not take the effect of hydrostatic stress into account caused by lithium-ion insertion. When a relatively high charging rate is required, it could be essential to consider the effect of hydrostatic pressure inside anodes in order to precisely describing the diffusion behavior of lithium-ion because the surface-locking problem occurs. In this study, we propose a fully coupled diffusion-stress FE model that simultaneously solves the diffusion process and stress equilibrium equations by taking into consideration the mechanical contact between the tin dioxide active layer and the copper foam, and the hydrostatic stress gradient inside the active layer as well as the complex geometric structures. The numerically calculated strain values are compared with the experimental strain data collected in operando using X-ray diffraction. The computed strain histories are in good agreement with the measured data, enabling the possibility of precisely predicting the mechanical behavior and performance of anode materials in LIBs by only using numerical FE modeling.

Phase field modeling is a powerful and versatile method used for the study of microstructure evolution, as it allows the prediction of arbitrarily complex microstructures of multi-phase systems using continuous variables such as the solute concentration, and the order parameters, which act as phase identifiers. In order to compare the phase field predictions to experimental observations, recent work has focused on developing quantitative models that accurately describe the chemical and elastic properties of the systems of interest. However, quantitatively describing the kinetic properties of a system in a phase field model is not straightforward. Although the diffusion coefficient D of the solute might be known experimentally, the mobility L of the order parameters cannot be directly measured. As a result, L is often arbitrarily fixed to describe either a diffusion-controlled process or a reaction-controlled process. In this study, a method is developed to define the value of the kinetic parameter L based on a comparison between experimental observations and simulation results. This method is applied to the study of zirconium hydride precipitation.

For an ICME (Integrated Computational Material Engineering) modeling framework used for the age-hardening aluminum alloy design and heat treatment parameters optimization, it is critical to take into account the geometric shape of precipitates, as it is tightly related to the precipitation kinetics and particles' hardening effect. The aim of this study is to present such an ICME modeling approach to describe the precipitation of disk-shaped hardening particles during aging treatment and predict the final yield strength. The classical Kampmann-Wagner Numerical (KWN) model is extended to consider the influence of disk-shaped particle morphology on growth kinetics. The extension consists of two correction factors to the growth rate equation and to the Gibbs-Thomson effect. The extended model, coupled with a metastable thermodynamic database, is applied to simulate precipitation kinetics of Al-Cu and Al-Mg-Zn alloys during aging treatment. The predicted microstructural features are in reasonable agreement with the reported experimental observations. Furthermore, a strengthening model for disk-shaped particles, which considers the size distributions of precipitates, is developed. The predicted yield strengths are compared with reported tensile test results and with predictions from other strength models. Unlike other models, the proposed strength model can reveal the strength contribution from disk-shaped precipitates without an additional tuning parameter for accounting for the impact of the mean particle spacing in the slip plane.

For classical lens design, a usual problem is to find two surfaces so that the region bounded between them, filled with a homogeneous material, refracts light in a prescribed way. For the design of metalenses, usually a surface is given and the question is to find a phase discontinuity (a function defined on the surface) such that the surface and phase discontinuity refract light in a desired way.

We provide a mathematical approach to study metasurfaces in nonflat geometries. Analytical conditions between the curvature of the surface and the set of refracted directions are introduced to guarantee the existence of phase discontinuities. Both the far field and near field cases are considered. The starting point is a vector form of Snell's law in the presence of discontinuities on interfaces.

Our results are new since metasurfaces considered in the literature are normally flat, and the incoming radiation is directly normal to the surface. Our analysis allows for the possibility of designing non-flat metasufaces for optical applications. In particular, we allow for a variable set of directions in which we want to refract radiation. There is a compatibility condition between the variable field and the surface that must be satisfied; in the case that the field is constant, this condition is trivially satisfied.

Our general analysis includes planar and spherical lenses as special cases, and we recover previous work of F. Capasso et. al. when the surface is a plane or a sphere, as one would expect.

This is joint work with C. E. Gutierrez and L. Pallucchini (Temple University).

We present careful numerical convergence studies, using parameterized curves to reach very high resolutions in two dimensions, of a level set method for multiphase curvature motion known as the Voronoi implicit interface method. Our tests demonstrate that in the unequal, additive surface tension case, the Voronoi implicit interface method does not converge to the desired limit. We then present a variant that maintains the spirit of the original algorithm, and appears to fix the non-convergence.

Abnormal grain growth (AGG) occurs when an individual grain grows faster than the mean growth rate, and so achieves a size advantage that increases with time. AGG occurs in a wide variety of polycrystalline systems, including superalloys, nanocrystalline metals, and functional ceramics, and it can be desirable or detrimental. Despite its technological importance, the mechanisms and outcomes of AGG are not fully understood. In this presentation, we review AGG phenomenology and discuss the universal characteristics of AGG, namely the growth advantage and the persistence mechanism. Two approaches to modeling AGG are developed: A statistical model predicts the frequency of abnormal growth events as a function of microstructural characteristics and grain boundary properties, but does not capture the growth trajectories or morphologies of individual grains. In contrast, a local model based on grain boundary breakaway provides insight into the effects of grain characteristics and environment on the likelihood that a given grain will grow abnormally and on its outcome morphology. We will conclude by outlining the remaining knowledge gaps that may provide fruitful areas for further research.

Many transition metal (TM) dichalcogenide (TMD) trilayer sheets may be alloyed in the upper or lower chalcogen layers or within the TM layer. Since alloying elements typically are of different size than the atoms for which they substitute, alloying leads to localized strain which may cause bending in a free standing TMD sheet. On the other hand, bending the sheet creates spatially-localized site preference for the different alloying elements. This effect can be used to impart 3D sheet shapes through composition patterning or composition patterning through sheet bending. In this presentation we focus on shape and composition programming in MoSe_2xS_2(1-x); i.e., MoSeS. We develop a compositional strategy to do this through a quantitative, first principles-informed continuum model for the coupled mechanics and alloy thermodynamics. We apply this method to optimizing the MoSeS composition profile to reproduce the topography of a Benjamin Framklin half-dollar coin. Non-uniform spatial composition distributions give rise to electric dipoles that may lie in-plane or perpendicular to the plane of the sheet (Janus structures). The shape of such compostionally-modulated Janus sheets may be dynamically actuated through application of an electric field.

I will report on recent developments in a class of algorithms, known as threshold dynamics, for computing the motion of interfaces by mean curvature. These algorithms try to generate the desired interfacial motion by alternating two very simple operations: convolution, and thresholding. I will present a simplified version of the threshold dynamics algorithm given in the work of Esedoglu and Otto (2015) for the isotropic multiphase case that does not require the use of retardation functions. I will discuss the stability and convergence of the proposed algorithm, and threshold dynamics in general, which rely heavily on the positivity of the convolution kernel and its Fourier transform. Some counterexamples in which Γ-convergence fails in the very simple isotropic, multiphase case will also be presented. Finally, I will discuss recent results on Γ-convergence of the two phase, anisotropic case for sign changing kernels. Our contribution, which will also be discussed, is enlarging the class of admissible sign changing kernels since it is possible to construct interesting anisotropies not covered by the previous work. Joint work with Selim Esedoglu.

Epitaxial growth of metal films is an established and mature technology. The management of epitaxial strain is not a technical burden and difficulty, but opens new opportunities and functionalities in heterostructure materials. I am presenting here a case where epitaxial strain is exploited for the manipulation of the magnetic anisotropy, for example for applications in magnetic data storage. The materials engineer has only a limited set of materials and process parameters available, all of which shall be entered in an analytical mathematical model. In particular are these the chemistry of substrate including its lattice constant, and the magnetic paramters of the evaporation material including its lattice constant, plus the final film thickness. From these choices the elastic constants of the film material are immediately derived by standard materials databases.
Now, an energy balance approach is used to enter the mathematical modelling. The result is a closed mathematical expression for the thickness (critical reorientation thickness) for the epitaxial film beyond which the magnetization axis switches into the film plane. The approach requires virtually no computing power.
The model is tested for three different materials systems and delivers good match.
[1] Braun A, Feldmann B, Wuttig M: Strain-induced perpendicular magnetic anisotropy in ultrathin Ni films on Cu3Au(0 0 1). Journal of Magnetism and Magnetic Materials 1997, 171:16-28.
[2] Braun A: Conversion of Thickness Data of Thin Films with Variable Lattice Parameter from Monolayers to Angstroms: An Application of the Epitaxial Bain Path. Surface Review and Letters 2003, 10:889-894.
[3] Braun A: Quantitative model for anisotropy and reorientation thickness of the magnetic moment in thin epitaxially strained metal films. Physica B: Condensed Matter 2006, 373:346-354.

The modeling of dissipative phenomena is of particular complexity. Atomistic/particle models often offer the required physical fidelity; yet, their direct simulation beyond a few nm and ns has proven a challenge even with the most advanced computing platforms. For their part, continuum models can close the length and time scale gap needed in applications; yet, the models are primarily phenomenological. In this talk we will review some existing techniques for non-equilibrium coarse-graining and present a novel physically-based computational strategy to discover dissipative PDEs from particle data. The proposed approach, which is based on an infinite-dimensional fluctuation-dissipation relation, fully computes the dissipative operator, enabling macroscopic simulation of arbitrary initial data without the need for further particle simulations.

Distributions of aggregate sizes in various polymerization processes may be described by measures on partitions of integers and sets. We explicitly compute limit shapes for several grand canonical Gibbs ensembles and show that all possible limit shapes for these ensembles fall into distinct classes determined by the asymptotics of the internal energies of aggregates. Further on, we prove the large deviations principle and compute the corresponding free energies (rate functionals) which describe the deviations from the limit shape distributions.

The perpetual conversion of either internal chemical energy, or utilization of energy from the environment into directed motion is an integral process in active matter [1]. Its respective out-of-equilibrium nature is the origin of intriguing emerging structural and dynamical properties, which are absent in passive systems. This particularly applies to soft matter systems, e.g., comprised of filaments or polymers, which renders active soft matter a promising new class of materials. The spatiotemporal dynamics of motile active matte systems is controlled by the propulsion mechanism of the active agents in combination with various direct interactions, such as steric repulsion and hydrodynamics. These direct interactions are typically anisotropic, and emerge from different sources, such as spherical and elongated particle shapes, intrinsic flexibility, pusher and puller flow fields, etc. The combination of the various aspects leads to new emergent behavior, with a possible synergistically or antagonistically effect of the various interactions [2].

Our simulation studies of prolate spheroidal microswimmers---called squirmers---in quasi-two-dimensional confinement reveal a suppression of motility-induced phase separation (MIPS) by hydrodynamic interactions in contrast to MIPS in similar non-hydrodynamic active Brownian particles (ABPs) ensembles. The fundamental difference between ABPs and squirmers is attributed to an enhanced reorientational dynamics of squirmers by hydrodynamic torques during their collisions. In contrast, for elongated squirmers, hydrodynamics interactions enhance MIPS. The transition to a phase-separated state strongly depends on the nature of the swimmer’s flow field, with an increased tendency toward MIPS for pullers, a reduced tendency for pushers. Thus, hydrodynamic interactions show opposing effects on MIPS for spherical and elongated microswimmers, and details of the propulsion mechanism of biological microswimmers may be very important to determine their collective behavior.

Hydrodynamic interactions play a particular role in systems with a large number of internal degrees of freedom like in filamentous, polymer-like structures. As simulation and analytical studies show, activity leads to swelling of flexible polymers, shrinkage and reswelling of semiflexible polymers, and an enhanced dynamics [3]. In such systems, hydrodynamics enhances shrinkage, modifies swelling significantly, and changes the intramolecular dynamics. The shrinkage, even in the presence of excluded-volume interactions, results in an enhanced packing, which might be important for polymer organization in confinement.

[1] J. Elgeti, R. G. Winkler, G. Gompper, Physics of microswimmers---single particle motion and collective behavior: a review, Rep. Prog. Phys., 78, 056601 (2015)
[2] M. Theers, E. Westphal, K. Qi, R. G. Winkler, G. Gompper, Clustering of microswimmers: Interplay of shape and hydrodynamics, Soft Matter, 14, 8590 (2018)
[3] T. Eisenstecken, G. Gompper, R. G. Winkler, Internal dynamics of semiflexible polymers with active noise, J. Chem. Phys., 146, 154903 (2017)

While first-principles methods can compute activated state energies using the nudged-elastic band method, upscaling to mesoscale mobilities requires the solution of the master equation. The general solution for mass transport coefficients involves the inversion of rate matrix: the Green function. However, while the rate matrix is typically sparse, its inverse is not; moreover, numerical algorithms for the Green function are currently only available for cases of dilute defect concentration. I will present an alternative derivation of the transport coefficients, recast as a variational problem, where the true transport coefficients are minimal against variations in state position. This expression for transport coefficients involves a thermal average over the states, making it amenable to Monte Carlo algorithms rather than kinetic Monte Carlo. The result encompasses previous computational methods for diffusion, and provides a framework for the development of new approximation methods.

The Dean-Kawasaki model is a mass-conserving stochastic partial differential equation describing mesoscopic fluctuations of a Langevin particle system in a linear mobility regime. This equation features a deterministic gradient-flow-driven drift, as well as a non-Lipschitz noise term in spatial divergence form. This model is used and simulated in physics, e.g., in the description of thin liquid films, of nucleation for colloids and macromolecules, and of nonequilibrium bacterial dynamics. However, analytical and numerical well-posedness are still open problems. This is mainly due to the noise irregularity, which can be blamed on its conservative form and its lack of Lipschitz properties. We derive a regularised version of the Dean-Kawasaki model, based on a system of weakly interacting particles following underdamped McKean-Vlasov dynamics. We provide high-probability existence and uniqueness results for our regularised model. We also discuss work in progress concerning suitable deterministic drift corrections for both the original and the regularised Dean-Kawasaki model.

Third-order elastic constants allow to quantify and characterize the nonlinear elastic properties of a material, e.g. thermal expansion, Gruneisen parameters, specific heats, sound wave mixing and attenuation, and elastic phase transitions. Measuring these physical parameters is difficult and prone to large errors, whereas current methods based on fitting energy-density vs. strain curves work well only for high symmetry crystals.
In this paper, we present and illustrate a novel method to calculate from first principles the full set of third-order elastic constants of a material of arbitrary symmetry. The method relies on a periodic first-principles scheme to to calculate the Cauchy stress tensor of a material and numerical differentiation of the second Piola-Kirchhoff stress tensor to evaluate the elastic constants. Here we discuss an implementation of this method based on the use of a plane-wave density functional theory scheme to calculate the Cauchy stress tensor. In particular,
we show that finite difference formulas lead to a cancellation of the finite basis set errors, whereas we propose simple solutions to eliminate the numerical errors arising from the use of Fourier interpolation techniques. Then, we show applications of our method to diamond, silicon, aluminum, magnesium, graphene, and a graphane conformer. In all cases, our method is shown to give results
in excellent agreement with both the experiments and previous calculations based on fitting energy-density curves, thereby demonstrating both accuracy, validity, and generality of our new computational approch to investigate nonlinear
elastic behaviors of materials.

Atomistic (e.g., molecular dynamics or Monte Carlo) simulations are important for studying biomolecular processes, gas adsorption in porous materials, and many other materials. Force fields for these classical atomistic simulations are often parameterized semi-empirically. Recently, interest has grown for multiscale modeling methods in which force fields have parameters extracted from quantum chemistry simulations. This allows the accuracy of quantum chemistry simulations to be partially transferred to the classical atomistic simulations that can probe larger distance and time scales. Here, we present an overview of recent developments in our research group for automatically extracting force-field parameters from quantum chemistry calculations: net atomic charges[1], electron cloud parameters, polarizabilities[2], dispersion coefficients[2], etc. Our method for computing polarizabilities and dispersion coefficients has many interesting mathematical features.[2] Of particular interest, the required computational time and memory scale linearly with increasing number of atoms in the unit cell. This is achieved by reformulating the polarization and dispersion equations so large matrix inversions are completely avoided.[2] Our method combines screening increments, a Schulz iteration, and Richardson extrapolation to eliminate matrix inversions. The normal method of solving the polarization equations involves matrix inversion, which when using Gaussian elimination with partial pivoting (GEPP) yields a required computational time scaling cubically with increasing number of atoms in the unit cell. Thus, our method can be efficiently applied to much larger materials than conventional approaches. We demonstrate this with an example containing more than 250000 atoms in the periodic unit cell. We introduce the associated energy expression that is minimized to compute the induced atom-in-material dipole moments in polarizable force fields, and we discuss linearly scaling computational algorithms for solving it. Finally, we present simulation results for liquid-vapor co-existence curves to test the accuracy of different first-principles force-field parameterization schemes.

[1] N. Gabaldon Limas and T. A. Manz, "Introducing DDEC6 atomic population analysis: part 4. Efficient parallel computation of net atomic charges, atomic spin moments, bond orders, and more," RSC Advances, 8 (2018) 2678-2707. http://doi.org/10.1039/c7ra11829e

[2] T. A. Manz, T. Chen, D. J. Cole, N. Gabaldon Limas, and B. Fiszbein, "Introducing DDEC6 atomic population analysis: part 5. New method to compute polarizabilities and dispersion coefficients," ChemRxiv preprint 7205693, (2018) 1-57. http://doi.org/10.26434/chemrxiv.7205693

Technological advances in synthesis and preparation of aerogels have resulted in formulations that have the mechanical integrity (while retaining flexibility) to be utilized in a broad range of applications and have overcome the initial brittleness that this class of materials was once known for. Both structural and functional aerogels show a drop in performance when subjected to certain cyclic thermal or impact loading due to the wear and formation of cracks which reduces their lifespan. Here we present a computational tool set that connects the change in thermal profile (derived experimentally from remote sensing) to structural failure and degradation. In combination with an appropriate finite element (FEM) solver, we have developed a genetic algorithm that can reconstruct the size and shape of the defective region in silica aerogels given the temperatures from a sensor grid. Available experimental data was used to benchmark and refine the computational tool set.

Atom probe tomography is a powerful material analysis technique capable of reconstructing 3D positions of the atoms constituting the investigated sample. It relies on the analysis of mass-spectrometry data [1] and requires preliminary establishing of correct correspondence between the peaks in the spectrum and the particular atomic species. This peak identification procedure is crucial for accurate interpretation of the measured data [2, 3].

Despite rapid penetration of computer technologies practically in any area of data analysis, manual peak labeling still remains the main technique for mass spectrum peaks identification. The procedure is time-consuming and vulnerable to errors [3]. Also, it is extremely difficult to take into account all the available information about the sample, the measured data and the possible correlations between detection of different ion spices simultaneously during manual analysis.

Computer-assisted peak labeling has recently been shown to provide noticeable advantages over manual procedure [2]. However, large amount of work is still to be done in this area. First, it is desirable to have reliable physical and mathematical background, which the analysis algorithm could be based at to ensure its universality and the ability of adaptation to different measurement conditions. Also, the algorithm should be able to take into account all the information, known about the investigated sample a priori, which can be extremely useful for making the analysis results more reliable.

In this contribution, we propose a solution of the tasks listed above on the base of the Bayesian approach, which proved to be a powerful and useful tool for solving a very similar peak identification problem in powder X-ray diffraction [4]. The key idea of the proposed approach is to rank ions according to their probabilities, calculated from Bayes’ formula and conditioned by the measured data. More specifically, the initial information about the investigated specimen (e.g. presence or absence of some elements, correlations between certain ions’ detection, etc.) is encoded in prior probabilities of ions and their combinations. The expected inaccuracies of the measurement and peak search are quantified by the likelihood function. Based on the measured data, posterior probabilities are calculated for the specimen models (combinations of ions, possibly describing the actual measured data). The posterior probability for an ion is defined as the total probability of all the models, containing such ion, and can be estimated even without actual construction of all those models. If several ions have already been accepted as appropriate by the user, only the models containing those selected ions are taken into consideration. For each combination of ions, the corresponding optimal peak labeling is chosen by maximization of the measured mass spectrum likelihood.

The designed approach has been applied to a number of time-of-flight mass-spectrometry datasets, measured for inorganic samples (alloys and semiconductors). Each mass spectrum was iteratively analyzed by accepting one of the top-ranked ions at each step. A simple approach, when each ion was ranked according to the measured data and the already accepted ions only, was compared with “smart” peak labeling, when the higher probabilities of observing correlated ions groups (e.g. Te⁺ together with Te⁺⁺, Te₂⁺ and TeO⁺) were taken into account. As one would expect, including the additional information into the analysis increased the quality of the obtained results and enabled reliable construction of sample models, consistent with the results of manual analysis.

1. D.J. Larson et al., “Local Electrode Atom Probe Tomography”, (Springer, NY) (2013).
2. Haley D. et al., Ultramicroscopy, 159, 338 (2015).
3. Hudson D. et al., Ultramicroscopy, 111, 480 (2011).
4. Mikhalychev A., Ulyanenkov A., J.Appl. Cryst., 50, 776 (2017).

For the study of complex technological processes the mathematical modeling, allowing to reduce time and expenses for research, is used. In this paper, the approaches to the mathematical modeling of electromagnetic stirring of silicon melt used for silicon solar cell production are investigated.
The task complexity of calculating fluid magnetohydrodynamics of silicon melt stirring under the action of a traveling magnetic field is caused by the following factors:
silicon melt is simultaneously affected by the superposition of several magnetic forces induced by two or more inductors;
the superposition of magnetic fields changes in time; conductive medium itself is moving;
a series of closed contours locally weakening the field is inside the directional solidification unit;
the inductor fields interact with each other.
This paper presents the description of the mathematical modeling approaches of electromagnetic stirring of silicon melt. The paper describes a calculation model of electromagnetic stirring realized in the multiphysical modeling software (COMSOL Multiphysics). The model includes the Maxwell equations for describing the interaction of an electromagnetic field with current-conducting silicon melt, the equation for calculating the Lorentz force, hydrodynamic equations for an incompressible liquid, mass transfer equations, including directional mass transfer and diffusion for calculating the impurity transport inside the melt.
The model features axial symmetry, inhomogeneous triangulation computing mesh. The paper substantiates the approach of using the averaging of the Lorentz force effecting the melt during several periods. This provides a significant higher rate of carrying out the calculations. The possibility of using the proposed approaches for the calculation of electromagnetic stirring under the conditions of nonsinusoidal current waveforms in inductors and within the conditions of inductors supply by voltage sources is shown.
The research results are as follows:
the requirements to the geometry and local specification of the triangulation computing mesh are determined;
the approach based on the use of the Lorentz force averaging is justified for several periods;
the possibility of calculating electromagnetic stirring under the conditions of nonsinusoidal current waveforms in inductors is provided;
the possibility of calculating electromagnetic stirring within the conditions of inductors supply by voltage source is provided.
The described approaches are illustrated by calculations of electromagnetic stirring for large volumes of silicon (G5 crucible).

We examine how the character of the electron transport that occurs within zinc oxide is shaped by the crystal temperature, the doping concentration, and the non-parabolicity associated with the lowest energy conduction band valley. The aim of this analysis is to provide the emerging community of zinc oxide based device engineers with a detailed understanding of how the electron transport varies in response to changes in these parameters. While variations in the crystal temperature and doping concentration naturally arise from the range of device applications expected for zinc oxide based devices, the non-parabolicity coefficient associated with the lowest energy valley of the conduction band may be adjusted through the application of stress. An analysis of this effect will be presented. Both steady-state and transient electron transport responses will be considered within the scope of this analysis. The device implications of these results will be probed.

The robust self-assembly of proteins into virus capsids provides inspiration for designing synthetic building blocks that could aggregate into arbitrary structures, but the timescales of capsid assembly make atomistic modeling of such processes difficult. In this work we develop a simple simulation testbed for evaluating the assembly robustness of simple protein analogues and work towards modeling the conformational state changes that occur in real capsid proteins. We demonstrate how simple building blocks are difficult to design for robust assembly, and strategies for encoding desired structures into the design of individual particles. Finally, we discuss how such modeling techniques provide an accessible on-ramp for learning molecular simulations and intuitions for thermodynamic self-assembly.

Understanding the packing of the molecules used in organic photovoltaics (OPVs) is necessary for predicting charge transport and understanding structure-property relationships of these devices. In this work, we perform molecular dynamics simulations to understand how complex macromolecules ranging from representative asphaltenes, to new electron acceptors like ITIC aggregate. We develop computational tools for initializing simulations of these molecules in isolation, as mixtures, and with polydispersity in molecular weight when appropriate. We find that pi-stacking drives the aggregation of molecules, including the formation of symmetrical stacks for planar molecules with few side-chains. Such structures are beneficial for charge mobility, a good thing for efficient organic solar cells. This work demonstrates the potential to use the oil refining waste (asphaltenes) as an ingredient for solar energy production and has potential for identifying molecular features that correlate with robust assembly for charge mobility.

TiO2 has been widely studied as a photocatalytic material due to its non-toxicity, chemical inertness, and high photocatalytic activity. Here, we explore the operational behavior of a novel TiO2 micropillar array for treating recycled wastewater in space. A ray tracing model is used to model light adsorption, and the lattice Boltzmann method is used to simulate the water flow through the pillars and the evolution of hydroxyl radicals in the vicinity of the pillars.

The fragility of glass forming alloys in their molten state is found to correlate with the alloys’ ductility as a glass. Fragility is the deviation of the temperature dependence of a liquid's viscosity from Arrhenius behavior. It was recently shown that there are qualitative differences in the way the medium range order changes with temperature in strong vs fragile glass forming liquids. The relatively long range differences in atomic order hint at differences in cooperative packing of atoms in strong vs fragile alloys, and lead us to ask the question: is atomic size mismatch sufficient to cause fragility? The answer to this question made classical molecular dynamics models to be created for alloys of size mismatched, but chemically identical atoms. Simulations of these systems were performed across a range of temperatures and the viscosity and medium range order computed as a function of temperature.

Acknowledgments: The authors gratefully acknowledge the use of the National Science Foundation’s Extreme Science and Engineering Discovery Environment (XSEDE).

Electron backscatter diffraction (EBSD) is a technique that describes the orientation of crystals in a sample. A series of deep convolutional neural networks was built to determine the orientation of a polycrystalline sample based on its electron diffraction patterns. By establishing a fixed coordinate system, a mathematical model will take jpg files of a series of rotations as training data at a rate of approximately 250 seconds on each epoch. The networks determine the orientation of a sample image at substantially higher orders than traditional physics-based forward models. Statistical techniques are applied within the model to attain a loss in accuracy of only 5-8 degrees for rotations between 0 and 360 degrees. This simple and robust python software tool can be utilized in further offline analysis to index electron diffraction patterns much more efficiently while maintaining an average level of accuracy greater than or equal to 80%.

Electron energy loss spectroscopy (EELS) is an experimental technique providing structural and chemical information about a material. Recent advances in transmission electron microscopy (TEM) technology have improved the energy resolution of the electron beam which has realized the investigation of vibrational modes within a material [1]. When coupled with the exceptional spatial resolution of approximately 1 angstrom achievable with current aberration-corrected transmission electron microscopes, high resolution EELS mapping of the atomic structure, chemical and electronic properties, and now vibrational modes is possible. In the vibrational regime of EELS, electrons can excite vibrational modes within the sample through short-range “impact” scattering and “long-range” dipole scattering. We focus on vibrational modes dominantly excited by dipole scattering. In polar dielectric materials like SiO₂, as the electron beam is incident on the material, surface phonon polariton (SPhP) signals couple strongly at thicknesses employed in typical TEM sample geometries [2]. These surface coupling effects manifest as a signal peak in the “reststrahlen” band with thickness-dependent energy losses in the EELS spectra. In a SiO₂ truncated slab geometry, as the electron beam is moved towards the SiO₂/vacuum interface, there is a suppression in the intensity of the SPhP that is accompanied by an increase in the intensity of the SiO₂/vacuum interface phonon polariton [3].
Our work aims to suppress the SPhP signal through manipulation of material nanostructure. Utilizing classical dielectric theory, we can simulate dipole scattering in TEM samples of SiO₂ [2]. Commercial software COMSOL Multiphysics can accurately simulate dipole scattering through finite element method while also considering phonon-polariton interactions. Through finite element method, we can calculate the electron energy loss probability when an electron beam transmits the material, which is corroborated with experimental spectra. We can now vary the nanostructure and compute the corresponding electron energy-loss probability for a given model geometry with ease. Current calculations show significant surface coupling suppression using a 2x2 array of nanoholes surrounding the electron beam. We believe this phenomenon is the result of increased SiO₂/vacuum interface vibrational signal coupling, which will invariably reduce surface coupling intensity via surface peak reduction from competing vibrational modes. In practice, such nanoscale holes would be easily fabricated with conventional focused ion beam (FIB) lithography.

References:
[1] O.L. Krivanek et al., Nature, 514 (2014), p. 209.
[2] K. Venkatraman et al., Microscopy 67 (suppl_1) (2018), i14.
[3] A. Konecna et al., Physical Review B (accepted).
[4] The support from National Science Foundation CHE-1508667 and the use of (S)TEM at John M. Cowley Center at Arizona State University is gratefully acknowledged.

We present recent advancements of our open-source software package Epoxpy which utilizes coarse-grained molecular models to simulate epoxy curing. Our software package includes analysis tools that automate the calculation of various material properties including the glass transition temperature and the gelation point. We utilize a combination of dynamic data as well as structural data to find the glass transition temperature for simulations at an arbitrary degree of cure. We evaluate multiple curve fitting methods in detecting the glass transition temperature including segmented regression, hyperbolic fitting, and a combination of linear and power functions. We leverage the HOOMD-Blue molecular dynamics engine with our open source dynamic bonding plugin to capture the reaction kinetics of epoxy curing, enabling explorations of property-processing relationships. These tools enable simulation of 100 nm length scales, which are necessary to observe microphase separation. Additionally, our curing time scales are commensurate with experiments, demonstrating that we capture the relevant kinetic behaviors of epoxy curing. We show that atomistic representations resolved from the coarse-grained morphologies are preferable for calculating some material properties. Multi-scale modeling in this way enables material properties that depend on atomic resolution to be computed after less expensive models are employed to access long-time dynamics. We demonstrate these techniques with a crosslinking blend of quaternary amines and bifunctional epoxides.

Anomalous basal type stacking faults (SFs) called partial SFs (PSFs) have been observed extensively within {10-12} or extension twins for a variety of hexagonal close-packed (HCP) metals. The formation of PSFs can generate a stacking sequence in which every other basal plane has been displaced by a (1/3) <10-10> vector. Based on this special structural property of HCP metals, a new crystallographic model is developed for {10-12}<-1011> twinning system.
In the new model, PSFs are formed by (1/3) <101-0> displacement vectors (or relaxations) on every other basal plane and create a faulting plane on {101-2} plane. Subsequently, a zonal-twinning mechanism with a set of relatively simple shear and shear-shuffle atomic displacements having the smallest possible magnitudes is needed to accomplish {101-2}<1-011> twinning mode. Moreover, the pattern of atomic motions predicted in the new model is in excellent agreement with the ones obtained through molecular dynamic simulations in the literature.
The atomic displacement vectors associated with the zonal-twinning mechanism are derived analytically and expressed in a mathematical form generalized for the whole twinned domain of any HCP crystal with κ=c/a ratio within a practical range of 1.5<κ=c/a<1.9. This displacement field demonstrates that the atoms must move in a cooperative and coordinated manner during twinning.
For the special case of κ=√3, a pure shuffling mechanism is predicted by the current model in which the shear type atomic displacement vectors are zero, and shear-shuffle type atomic displacement vectors become net shuffle type atomic displacement vectors without varying magnitudes with respect to the height of the twinned domain.
For HCP metals with κ>√3, the twinning direction is reversed as reported in the literature previously. However, the new model reveals that the shear-shuffle type atomic displacement vectors are redirected toward the opposite of the twinning direction after some layers close to the twin boundary.
In contrast to the classical twinning mechanism, the new twinning mechanism suggests a simple and traceable pattern of atomic motions such that the associated absolute and relative magnitudes of the atomic displacement vectors are considerably smaller. Moreover, in the new mechanism, the resulted twinned lattice has a correct HCP structure, and the invariance of plain strain condition for {101-2} plane is preserved during twinning.
The new model describes a structural property of HCP crystals in which atoms can undergo spontaneous cooperative movements within these crystals, and consequently, substantial stress drops in the stress strain curves can be observed at the instance of twin formation. This macroscopic effect is expressed by a deformation gradient which is actually indicating that the macroscopic effect of the current mechanism is a simple shear process, in agreement with all experimental observations and simulation studies in the literature on the significance of shear during twinning.

Block copolyelectrolytes are one of the leading candidate materials for use as solid-state single-ion electrolytes (for ion transport) and ion-conductive membranes (for ion exchanging) in electrochemical energy storage systems, such as Li-ion batteries and fuel cells. They self-assemble into nanostructures that enable both ionic transport and maintain structural integrity. Ion transport in these charged block copolymers strongly depends on the alignment of specific anisotropic nanostructures. In this work, a modified dissipative particle dynamics (DPD) simulation framework has been developed to systematically study ion diffusion dynamics considering various experimentally controllable factors, such as block volume fraction, Flory-Huggins parameter, block charge fraction (or ion concentration), and dielectric constant. Using a novel “diffusivity tensor” approach, we predict the phase-dependent ion diffusivity along the principal microdomain orientations, and find that the degree of anisotropy in ion diffusivity strongly correlates with that of the polymer microdomains. Ion conductivity in block copolyelectrolytes strongly depends on ion concentrations and temperatures, but weakly on other experimentally controllable parameters. Surprisingly, we discover that the inverse topology gyroid and cylindrical phases are the optimal candidates for single-ion conductors with high-flux ion conductivity, well-percolated isotropic diffusion pathways, and mechanical robustness. Finally, we find that higher ion diffusivity can be achieved by increasing the dielectric constant, which facilitates ion diffusion across block microdomain interfaces. This work significantly motivates future efforts in exploring inverse phases without removing grain boundaries in order to enhance ion transport.

As electronic propulsion continues to be a lucrative approach for space travel, the uncertainty of thruster and spacecraft component lifetimes solicit a need for a reduction in harmful effects. The recent developments of micro-architected surfaces show that these surfaces can partially negate the effects of plasma-induced sputtering and erosion. Due to the recent developments, the sparse collection of sputter yield data on micro-architected surfaces and materials remain an inhibiting obstacle.
Therefore, we develop a 3D computational model to simulate the surface morphology evolution of micro-architected tungsten (W) samples and its instantaneous sputtering (erosion) rate under xenon (Xe) plasma exposure. We quantify the dynamic evolution of micro features of surfaces with built-in architecture and the corresponding variation in the sputtering yield through a Monte Carlo based ray-tracing method and adaptive meshing techniques. From our simulation, the surface profiles show the deposition of sputtered atoms as a result of geometric re-trapping and re-growth of surface layers. We find that micro-architected W samples exposed to 300 eV Xe plasma result in sputtering yield reductions of over 50% of the corresponding value for flat samples, indicating a self-healing mechanism of micro-architected surfaces under plasma exposure. Our model provides a framework allowing rapid expansion of its capabilities, such that the problem of sputtering redeposition from complex surfaces can be easily separated into a “physics” part (in idealized scenarios) and a “geometry” part (discretized to take advantage of the physical calculations).

Surface erosion and secondary electron emission (SEE) of the channel walls of Hall-effect thrusters for space propulsion have been identified as the most critical life-limiting factors. Recently, new wall concepts based on micro-architected surfaces have been proposed to mitigate surface erosion and SEE. The idea behind these designs is to take advantage of very-high surface-to-volume ratios to reduce SEE and ion erosion by internal trapping and redeposition. This has resulted in renewed interest to study electron-electron processes in relevant thruster wall materials.
In this work, we present calculations of SEE yields in micro-porous W and BN surfaces using discrete event simulations. First, we study SEE as a function of primary electron energies between 100 eV and 1 keV and incidence angles between 0 and 90 degrees. For this, we use Monte Carlo simulations of electron multi-scattering processes, tracking the trajectories of primary electrons as they move through the materials and accounting for elastic and inelastic energy losses, as well as plasmonic/polaronic interactions. The results are then used to present the source term in ray-tracing Monte Carlo simulations of arbitrary geometry surfaces.
We find that SEE yields in the energy range expected in Hall thruster conditions, are reduced several fold when using the new surface designs compared to ideal flat surfaces. This points to the suitability of these micro-architected surface concepts to mitigate SEE-related issues in compact electric propulsion devices.

Proton-conducting ceramic fuel cells (PCFCs) have been actively studied as a next generation alternative power source, thanks to their low operation temperatures (~ 800 K) [1]. BaZr_1-xY_xO₃ (BZYO) and BaCe_1-xY_xO₃ (BCYO) perovskites have been the two mostly studied electrolyte materials. Recently, BaZr_1-x-yCe_xY_yO₃ was introduced by mixing the two materials to overcome the weak points of the two materials; high sintering temperature and low grain-boundary conductivity in BZYO, and poor chemical stability in BCYO [2].
In this study, we used lattice statics and genetic algorithm [3] to find energetically stable lattice structures among many possible structures in BaZr_1-x-yCe_xY_yO_3-y/2 (y=0.15) bulk as a function of Ce content. The lattice volume increases with the Ce content, because Ce atom is bigger than Zr atom. Despite the increase in volume, the distance among Y atoms diminishes, indicating that the Y atoms are agglomerated with the Ce content. This behavior could be further explained by the interaction between Y atoms and oxygen vacancies. We also applied this procedure to study grain-boundary structures (sigma 3 and 5) [4].

References
[1] C. Duan, J. Tong, M. Shang, S. Nikodemski, M. Sanders, S. Ricote and A. Almansoo, Science, 2015, 349, 1321-1326.
[2] N. Nasani, D. Ramasamy, S. Mikhalev, A. V. Kovalevsky, and D. P. Fagg, J. Power Sources, 2014, 278, 582-589
[3] J. Kim, D.-H. Kim, J.-S. Kim, and Y.-C. Kim, Comput. Mat. Sci., 2017, 138, 219-224.
[4] J. Cheng, J. Luo and K. Yang , Comput. Mat. Sci., 2018, 155, 92-103.

Metallic fuels have regained interest for use in advanced reactors due to their passive safety during accident conditions and excellent thermomechanical properties. Uranium-plutonium-zirconium alloys are one important class of metallic fuels used in fast reactors. However, the thermal expansion coefficients of orthorhombic alpha uranium, which exists between -40°C and 633°C, are highly anisotropic. It has been found that some grains become stressed beyond the yield point during thermal cycling due to this high thermal expansion anisotropy. The ductility of alpha uranium is also found to decrease with decreasing temperature, while the tensile strength increases. The fracture behavior is affected by the ductility changes, transitioning from ductile failure above 350°C to mixed ductile and intergranular failure between 25°C and 350°C and finally, cleavage failure below 25°C. The anisotropic thermal expansion results in shape changes in the polycrystalline alpha uranium, inducing mismatched strains between grains. The present study models the anisotropic thermal expansion of polycrystalline alpha uranium at room temperature and predicts the intergranular fracture behavior in combination with a phase field brittle fracture model. This study will be helpful in predicting the stress at which fracture initiates for different temperature changes and the fracture propagation paths through alpha uranium polycrystals.

We present a novel numerical method to determine continuum stress intensity factors (SIF) from molecular dynamics (MD) simulations. Critical values of the SIF define an intrinsic measure of the resistance of a material to propagate a crack. At atomic scales, however, fracture occurs as a series of breaking bonds, differing from the continuum description. As a consequence, a formal analog of the continuum SIFs calculated from atomistic simulations can have spatially localized, microstructural contributions that originate from varying bond configurations. The ability to characterize fracture at the atomic scale in terms of the SIFs offers both an opportunity to probe the effects of chemistry, as well as how the addition of a microstructural component affects the accuracy. A significant difficulty at atomic length scales is locating the crack tip due to blunting and atomic reconfiguration. We have developed a method to locate the effective crack tip based on separability of the radial and angular terms in the Williams expansion. The accuracy of our numerical approach is first examined for a simple model, and then applied to atomistic simulations of fracture in amorphous silica. MD simulations provide time and spatially dependent SIFs, with results that are shown to be in good agreement with experimental values for fracture toughness in silica glass.

Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.

A thin elastic sheet wrinkles to avoid compression. When the compression is uniaxial the wrinkling direction is clear, and the main challenge is to understand its local length scale. When the compression is biaxial the situation is more complex, since even the direction of the wrinkling is hard to predict. A case of particular interest arises when the compression is due to geometric incompatibility, for example when a flat sheet is wrapped around a sphere. In all these settings, the wrinkling patterns are local minima of an "energy" functional representing the elastic energy of the sheet plus the work done by relevant loads. It is fruitful to investigate how the minimum energy depends on the sheet thickness and other relevant parameters. Identification of an "energy scaling law" requires finding both a good model of the minimum-energy pattern and an ansatz-free lower bound. I will discuss some recently-understood examples where arguments of this type provide useful insight about the energetically-preferred wrinkling patterns.

Consideration of the motion of a thin system of, say, metallic grains whose exterior surface is exposed to an ambient vapor or vacuum can lead to quite complicated models, involving a large variety of physical effects. To gain intuition, we focus on a particularly simple but relatively classical model in which the exterior surface is taken to evolve by motion by surface diffusion and the grain boundaries, which separate between neighboring
grains, are taken to evolve by mean curvature motion. Even within this simplistic context, the dynamics can be quite complicated. To emphasize the wealth of possible behaviors, we consider a number of special 3D geometries and demonstrate a severe lack of uniqueness even within the class of possible steady state solutions. Recent results reflect in part joint work with V. Derkach and J. McCuan.

The dynamic fracture of brittle solids is a particularly interesting collective interaction connecting both large and small length scales. Apply enough stress or strain to a sample of brittle material and one eventually snaps bonds at the atomistic scale leading to fracture of the macroscopic specimen. We discuss a nonlocal model for calculating dynamic fracture. The force interaction is derived from a double well strain energy density function, resulting in a non- monotonic material model. The material properties change in response to evolving internal forces eliminating the need for a separate phase field to model the fracture set. The model can be viewed as a regularized fracture model. In the limit of zero nonlocal interaction, the model recovers a sharp interface evolution characterized by the classic Griffith free energy of brittle fracture with elastic deformation satisfying the linear elastic wave equation off the crack set. We conclude with a numerical analysis of the model which is joint work with Prashant Jha.

A variational model for the delamination of polymer gel thin films from rigid substrates is presented. A formal asymptotic analysis of a simplified 2D version of the underlying governing equations show that, as the film grows thinner, the absortion of the moisture of its surroundings tends to produce a homogeneous vertical stretch in the part of the film that remains bonded to the substrate and a state of free swelling in the debonded part. A transition layer, with a width comparable to the film thickness, is developed in order to connect the two swelling modes. TThis work is joint with Carme Calderer (U. Minnesota), Carlos Garavito-Garzon (U. Minnesota) and Suping Lyu (Medtronic Inc.).

Thin elastic membranes readily take on shapes wildly different from their own. Motivated by recent experiments on wrinkling patterns observed in thin floating shells, we obtain a fully rigorous asymptotic expansion of the energy (elastic and otherwise) valid in the highly-wrinkled limit. After renormalizing by the typical energy of wrinkling, we derive a coarse-grained model in which an elastically compatible pattern is assigned energy proportional to the difference between its intrinsic undeformed area and its projected area in the plane. Energetically optimal patterns therefore maximize their projected area. Surprisingly, this limiting model turns out to be explicitly solvable in a large variety of cases, including for shells whose (possibly non-constant) curvature is of a known sign. We demonstrate our methods with concrete examples and offer comparisons with experiments. What results is an ansatz-free explanation for the geometry of wrinkle patterns seen in confined thin elastic shells.

Elongated molecules with mirror symmetry have a tendency to form nematic phases, where the ground state is constant in space and molecules roughly point in the same direction. However, the introduction of a small quantity of chiral molecules can be enough to cause a helically twisted ground state, with periodicity comparible to the wavelength visible light, making them ideal for many applications related to optics. The Helical Twisting Power (HTP) is the constant of proportionality between the concentration of dopant and the frequency of the ground state in dilute systems, and quantifies the ability of the dopant to twist the hosts ground state. Typically HTP depends highly on both the dopant and host, and may increase or decrease with temperature. In this work we attempt to understand the HTP by considering a mean-field model, which is essentially a non-local Landau-de Gennes type model describing coupling between the dopant and host species, accounting for the symmetries of their long-range attractive pairwise interactions. The model contains many terms and is illusive, but by taking a simultaneous large domain and dilute dopant limit in a periodic domain, we can obtain the Oseen-Frank free energy for cholesteric liquid crystals, and can relate the macroscopic HTP with certain scalar quantities derived from the pairwise molecular interactions. In particular, we see that in essence the HTP is heavily related to the ratio of quantities relating to the strength of achiral Host-Host interaction and chiral Host-Dopant interaction, and the qualitative temperature dependence of HTP changes depending on their relative strength.

In the last decade, tremendous advances have been made in the design and fabrication of 2D materials with novel electronic structures. Celebrated examples of such materials include graphene, black phosphorus and van der Waals heterostructure. The surface conductivity in the infrared regime in some of these systems may allow for the propagation of fine-scale electromagnetic waves, called surface plasmon-polaritons (SPPs), which can defy the usual diffraction limit. In this talk, I will discuss macroscopic consequences of the optical conductivity of 2D materials via solutions of classical Maxwell's equations. In particular, I will formally discuss how geometry, i.e., the presence of material edges and curvature, can influence SPP dispersion. In addition, I will show how slowly varying waves with nearly no phase delay may possibly propagate through suitably designed plasmonic crystals.

Vertical stacking of monolayers via van der Waals (vdW) interaction opens promising routes toward engineering physical properties of two-dimensional (2D) materials and designing atomically thin devices. However, due to the lack of mechanistic understanding, challenges remain in the controlled fabrication of these structures via scalable methods such as chemical vapor deposition (CVD) onto substrates. In this talk, we develop a general multiscale model to describe the size evolution of 2D layers and predict the necessary growth conditions for vertical (initial + subsequent layers) versus in-plane lateral (monolayer) growth. An analytic thermodynamic criterion is established for subsequent layer growth that depends on the sizes of both layers, the vdW interaction energies, and the edge energy of 2D layers. We then explore the effect of nonlinearity using a new, 2^nd order accurate diffuse domain formulation of the governing equations that places no constraints on the morphologies of the islands. We find that the analytic criterion accurately captures the conditions for which the 2^nd layer grows even when the morphologies are complex (e.g., highly anistropic, have negative curvature, etc.). Our model agrees with experimental observations of monolayer and bilayer growth of graphene and other 2D materials by CVD.

We derive a continuum variational model for a two-dimensional deformable lattice of atoms interacting with a two-dimensional rigid lattice. The starting point is a discrete atomistic model for the two lattices, which are assumed to have slightly different lattice parameters and, possibly, a small relative rotation. This is a prototypical example of a three-dimensional system consisting of a graphene sheet suspended over a substrate. A discrete-to-continuum procedure is used to obtain a continuum model that recovers both qualitatively and quantitatively the behavior observed in the corresponding discrete model. The continuum model predicts that the deformable lattice develops a network of domain walls characterized by large shearing, stretching, and bending deformation that accommodates the misalignment and/or mismatch between the deformable and rigid lattices. Two integer-valued parameters, which can be identified with the components of a Burgers vector, describe the mismatch between the lattices and determine the geometry and the details of the deformation associated with the domain walls.

Crystalline materials containing easily reorientable molecular species are an increasingly important category in the search for modern functional materials. Such materials include hybrid organic–inorganic halide perovskites, a record-breaking class of photovoltaic materials. Despite current interest in these materials, orientational fluctuations of arbitrary molecules are challenging to describe mathematically and do not fall within the conventional lattice models of solid state physics. The presence of orientable molecular species thus constitute a significant hurdle to truly predictive simulation of thermodynamic and kinetic phenomena in these promising materials.

In this talk, we present a mathematical framework for constructing accurate anharmonic model Hamiltonians that express the crystal potential energy in terms of collective orientations of interacting molecular ions, which we abstract as rigid rotors. The approach is motivated by the cluster expansion framework from alloy theory and is appropriate for describing both chiral and achiral molecules. By utilizing a quaternion parameterization of molecular orientation, along with a hyperspherical harmonic basis set, we are able to construct Hamiltonian expressions that reflect all symmetries of both the host crystal and the molecular species. The resulting Hamiltonian is compact, systematically improvable, and suitable for Monte Carlo or molecular dynamics simulation. Moreover, the basis functions employed in constructing the Hamiltonian are useful as order parameters for characterizing orientational order and phase transitions. We describe results obtained via this method for the CH₃NH₃PbX₃ lead-halide hybrid perovskites, including symmetry-adapted order parameters and Hamiltonian basis functions.

The Kohn–Sham scheme of density functional theory (KS-DFT) is popularly used to solve electronic structure problems in materials for applications, like catalysis, hydrogen storage, electronic and optical applications, mechanical strength and phase stability. In this talk, I will present a machine learned DFT (ML-DFT) framework that can be used to predict the energy and forces for a given atomic structure at significantly reduced computation cost. This ML-DFT framework relies on two different machine learning models, one that allows us to predict the ground state charge density of a structure and another that allows us to predict the total energy and forces acting on each atom from the ground state charge density. We believe such a framework can be used to handle multi-component systems, systems that involve charge transfer during evolution and in principle can also be trained using a database derived from hybrid functional calculations. This ML-DFT framework requires fewer computing resources than direct KS-DFT calculations, meaning that a system can be evolved for longer durations using molecular dynamics with ML-DFT. In addition, the memory requirements in ML-DFT are limited to storing the charge densities of structures in the database. This fact suggests that system sizes (up to 4000 atoms) can be studied that are not easily accessible in KS-DFT.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Mathematical modeling of materials is emerging as an independent branch of materials science. A mathematical model of a material plays a major role in the characterization and production of real materials as it provides a prototype for virtual experimentation. It is particularly useful for modern 2D (two-dimensional) materials because of the obvious expense and difficulty in doing real experiments on actual materials. However, for such applications, it is vital that robust and computationally efficient techniques are available for solving relevant equations. Green’s function (GF) is known to be a powerful mathematical technique for materials modeling, as it provides an efficient technique for solving operator equations. Our primary thesis goes beyond this mathematical fact. We show that the GF is not merely a mathematical artifact but a physical entity that can be experimentally measured. This link makes GF an even more useful technique for modeling of 2D materials by reducing the uncertainty, which is often inevitable in a purely theoretical modeling.

The connection between the GF and the process of measurement arises from the fact that the GF corresponding to an operator is the inverse of that operator. In Physics, the process of measurement of a quantity for a system is represented by an operator equation, which gives the response of the system to an applied probe. The GF is independent of the probe. Operating on the probe, GF gives the result of that measurement. It contains, in principle, the complete information and the total contribution of all the physical processes that are included in that specific operator. If, instead of calculating the GF through various physical assumptions, mathematical steps, and computational approximations, the GF could be directly measured, it would prove to be a valuable tool for physical characterization as well as mathematical modeling of materials. The objective of this work is to advance this thesis of possible measurement of the GF for 2D semiconductors – graphene and beyond graphene.

On the experimental side, in a typical scanning probe microscopy (SPM) experiment, one measures the response of a sample to a point probe or a distribution of point probes. Interestingly, that is exactly the mathematical definition of the Green’s function – response of a material to a point probe. With the recent developments in the science and technology of SPM, it should be possible, at least in principle, to measure the GF of a material system for the SPM measurement. It would prove to be an especially useful technique for the new 2D semiconductors. In these systems, it is still a challenge to understand the relationship between their structure and their physical properties. Understanding such relationships is vital for their engineering industrial applications.

In this talk, I will describe a new and computationally powerful technique for calculating GFs for the Laplace/Poisson equation for the thermal/electrostatic response of a system of 2D composite materials. The technique is based upon a partial Fourier representation of the GF. In this representation, integration over one variable can be obtained analytically, which results in a substantial economy of computational time and effort. We have used the technique for calculating the response map of anisotropic 2D phosphorene containing a metallic inclusion. Phosphorene and its composites are materials of strong topical interest because of their potential applications in electronic and other devices. We have also used our formulation to calculate an effective value of the conductivity of the composite, which can be a useful single parameter for preliminary design of devices. Numerical results will be presented and discussed.

Carbon nanotube (CNT) yarns are promising synthetic nanomaterials for a variety of applications. To better understand the mechanism for yarn formation process, we performed mesoscopic scale distinct element method (mDEM) simulations for the stretching of CNT networks. With parameters input from full atomistic simulations, mDEM bridges across the atomistic and mesoscopic length scales. Our model predicts accurately the mechanical response of the network over a large deformation range. At small and moderate deformations, zipping relaxations along the applied strain direction dictates the microstructural evolution. At larger deformations, the occurrence of energetic elasticity promotes yarn densification, by relaxing CNT waviness and eliminating squashed pores. Next, we reveal that the mesoscale friction and film morphology are key factors for the yarn formation: While lack of friction compromises the strain-induced alignment process, phononic and polymeric friction promote CNT alignment by enabling load transfer and directed zipping relaxations, especially in networks containing long and entangled CNTs. Yarns drawn from cellular networks instead maintain high porosity, even with enhanced polymeric friction.

Cobalt-based superalloys with γ/γ’ microstructures offer great promise as candidates for next-generation high-temperature alloys for applications in, e.g. turbine blades. It is essential to understand the thermodynamic and kinetic factors that influence the microstructural evolution of these alloys in order to optimize the alloy compositions and processing steps with a goal to improve their coarsening, creep and rafting behavior. We are developing a phase field approach to study the diffusion process and to predict the equilibrium shapes of Co-Al-W γ’ precipitates. In order to obtain quantitatively predictive capabilities we use a combination of molecular dynamics simulations, 3D atom probe tomography, and thermodynamic data to obtain bulk and interfacial energies, misfit strain, elastic constants and mobilities. A particular focus of our study is to understand how different energy balances, misfit strain and kinetics affect the rafting behavior and necking of γ’ precipitates, building on a previous model [1]. Apart from the shape bifurcation behavior reported in the previous work, we observe other characteristic microstructures as a result of the competition between the bulk free energy, interfacial energy, and elastic energy. The results indicate that the form of the bulk free energy may also have an influence on the equilibrium shape of the γ’ precipitates. Bulk free energy properties such as the separation between free energy minima and the height of the barrier between them may indirectly affect the morphology of the precipitates, as the total interfacial free energy also depends on this barrier. We are examining how the form of the bulk free energy affects the coarsening and rafting behavior of superalloy models, especially when the composition changes across interfaces.

[1] A. M. Jokisaari et al., Predicting the Morphologies of γ’ Precipitates in Cobalt-Based Superalloys, Acta Materialia 141, 273-284 (2017).

Most of the studies done around the plasmonic properties of metallic nano-objects are based on a static regime. Nevertheless, nano-objects are continuously changing in size and shape driven by the vibrations of their lattice. Therefore, the plasmon properties of the nano-objects are being modulated by the lattice vibrations, and vice versa. We present the study of the interaction between the acoustic vibrations and the localized surface plasmon resonance (LSPR) which give rise to a low-frequency Raman scattering response.
The resonant Raman scattering process is described using a new approach on a single acousto-plasmonic interaction step described by a Fermi golden rule. In this effort, we introduce the concept of Raman energy density (RED) which represents the electromagnetic energy density excited by the Raman probe and modulated by the acoustic vibrations of the nano-object.
To illustrate this new concept, RED spatial distributions and the Raman spectra are calculated for a gold nanoparticle with 5 nm in radius excited at its LSPR which give a clear picture of the coupling between the plasmon and the first isotropic and anisotropic acoustic vibrations (breathing and quadrupolar modes) and these are in good agreement with the experimental reports.

The Frank-Read (FR) source is one of the main dislocation multiplication mechanisms in crystalline materials. A FR source can produce multiple dislocations which may pile up against a grain boundary or other material interfaces. In a typical polycrystalline engineering material, the diameter of a dislocation loop generated from a FR source may be hundreds of nanometers and even above, which is beyond the reach of traditional molecular dynamics (MD) simulations using a modest computational resource. Moreover, the artificial force from free surface or periodic image on FR source in an atomistic model is often significant due to its length-scale limit. In this work, using a recently developed concurrent atomistic-continuum (CAC) computational tool, which is based on a theoretical framework that unifies the atomistic and continuum description of materials, we simulate growth of a micron-sized dislocation loop through a FR source from the bottom up. At a fraction of the cost of full MD, CAC overcomes the length-scale limitation of MD and scale up the model up to several microns and even above, in which the image force effect becomes negligible. It enables us to systematically calibrate the relationship between FR source length and the critical shear stress for dislocation bowing out from the atomistic to the micrometer level. This work demonstrates the applicability of CAC in predicting the dynamics of dislocations at the mesoscopic level. It also opens up the possibilities of understanding a variety of dislocation-interface reactions which require a full atomistic resolution near an interface and a coarse-grained material description of the long-range interactions between dislocations away from the interfaces in materials.

The deformation behavior of crystalline solids in the plastic regime is primarily dictated by the motion and interaction of dislocations, which are linear defects in an otherwise perfect lattice. A number of macroscale and mesoscale plastic localization phenomena can then be directly attributed to the collective behavior of interacting dislocation populations which form characteristic patterns under different regimes. In relation to macroscopic measurements, different stages of the stress-strain curve and the transitions between them can be correlated with rearrangements in the dislocation microstructure. This justifies the need to develop a physics-based modeling approach to capture different dislocations patterns intricately coupled with the mechanics of the underlying material.

In this direction, reaction-advection-diffusion(RAD) systems are pointed out as a possible recipe to model dislocation dynamics. The approach of Discrete Dislocation Dynamics(DDD) although primarily founded on fundamental physical laws, have been unable to predict dislocation patterning during fatigue. They are prohibitively time-consuming for large number of fatigue cycles, computationally expensive for practical sample sizes, and are only relevant to very high rates of deformation. The idea to motivate the dislocation system as a RAD system is motivated by drawing analogies between interacting dislocation populations and similar phenomena occurring in chemical and biological processes, where reaction-diffusion models have played a vital role. We consider the problem of dislocation pattern formation in the fatigue of copper single crystals - dislocation veins, PSB(Persistent Slip Band) ladders, extended ladders, dislocations cells and labyrinth patterns - all appearing for different crystal orientations and strain amplitudes. More specifically, we attempt to model the formation of dislocations veins under low strain amplitudes and subsequent appearance of PSBs with the ladder structure for increased strain amplitudes. A three-dimensional generalized version of the Walgraef-Aifantis model is taken as the starting point for describing the evolution of dislocation densities. The general features of the model involve: (i) Nonlinear gradient terms appearing from a stress-velocity constitutive law for dislocations describing their mobility through space, (ii) Nonlinear reaction terms describing local dislocation production and reaction rates and (ii) Separating the dislocations into immobile, mobile edge and mobile screw dislocations resulting in a three-species reaction-advection-diffusion system. Phenomenological contributions to the reactive part of kinetic equations are proposed based on dominant operating mechanisms - immobile dislocation production, mobile source activation, dipole annihilation, mobile-immobile conversion and dipole disintegration. The system of partial differential equations subject to homogeneous Neumann boundary conditions is discretized and solved using a semi-implicit finite element scheme to obtain steady state spatial patterns. Parametric studies are conducted to explore their effect on the nature of patterns and transitions from one pattern to another, and justify their appearance based on changes in the experimental conditions and operating physical mechanisms.

Due to its outstanding electromechanical properties, lead titanate (PbTiO₃) and it's solid solutions, such as Pb(Zr, Ti)O₃ (PZT) and Pb(Mg_1/3, Nb_2/3)O₃-PbTiO₃ (PMNPT), have been widely applied in technological devices like non-volatile random access memories (FeRAMs), micro-actuator and sensors. For such materials, a slight atomic displacement strongly affects both mechanical and ferroelectric (FE) properties. Thus, a proper understanding of the atomic behavior of ferroelectrics is required. Computational simulation is a promising tool for analysis of the properties above. In this study, we have successfully reparameterized the core-shell model to model PbTiO₃. The new parameters considered the crystal structures, elastic properties and phonon frequencies which calculated by ab-initio density functional theory calculations. In addition to material mechanical properties, it can reproduce the thermodynamic properties like isometric heat capacity, adiabatic compressibility.Moreover, the new parameters were also tested by molecular dynamics simulations to examine the predictive ability at finite temperature. The phase transition temperature of 490 °C from the tetragonal to the cubic phase was obtained, which is consistent with experimental results.Furthermore, the new model was then used to study the interfacial effect on ferroelectricity at the nanoscale. The observed morphology changes of the ferroelectric domains at various temperatures, with or without an electric field in bulk, thin film and nano-rods match with the experiments. From the simulated models, a deeper understanding of the domain motion in atomic scale was achieved.

In this work, the mobility of dislocations in f.c.c. and b.c.c. metals under deformation is quantified using a recently developed coarse-grained (CG) atomistic model. Fundamental to the CG method is an atomistic field formulation that unifies the atomistic and continuum description of materials through an Irving-Kirkwood procedure in statistical mechanics. At a fraction of the cost of full molecular dynamics (MD), the CG model is shown to be applicable for calibrating the dislocation mobility law from the bottom up. In particular, using a modest computational resource, the mobility of extremely long dislocation lines (a length up to one micron and even above) in billion-atom material samples under deformation is predicted from the atomistic to the microscale. Results show that the dislocation mobility law in f.c.c materials is insensitive to its length. In contrast, in b.c.c. materials, when the length of a dislocation is at micrometer level and above, its mobility is found to be proportional to the line length, i.e., the longer the dislocation line, the faster it moves. With the atomistic information being retained, the CG simulations also provide an explanation for the observed different line length-dependent dislocation mobility in f.c.c. and that in b.c.c. materials: both short and long dislocation lines in f.c.c. materials migrate through gliding, while dislocation motion is mainly controlled by kink pair nucleation and diffusion in b.c.c. materials, and the length will be the limiting factor for kink pair density. In parallel, in order to validate the CG model, a limited set of full MD simulations has been also performed. The one-to-one comparison between CG and MD demonstrates that the CG simulation does retain the atomistic nature of dislocation dynamics, such as (i) the dislocation velocity-dependent line energy and core structure; and (ii) the phonon wave emission from a moving dislocation. In addition, the limitation, the future development, and the potential applications of CG will be also discussed in this presentation.

The discovery and optimization of phase-change materials is often a long and expensive process. Here we propose a simple computational method to determine the ideal phase change material for a given chemical composition. Using first-principles calculations, within a high-throughput computational framework, we are able to determine the ideal composition of a phase-change material between any two assumed phases. This ideal composition is determined by a calculation of the cofactor conditions which, when equal to 1, ensure the lowest possible interface energy between the phases during the phase transition. This ideal composition can then be targeted to produced materials with low mechanical failure rates due to the low interface energy. Here we will provide evidence of the effectiveness of our calculations for two well-known phase change materials in which we predict the ideal composition and compare it to experimental results.

Density Functional Theory (DFT) has the highest accuracy to cost ratio among all electronic structure methods and provides valuable insight in understanding and predicting a wide range of materials properties. Defects in crystalline solids, play crucial role in determining macroscopic properties of materials. The profound significance of defects underlies from the coupling between the discrete effects of the lattice, chemical effects of the core and the long-range effects of the elastic field. While DFT is capable of accurately describing the chemistry of the defect core, but are too complicated and expensive for defects with long range fields, methods capable for describing long-range fields rely on empiricism and lack fidelity. This poses an outstanding dual challenge to simulate defects from first principles. To overcome this, we develop a novel coarse-grained formulation of DFT. We employ Linear Scaling Spectral Quadrature method to solve for the electronic fields and develop a coarse-graining strategy based on updated Lagrangian method to describe the long-range elastic fields. We discuss a real space formulation based on high-order finite differences, local reformulation of electrostatics, reformulation of the atomic forces and a parallelization strategy based on domain decomposition method. The developed formulation is sublinear scaling with respect to the system size and can efficiently simulate defects in crystalline materials from first principles.
To demonstrate our framework, we use magnesium, the lightest among all commonly used structural metals, and has a high strength to weight ratio, and therefore has the potential to be significantly integrated into a wide range of applications. Existing magnesium alloys show low spall strength, which inhibits their applicability in dynamic environments. Magnesium alloys achieve their desirable mechanical properties via age hardening, and controlling the size of precipitates is necessary for spall strength. We use first principles calculations to understand the energetics of precipitation in magnesium-aluminum binary alloys and report the influence of deformations (which arise during processing) on the energetics of the precipitation process. We show strong dependence of concentration and volumetric strains on the cohesive energies of binary solid solutions and report a power law for these mixtures. We also show that the formation enthalpies of the solid solutions are dependent on the concentration and hydrostatic stress. The free energies of the solid solutions reported in this study are used to draw insight into the precipitation process. Since vacancies serve as nucleation sites for voids, which lead to spall, we also discuss the interaction of point defects (vacancy-solute, solute-solute and vacancy-vacancy) in magnesium-aluminum alloys using first principles calculations.