Symposium MT01—Advanced Atomistic Algorithms in Materials Science

Simulations based on modern orbital-free density (OF) functional theory (DFT) can be orders of magnitude faster than those performed with conventional Kohn-Sham DFT, largely because OF-DFT does not require any wavefunction operations. However, while OF-DFT is sufficiently accurate for many materials of practical relevance, it cannot yet replace wavefunction-based KS-DFT entirely, a longstanding goal. We report recent progress towards this goal, beginning with advances in the construction of the local pseudopotentials (LPSs) that are commonly employed to account for the electron-ion (screened nucleus) interaction. In particular, we highlight a method for constructing LPSs appropriate for both solid and liquid phases, commenting on OF-DFT simulations that recently confirmed the existence of a second excitation mode in the collective dynamics of liquid Sn. We also discuss efforts to overcome the main obstacle to universally accurate OF-DFT, which is determining the (noninteracting) electron kinetic energy without recourse to wavefunctions. These efforts range from careful study of the local, gradient-based kinetic energy density of nearly free electrons, including the derivation of new response functions for this purpose, to hybrid techniques that incorporate localized atom-centered density matrices within the usual OF formalism. We conclude by outlining recent updates to our widely used OF-DFT code, PROFESS.

DFT-based computational materials science literature is dominated by the use of plane wave codes, with a single such code being most widely used. In the modeling of bulk materials and interfaces, GGA functionals remain most widely used. There is a well-established narrative according to which such functionals (e.g. the popular PBE functional) cannot reproduce localized states, resulting from e.g. doping of oxides with alkali and alkaline earth atoms, and cannot correctly reproduce the qualitative nature of the bandstructure, e.g. ordering of p and d bands in some oxides. This is typically remedied by the use of the Hubbard correction (DFT+U). Such correction, however, often distorts the geometry and worsens the issue of local minima in density optimization. The perceived need to employ the Hubbard correction for semiconductor but not for metallic phases is also an issue, as it may be impossible to correctly model different phases with a single computational setup; this in turn makes more difficult the modeling of phenomena involving phase transitions.
We show that it is possible to obtain qualitatively correct bandstructures and localization of states with the PBE functional without the Hubbard correction when using localized basis sets in cases where +U is necessary with plane waves. This was observed in different codes, using pseudopotentials or all-electron calculations. Correct band ordering can be obtained (such as O(p) dominated valence band maximum in NiO) and localized states and spin ladders due to interstitial n-doping can be obtained in oxides (titania and vanadia). Moreover, without +U, qualitatively correct bandstructures can be obtained with PBE for metallic and semiconductor phases of VO2 simultaneously, which we could not obtain in VASP without or with +U with commonly used U values. Formation energies computed in such systems without +U are in good agreement with experiments (as probed by insertion voltages with oxide cathodes in metal ion batteries).
We highlight the fact that with plane wave basis sets, expansion errors are higher for localized than for delocalized states, which could be contributing to artificial delocalization. To fully nivelate this issue, unrealistically high plane wave cutoffs could be required.

Large-scale DFT calculations are very crucial to improve the predictive capabilities of computation based design of new materials in a variety of application areas. For example, accurate determination of dislocation core properties in metals and semiconductors, understanding ion conduction mechanisms and computing diffusivities in solid-state electrolytes, and large-scale bio-molecular simulations, all require the ability to conduct accurate and computationally efficient DFT calculations involving many thousands of atoms on both metallic and insulating systems. However, the solution to governing equations in DFT demands significant computational resources, and accurate DFT calculations are routinely limited to materials systems with at most a few thousands of electrons restricting the system sizes to a few hundred atoms. This computational complexity is of even bigger concern in the context of ab-initio molecular dynamics with long time-scales, and geometry relaxation of large-scale systems needing large number of relaxation steps.

To overcome the above limitations, an open-source code DFT-FE [1] based on adaptive finite-element discretization of DFT has recently been developed. This talk will present the methodological and algorithmic advances incorporated into DFT-FE, which has enabled systematically convergent and computationally efficient large-scale DFT calculations on materials systems with tens of thousands of electrons for both metallic and insulating systems, on massively parallel computing architectures, while allowing for arbitrary boundary conditions and complex geometries. Key computational ideas involved in the development of DFT-FE will be discussed, which include: (i) a unified real-space formulation for DFT treating both pseudopotential and all-electron calculations in the same framework [2]; (ii) mesh adaption for efficient finite-element discretization [2]; (iii) adaptive higher-order spectral finite-element framework in conjunction with Chebyshev acceleration [3], spectrum splitting based Rayleigh procedure combined with mixed precision arithmetic; (iv) configurational force approach for geometry optimization [4]. The numerical investigations conducted with DFT-FE on representative benchmark examples will be presented alongside performance comparisons with other widely used DFT codes like Quantum Espresso and ABINIT.


References

[1] GitHub Page: https://github.com/dftfeDevelopers/dftfe
[2] Motamarri P, Das S, Rudraraju S, Ghosh K, Davydov D, Gavini V. DFT-FE: A massively parallel adaptive finite-element code for large-scale density functional theory calculations. arXiv:1903.10959.
[3] Motamarri P, Nowak MR, Leiter K, Knap J, Gavini V. Higher-order adaptive finite-element methods for Kohn-Sham density functional theory, Journal of Computational Physics 2013;253:308-343.
[4] Motamarri P, Gavini V. Configurational forces in electronic structure calculations using Kohn-Sham density functional theory, Physical. Review B 2018;97:165132.

Two-dimensional (2D) van der Waals (vdW) heterostructures offer a fascinating platform to pursue fundamental science and novel device applications. Owing to quantum size confinement and reduced dielectric screening, electron-electron interaction and excitonic effect are prominent in 2D materials, which could dominate their electron and exciton dynamics. Recent experiments have revealed ultrafast charge (~100 fs) and energy (~1 ps) transfer dynamics in transition metal dichalcogenides (TMDs) heterostructures, which is surprising giving the strong electron-hole binding of the interlayer excitons. Most theoretical calculations to understand the experiments were performed based on local and semi-local exchange-correlation functionals, thus cannot capture the excitonic effect accurately. In this talk, I will introduce a first-principles method that combines non-adiabatic molecular dynamics (with the fewest switch surface hopping algorithm) and linear-response time-dependent density functional theory. Importantly, the method is formulated in the planewave bases and PAW pseudopotentials with range-separated hybrid functionals. As a result, the method can capture the excitonic effect accurately in extended 2D heterostructures. Using this method, we can shed light on the ultrafast charge transfer dynamics in TMD vdW heterostructures and elucidate the role of “hot” excitons in promoting ultrafast charge transfer despite the momentum mismatch in twisted heterostructures. We can also examine the properties of interlayer and intralayer excitons trapped by the moiré superlattices in the twisted heterostructures, and tune their properties and lifetimes by changing the twist angle, pressure and electric field.

Vibrational Raman spectroscopy is a widely used analytical tool in gas-phase and condensed-matter materials. Interpretation of the Raman spectra and detailed deduction of structures and/or conformations require theoretical simulations. However, the application of the present established computational Raman approaches to large-scale systems remains challenging, in a large part owing to their high computational demand. Here we propose a method that combine plane-wave pseudopotential method, time-dependent density functional perturbation theory and an elegant Lagrangian approach to circumvent key bottlenecks of Raman spectra simulations for large-scale systems. The method is validated for several molecular and solid systems, and exhibits excellent agreements with the experimental Raman spectra. We show that the method has a quadratic computational scaling with the number of electrons which can treat systems with more than 4000 electrons, opening doors for many large-scale Raman computations that are beyond the reach of previous approaches.

High-Throughput Quantum-Mechanics computation of materials properties by ab initio methods has become the foundation of an effective approach to materials design, discovery and characterization. This data driven approach to materials science currently presents the most promising path to the development of advanced technological materials that could solve or mitigate important social and economic challenges of the 21st century. In particular, the rapid proliferation of computational data on materials properties presents the possibility to complement and extend materials property databases where the experimental data is lacking and difficult to obtain.
The practical realization of these opportunities depends almost exclusively on the design of efficient algorithms for electronic structure simulations of realistic material systems beyond the limitations of the current standard theories. In this talk, I will review recent progress in theoretical and computational tools for data generation and advanced characterization, and in particular, discuss the development and validation of novel functionals within Density Functional Theory and of local basis representations for effective ab-initio tight-binding schemes.

Applications of these methods to a variety of systems (anomalous transport properties of 2D materials, topological insulators, minerals as novel thermoelectrics, etc.) will be showcased.

Magnesium alloys are promising as lightweight structural materials, electrodes in Mg batteries, and for use in medical implants.^1-3 As use of Mg is limited largely due to poor corrosion resistance a better understanding of these processes is required. Corrosion is an electrochemical process, involving transfer of electrons, thus a quantum mechanical theory is needed to model it. Density functional theory (DFT) has been used to study aspects of Mg corrosion, but its high cost limits the length and time scales which can be probed.^4,5 Density functional tight binding (DFTB) represents a cheaper, though approximate, alternative to full DFT that maintains a fully quantum mechanical description. In this work we present a DFTB model for Mg which can describe the many properties of the metallic system at a level similar to DFT

Compared to plane-wave DFT, our model introduces three approximations: we use the Harris-Foulkes functional up to first order (no self-consistency) ; we use a small (single-zeta) spd basis set; we treat electron exchange and correlation in an approximate way using the McWeda method.⁶

The validity of our DFTB model is demonstrated by comparison to plane-wave results for structural properties (e.g. volumes and bulk-moduli), electronic properties (density of states) and defect energies (self-interstitials and vacancies). Future work will focus on adding the necessary elements to investigate aqueous magnesium corrosion.


References
1. Li, N et al., J. Mater. Sci. Technol. 29, 489–502 (2013).
2. Saha, P. et al., Prog. Mater. Sci. 66, 1–86 (2014).
3. Abbott, T. B, CORROSION 71, 120–127 (2015).
4. Surendralal, S. el al., Phys. Rev. Lett., 120, 246801 (2018)
5. Yuwono, J. A. et al., J. Phys. Chem. C 120, 26922–26933 (2016).
6. Jelínek, P. et al., Phys. Rev. B., 71, 1–9 (2005)

Semi-empirical tight binding models are a powerful tool for the atomistic and molecular dynamics simulation of molecules and materials. The models can be derived from Kohn-Sham density functional theory via a hierarchy of approximations. Tight binding models with a non-orthogonal basis are relatively straightforward to parameterize because good initial guesses for the bond and overlap integrals can be obtained directly from density functional theory. The pairwise repulsive term, which is the crudest approximation in tight binding theory, can then be parameterized to reproduce total energies, forces, and stresses from either experiment or accurate density functional theory calculations. In practice, the pairwise term is used to capture all of the contributions to the potential energy and interatomic forces that are not explicitly included in the bond terms. I will describe our recent work on the development of accurate and transferable non-orthogonal tight binding parameterizations for organic molecules containing C, H, N, and O by the iterative optimization of the adjustable parameters in the models so that the errors in the binding energies and forces with respect to high-quality density functional theory calculations are minimized. We have applied the same general approach to the development of tight binding models for transition metals and alloys. The transferability of these models is generally good but with some exceptions. For instance, we find that the predicted energies of point defects are systematically too low with respect to density functional theory. Density functional theory calculations at several levels of theory suggest that the limited transferability of non-orthogonal tight binding models can be attributed to the use of pairwise repulsive interactions and a minimal basis of atomic orbitals.

A major challenge in predicting the properties of materials at realistic conditions is the accurate inclusion of finite temperature effects. Doing this on an ab initio level requires complex simulation protocols. These complex protocols, which often couple several specialised codes, make a quantitative description of error propagation and uncertainty quantification a critical issue.

To handle this high level of complexity we have developed the integrated development environment (IDE) pyiron – http://pyiron.org. The pyiron IDE combines a web based source code editor, a job management system for build automation, and a hierarchical data management solution. The core components of the pyiron IDE are pyiron objects based on an abstract class, which links application structures such as atomistic structures, projects, jobs, simulation protocols and computing resources with persistent storage and an interactive user environment. The simulation protocols within the pyiron IDE are constructed using the Python programming language.

In order to analyse the delicate interplay of uncertainties in our complex protocols, we introduce the concept of uncertainty phase diagrams. Based on this concept we model the convergence gradients of the contributing errors. It also provides a firm basis to automate the convergence process. Employing this approach we find that already the determination of the equilibrium lattice constant and bulk modulus require a careful analysis of the fitting of energy-volume curves, going beyond the consideration of standard convergence parameters like energy cutoff and k-points. Our investigations revealed that commonly used rules of the thumb for fitting ground state materials properties become invalid for high precision calculations, as the dominating sources of error change.

Tight binding models offer a compromise between the accuracy and transferability of DFT and the speed of empirical force field models. By explicitly computing the electronic structure of a given system, they can improve the understanding of the bonding in a given material. The use of a limited basis set of atom-centric orbitals and further approximations allow computationally expensive numerical integration to be replaced by a fast interpolation between pre-calculated values.
These models are promising when it comes to modeling corrosion or large scale defects (e.g. dislocation cores), where a good description of bonding plays a major role, but the system sizes needed are beyond the reach of conventional DFT. In particular, we would like to gain mechanistic insights about the effect of irradiation on the corrosion of zirconium alloy nuclear fuel cladding. To model the interaction of alloying and impurity elements with the damage resulting from neutron irradiation, a tight binding model capable of representing a range of impurity and alloying elements and atomic environments with good accuracy and transferability is needed. Here, we concentrate on the case of zirconium, however, the results can be generalised to most close-packed metals.
The majority of current methods for building tight binding models are heavily dependent on optimisation, making the addition of new elements challenging and the transferability limited [1,2].
In this work, we have followed a systematic, parameter-free method, starting from an appropriate choice of basis set and evaluating the effect of further approximations on the model. We have explored the two-centre approximation for both the pair potential and the components of the Hamiltonian matrix: crystal field terms and hopping integrals. While transferable and reasonably accurate corrections to the first two components are relatively straightforward, an equivalent correction for the hopping integrals is both theoretically and computationally challenging. We show that keeping the hopping integrals two-centre but adjusting their distance dependence causes significant transferability problems not only when describing other crystal structures, point defects or interstitial impurities, but also when simply considering one crystal structure in tension or compression.
The inclusion of three-centre integrals significantly improves the transferability of the model.

References
1. Schnell I, Jones MD, Rudin SP, Albers RC. Tight-binding calculations of the elastic constants and phonons of hcp Zr: Complications due to anisotropic stress and long-range forces. Physical Review B. 2006 Aug 8;74(5):054104.
2. Dufresne A, Ribeiro F, Tréglia G. How to derive tight-binding spd potentials? Application to zirconium. Journal of Physics: Condensed Matter. 2015 Aug 3;27(33):336301.

Understanding phonon-mediated heat transfer at the nanoscale remains an open and technologically-important challenge, with diverse applications across thermoelectrics, nanoelectronics, catalytic cells, and nanotheranostics. Such systems are often spatially inhomogeneous, necessitating theoretical descriptions with spatial resolution¹, and consist of non-crystalline disordered solids. Ab initio mathematical treatments developed for studies of lattice vibrations in perfectly ordered crystals, however, offer limited utility in investigations of disordered solids. These limitations include: (i) inseparability of diagonal (substitutional), off-diagonal, and ‘environmental’ disorder², and (ii) inconsistency in using perfect lattice states (i.e. phonons) as basis functions for the momentum representation of collective excitations³. Numerical calculations, using lattice and molecular dynamics, on simulated disordered solids have corroborated the breakdown of phonons as well-defined excitations (i.e. with mean free paths much larger than their wavelength), and have suggested an alternate classification of ‘propagons’, ‘diffusons’, and ‘locons’^4-5.

While these advances have provided insight in spatially homogeneous (bulk) thermal transport in disordered solids, they depart from the intuitive ‘kinetic theory’ picture and as such offer limited use in atomistic-to-continuum modeling of spatially-resolved transport. We treat the medium as a continuum capable of supporting approximate plane-wave excitations, experiencing scattering due to the atomistic disorder, in much the same way external plane waves (e.g. x-rays) would³. We compute these scattering rates, first proposed by Ziman and Morgan⁶, from first principles using Fermi’s golden rule. In particular, the formalism requires access to the real-space harmonic and anharmonic interatomic force constants which, we will show in this talk, are extracted at finite temperature by sampling the canonical phase-space using Born-Oppenheimer ab initio molecular dynamics (BO AIMD)⁷. We demonstrate this formalism by computing the spatially-resolved thermal transport across crystalline/amorphous silicon (Si/α-Si) heterostructures with micron dimensions, length scales typically not accessible to (classical) molecular dynamics.

¹ G. Varnavides, A. S. Jermyn, P. Anikeeva, and P. Narang, arXiv:1811.01059v1 (2018).
² S. Ghosh, ,P.L. Leath, and M.H. Cohen, Physical Review B, 66 (2002), 214206.
³ J. M. Ziman, Cambridge University Press Archive, 1979.
⁴ P. B. Allen, J. L. Feldman, J. Fabian, and F. Wooten, Philosophical Magazine B 79, 1715 (1999).
⁵ H.R. Seyf et al, npj Computational Materials 3, 49 (2017).
⁶ G. J. Morgan, Journal of Physics C: Solid State Physics 1, 347 (1968).
⁷ O. Hellman and I. A. Abrikosov, Physical Review B 88 (2013), 144301.

Computing the density matrix for a large Hamiltonian is the most expensive part of many electronic structure calculations. At zero temperature, the density matrix is a projection matrix constructed from the outer product of the lowest occupied eigenstates, satisfying idempotency and trace conditions. In this talk, we present a method combining two well-known algorithms – the Kernel Polynomial Expansion (KPE) method and Spectral Projection Purification – which allows us to obtain idempotent density matrices of large systems while relying only on sparse matrix-vector multiplication kernels. An implicit representation of the density matrix is obtained by combining the two aforementioned algorithms, which is then sampled using random vectors. This scheme allows us to compute the density matrix at linear computational cost in terms of floating-point operations. Additionally, the method is independent of the underlying basis set, and does not require the density matrix to be sparse. Within the context of density-functional based tight-binding (DFTB), we present results of self-consistent density matrix calculations of systems containing 10^7 atoms. As an application of the method, we present results on time-dependent electron emission from dielectric surfaces under AC voltages in the kHz frequency range. The added challenge here is that while oscillations of the electronic system are in the THz or higher range, added perturbations due to AC voltage are in the kHz frequency range. As such, direct integration of the evolution of the electronic system requires long time-scales. We solve this challenge first using a semi-classical approximation in order to avoid temporal integration. Additionally, we show that a kernel polynomial expansion of the time-evolution operator can be used to allow direct integration of the quantum-mechanical evolution under an external AC voltage. This demonstration of KPE for both the electronic structure and the time integration of the system provides a computationally efficient framework for an inherently multiscale electronic structure problem, with implications for a variety of challenging problems.

Strong (10¹⁰ V/m) electric fields can be used to trigger chemical processes with extreme precision by selectively stabilizing or weakening bonds to initiate reactions which are otherwise slow or do not proceed at all. The ability to manipulate electric fields to tailor and stimulate bond-breaking events is therefore a powerful experimental control knob, but one whose effects are difficult to predict due to a lack of suitable tools to probe its associated atomic-scale mechanisms. Atom Probe Tomography (APT) directly probes field-induced bond-breaking mechanisms, but its results must be interpreted through theoretical models, which are currently overly simplistic. Here we introduce a novel approach, which we term the Generalized Dipole Correction (GDC), that enables the direct study of ultra-high fields effecting bond-breaking and desorption at the level of single atoms using Density Functional Theory (DFT). As a prototype application, we consider field evaporation from a kinked W (110) surface. We reveal two qualitatively different competing mechanisms that can be switched by the applied field.

The Einstein-de Haas (EdH) effect, where the spin angular momentum of electrons is transferred to the mechanical angular momentum of atoms, was established experimentally in 1915. While a semi-classical explanation of the effect exists, modern electronic structure methods have not yet been applied to study this phenomenon. In this talk we investigate its microscopic origins by means of a non-collinear tight-binding model of an O₂ dimer as implemented in Ref. [1]. This includes the effects of spin-orbit coupling, coupling to an external magnetic field, and non-collinear Stoner exchange. By varying an external magnetic field in the presence of spin-orbit coupling, a torque can be generated on the dimer, validating the presence of the EdH effect. Avoided energy level crossings and the rate of change of magnetic field determine the evolution of the spin. We find also that the torque exerted on the nuclei by the electrons in a time-varying B field is not only due to the EdH effect: other contributions arise from field-induced changes in the electronic orbital angular momentum and from the direct action of the Faraday electric field associated with the time-varying magnetic field.

[1] T. Wells, A. P. Horsfield, W. M. C. Foulkes, and S. L. Dudarev. The microscopic Einstein-de Haas effect. arXiv:1905.10449, May 2019.

Computational materials research has become increasingly vital in probing the properties of crystalline materials, especially in screening materials at a large scale to accelerate material discoveries for a wide range of applications.^[1-2] A routine operation for calculating the properties of crystalline materials is the evaluation of integrations over the Brillouin zone, which can be approximated by sampling the Brillouin zone at a set of points known as k-points. Many popular computational materials software packages generate k-points using the traditional Monkhorst-Pack (MP) scheme,^[3] which creates k-point grids that are regular and aligned with the reciprocal lattice vectors. We have previously demonstrated that much more efficient lattices can be generated if the Monkhorst Pack approach is generalized by relaxing the requirement that the grids be aligned with the reciprocal lattice vectors,^[4] which has since been confirmed by other researchers.^[5] This approach can roughly double the speed of well-converged calculations on crystalline materials, providing strong incentives to find an efficient algorithm for fast generation of generalised Monkhorst-Pack grids.

Here, we present algorithms for rapidly generating generalized Monkhorst-Pack grids. To better serve the computational materials community, we provide three implementations of these algorithms in various languages and formats:
1. K-Point Grid Server:^[4] a ready-to-use online application with a pre-generated database for fast search.
2. K-Point Grid Generator: an open-source, stand-alone application for runtime environments that might not have an internet connection.
3. kplib: a lightweight open-source C++ library that does not make use of databases, for integration with computational materials software packages.

Benchmark results are provided to demonstrate the speed of these algorithms and the quality of generated grids.

References
1. A. Jain, S.P. Ong, G. Hautier, W. Chen, W.D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K.A. Persson, APL Materials 1, 011002 (2013).
2. S. Kirklin, J.E. Saal, B. Meredig, A. Thompson, J.W. Doak, M. Aykol, S. Rühl, and C. Wolverton, Npj Computational Materials 1, (2015).
3. H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
4. P. Wisesa, K.A. Mcgill, and T. Mueller, Physical Review B 93, (2016).
5. W.S. Morgan, J.J. Jorgensen, B.C. Hess, and G.L. Hart, Computational Materials Science 153, 424 (2018).

High-accuracy quantum simulations of materials is becoming increasingly important for the development of new devices for electronic, thermoelectric, and energy applications. While these simulations tend to scale exponentially with the size of the system on classical computers, it has been shown that they can be made to scale polynomially on quantum computers. Performing dynamic simulations on quantum computers, therefore, can open new doors by allowing systems of larger (and more relevant) sizes to be simulated. In this work, we simulate photonic control of emergent magnetism in a two-dimensional material on Rigetti’s Aspen quantum processor. Results from the quantum computer are found to be in good qualitative agreement with theoretical ground-truth results, and in even better agreement with results from classically simulated noisy qubits. Furthermore, we map out the feasibility bounds of such dynamical simulations on state-of-the-art, near-term and near-future quantum computers.

This work was supported as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0014607.

Determining the atomistic details of nanoscale structures is a fundamental problem. Although there are both experimental and computational methods to determine these nanoscale structures, they both possess limitations. We develop the FANTASTX code (Fully Automated Nanoscale To Atomistic Structure from Theory and eXperiment) to overcome the limitations of either by combining both experimental and computational data. We demonstrate the effectiveness of FANTASTX by determining the structures of nanoparticles and solid interfaces from x-ray and electron microscopy data combined with atomistic and first principles energies, using multi-objective optimization and a variety of canonical and grand canonical sampling algorithms.

The properties of molecules and materials are ultimately determined at the atomistic level. On the other side, their behavior at macroscopic scales is far more relevant to model for example chemical reactor processes and to determine promising changes to optimize catalysts, additives, and operating conditions in industrial processing plants. The multi-scale simulation platform ReaxPro aims to provide a customer-ready solution by combining several well-established tools and simulation programs such as the Amsterdam Modeling Suite, kinetic Monte Carlo methods, and computational fluid dynamics to bridge the many orders of magnitude in scale between atomistic quantum chemistry models and macroscopic chemical reactor design.
Such multi-scale simulations will not only significantly reduce experimentation costs but also allow for an in-silico high-throughput screening and predictions to guide experiments. In addition, valuable insights into the mechanisms at each scale can be gained in order to help researchers to improve reactor designs at different levels. The multi-scale material and process optimization of ReaxPro can also be applied in other fields of chemical industries such as the application-specific design of new polymers, photovoltaics, and batteries.

We have developed a module for calculating the mobility and conductivity of semiconducting materials based on the iterative method as proposed by Rode [1]. Various physical properties calculated through the density function theory (DFT) are used as inputs. We show the use of this module for Cadmium Sulfide (CdS) and Zinc Selenide (ZnSe). Since this strategy can take into consideration the inelastic scattering of the carriers due to the phonons, this module should serve a better objective than those based on the relaxation time approximation.
References
[1] D. L. Rode, Semiconductors and Semimetals (Academic Press, New York, 1975), Chapter 1.

We introduce an orbital free electron density functional approximation based on alchemical perturbation theory[1]. Given convergent perturbations of a suitable reference system, the accuracy of popular self-consistent Kohn-Sham density functional estimates of properties of new molecules can be systematically surpassed - at negligible cost.

For example, using the CCSD solution for N₂, APDFT calculated properties of CO are more accurate than PBE already at 1^st order (energies and dipole moments) and 2^nd order (quadrupole-moments and forces).

The associated energy functional is an approximation to the integrated energy derivative, requiring only perturbed reference electron densities: No self-consistent field equations are necessary to estimate energies and electron densities. Instead, our approach relies on the electron density response w.r.t nuclear charges and treats changes of nuclear charges at any sites as perturbations to the system. We show that the resulting expansion in perturbation orders converges quickly by analytical proof for the hydrogenic atom and for any free atom. Numerical convergence is shown for alchemical perturbations of H₂, N₂, and benzene[1, 2].

APDFT based estimates of the electron density of a target molecule are obtained for the same perturbations.
Estimated electronic ground state properties considered include covalent bonding potentials, atomic forces, as well as dipole and quadrupole moments.

We demonstrate the applicability of APDFT to scan a large chemical space by using 13 CCSD single point calculations to predict changes in energy and forces by BN-doping a C₂₀ fullerene. This allows us to give total energy estimates for all 3.1 million uniquely BN-doped fullerenes.

APDFT is widely applicable to any level of theory that makes electron densities available and allows to assess a combinatorial number of molecules with one fixed set of calculations rather than calculating molecules one-by-one.

If the perturbation series converges and if the reference level of theory is of sufficient quality, APDFT represents a systematically improving DFT approximation of hierarchies of accuracy.


[1] G. F. von Rudorff, O. A. von Lilienfeld, Alchemical perturbation density functional theory (APDFT), arXiv:1809.01647, 2019.
[2] G. F. von Rudorff, O. A. von Lilienfeld, to be submitted, 2019.

Density Functional Theory (DFT) has become a ubiquitous and vital tool in a wide range of materials physics problems. It occupies an important space in the tradeoff between accuracy and computational cost - able to achieve high levels of accuracy and maintain computational feasibility. Within the DFT formulation this tradeoff is adjusted through the exchange-correlation energy functional. Advanced functionals can be substantially more accurate but can also be far more expensive computationally. Recent work has attempted to push back the optimal performance tradeoff curve through machine learned functionals trained on high level exchange-correlation energies, such as those from hybrid functionals or CCSD(T) [1,2]. A common theme among these is to assume a local model, where the exchange-correlation energy is a sum of local contributions.

While recent efforts have been primarily driven by neural network models to predict these contributions, we consider an alternative procedure wherein the electronic densities are projected onto a Gaussian basis to drastically reduce their dimensionality. An ordinary least squares problem is then solved to regress the exchange-correlation energy. In this context, we examine the limitations of a local model: we demonstrate that our model can accurately reproduce LDA functional data and explore effects on accuracy when trained on hybrid functionals. We further consider the implications of this model in its capacity to accurately and efficiently predict material properties.


[1] Lei, Xiangyun, and Andrew J. Medford. "Design and Analysis of Machine Learning Exchange-Correlation Functionals via Rotationally Invariant Convolutional Descriptors." arXiv preprint arXiv:1901.10822 (2019).
[2] Kolb, Brian, Levi C. Lentz, and Alexie M. Kolpak. "Discovering charge density functionals and structure-property relationships with PROPhet: A general framework for coupling machine learning and first-principles methods." Scientific reports 7, no. 1 (2017): 1192.

The development of reliable, yet computationally-efficient interatomic forcefields is key to ensure the accuracy of classical molecular dynamics simulations. However, parameterizing new forcefields is notoriously challenging, since the high number of parameters renders traditional optimization methods inefficient or subject to bias. Here, we present a new parametrization method based on machine learning, which combines ab initio molecular dynamics simulations and Bayesian optimization. By taking the examples of silicate and chalcogenide glass, we show that our method yields new some empirical interatomic forcefields that offer an unprecedented level of agreement with ab initio simulations. This method offers a new route to efficiently parametrize new interatomic forcefields for disordered solids in a non-biased fashion.

Recent developments in neural network interatomic potentials have provided a pathway for reliable large-length and long-time scales atomistic simulations. However, their transferability is questionable because of incompleteness of the training dataset that relies on structure-energy (and/or forces) calculated using density functional theory (DFT), leading to inaccurate extrapolations or interpolations of the neural network. In this talk, we present two efficient methods for generating configurational spaces for training dataset: (1) molecular dynamics simulations with large time steps and (2) random walks on uniform grids. Both methods create chaotic displacements of atoms that move them far from equilibrium positions or trajectories, effectively enlarging the configurational space. We will demonstrate the effectiveness of the methods using examples of transferable neural network interatomic potential for: 1) the Si(001) surface and 2) defect-laden, two dimensional, hexagonal boron nitride (dh-BN). Both potentials are able to predict total energy within few meV of that obtained from DFT. Molecular dynamics simulations of the Si(001) surface for a wide range of temperature that show that the asymmetric, buckled structure of the surface dimer exists at all temperatures, but their increased flipping rate leads to a larger occurrence of the symmetric configuration, shedding light on the semiconducting-conducting transition of the surface. In the case of dh-BN the method allows simulations of grain-boundary structures consisting of thousands of atoms whose local configuration is not included in the training set. We compare our results with available data and discuss the pros and cons of the approach.

* Work supported in part by DOE Grant No. DOE-DE-FG02-07ER46354.

High quality force field models are not available for many phase change materials such as Ge₂Sb₂Te₅, which has diverse applications in neuromorphic computing devices and non-volatile memory. Designing semi-empirical force fields for these systems is nontrivial and also doing ab-initio molecular dynamics (MD) simulations is not scalable for large systems, which is important to amorphization and recrystallization in these systems. We have developed a neural network (NN) based force field model for GeSe₂ and Ge₂Sb₂Te₅, which is trained using ab-initio MD simulations data to directly produce atomic forces at density functional theory (DFT) level accuracy. The neural-network architecture and input parameters of the feature vector for each atom type are optimized for training using Bayesian optimization so as to design a model with higher accuracy but with minimum complexity. Special care is taken to design radial and angular feature vectors that obey translational and permutational invariance and rotational covariance for forces. The accuracy of the NN force field model is validated by doing a MD simulation using NN force field model for up to 100,000 atoms and by computing various structural and dynamical properties such as radial distribution function, bond angle distribution, neutron structure factor and phonon density-of-states.

Acknowledgment
This work was supported as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under award number DE-SC0014607. The simulations were performed at the Argonne Leadership Computing Facility under the DOE INCITE and Aurora ESP programs and at the Center for High Performance Computing of the University of Southern California.

Deep Learning has emerged as a powerful tool for accelerating simulations of molecular and condensed matter by learning complex representations of the potential energy directly from atomic positions without the need for handcrafted descriptors. We present a deep learning approach aimed at molecular dynamics that directly predicts atomic forces instead of total energy. The proposed force-field is rotationally covariant, handles multiple species naturally, and bypasses the need for computing expensive derivatives of an energy-based model by predicting atomic forces directly from atomic position coordinates. The proposed algorithm achieves very high accuracy on a data set of organic molecules while outperforming energy-based models with respect to computational efficiency, both at time of training and inference. We show that our proposed method can accurately reproduce thermodynamic and kinetic materials properties from Molecular Dynamics simulations and benchmark the efficiency of the direct force prediction against energy-based models.

We present a general framework to generate machine learned interatomic poten-
tials from ab initio density functional theory (DFT) databases. In our framework,
atomic structure information residing in the configuration space is projected on to
a space of descriptors that can systematically capture local many-body correlations.
Specifically, we have included two-body, and a variety of three and four-body cor-
relations. These correlations are realised using Chebyshev basis functions, but the
framework is general and other basis functions can be used. To generate the inter-
atomic potential, we have explored linear as well as non-linear regression techniques
with different regularisations. The different correlation terms were grouped accord-
ing to the number of atoms, and according to the bonds present in them. These
variables control the dimensionality of the parameter set of the interatomic poten-
tial. We applied this framework to generate an interatomic potential for germanium
using DFT (hybrid functional) database consisting of about 1.5 million atoms en-
vironments. This database includes structures with different crystal symmetries as
well as liquid and amorphous structures. Different potential models were generated
by using different types of correlation descriptors and the effect of different correla-
tions were compared by analysing root mean square error as well as with histogram
of errors between the model and DFT. In summary, we have developed a framework
for building potentials with tuneable accuracy and shown its application to the case
of germanium. The model is implemented as an extension to LAMMPS molecular
dynamics code and can be used in large-scale simulations.
This work performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

Neural-network force field (NN-FF) has decome an intensedly studied topic in recent years. In this talk, I will present our recent work in developing NN-FF and using it in molecular dynamics simulations. In our approach to train the NN-FF, not only we use the atomic forces from ab initio calculations, we also use decomposed atomic energies from density functional theory (DFT) calculations. This increases the number of data set, also makes the training process easier. We have developed Si and Na NN-FF in this approach. We also developed multiple element force fields like Al-H and Fe-H. We have used the NN-FF to study the crystal Si growth from melt Si. This NN-FF shows similar behavior as that of the direct DFT calculations, but it allows us to study much larger systems. I will also present our approaches to deal with long range Coulomb interactions for systems with ionic charges.

Defining interatomic potentials as a high-dimensional fit of the reference (usually density functional theory) potential energy surface by taking advantage of regression methods from machine learning has recently become a powerful approach. However, because of their variational freedom, such potentials require large fitting datasets, which are typically developed using large amounts of manual selection and tuning of configurations by the researcher. To simplify this process, we present a highly automated iterative method where a preliminary Gaussian approximation potential (GAP) is used to carry out random-structure searches, and selected configurations from the searches are used to fit the next iteration's GAP. The resulting set of configurations constitutes an effective fitting database, spanning a wide range of geometries with a minimal number of reference energy evaluations. We test the method on several elemental and multi-component materials with different bonding types, from insulators to metals, including boron with its complex icosahedral α-B₁₂ structure and the phase-change material GeTe. We show how the process converges in a few iterations, and how the resulting potentials reproduce the reference DFT values on a number of bulk and defect properties.

MRI contrast agents based on iron oxide nanoparticles are advantageous alternatives to Gadolinium-based constrast agents in terms of toxicity, safety, and biocompatibility. However, it remains necessary to develop novel iron oxide nanoparticles able induce better T₁- and T₂-weighted relaxivities. Specifically, the manipulation of the iron oxide particle geometry should help increase the effective surface while provoking larger local disturbances in the magnetic field with a concomitant enhancement in relaxivities. Thereby, we synthesized spherical, cubic, hexagonal, and star-shaped iron oxide nanoparticles, characterized their structural and magnetic properties, and measured their T₁- and T₂-weighted relaxivities. The results show a significant relaxivity enhancement for the non-spherical particles, thus confirming the hypothesis that geometry has a large impact on the nanoparticles ability to locally disturb the the random arrangement of the spins. We will discuss the nanoparticles property-function relationships and their potential as MRI dual T₁-T₂ contrast agents.

Strongly anharmonic solids have been widely sought, ranging from thermoelectrics to multiferroics and to shape-memory materials, but the effects of strong anharmonicity to their thermodynamic and transport properties are less understood. In this work, we extend our earlier crystal graph convolutional neural network framework, and explore an alternative Bayesian description for the vibrational modes and their transport in strongly anharmonic solids exhibiting soft modes. We demonstrate our theory on a one-dimensional toy lattice, including a probabilistic extension to the conventional phononic dispersion relationship, as well as its consequences to lattice conductivity, other thermal properties, and phase transition. This Bayesian framework also provides naturally an uncertainty quantifying scheme for these physical quantities, including those uncertainties inherent for meta-stable and partial-crystal-partial-glass solids. Our theory is then extended to two halide perovskites: all inorganic CsPbI₃ and hybrid inorganic-organic MAPbI₃, and our Bayesian predictions are compared with existing experiments of heat capacity and lattice conductivity. More interestingly, some anomalous phononics observed from experiments could be explained within this framework.

Modeling the thermodynamic and kinetic properties of materials commonly requires taking into account both vibrational and compositional degrees of freedom.
While electronic structure calculations strictly speaking can provide the necessary information, they are computationally too demanding for many practical purposes.
One therefore requires models that can reproduce the relevant parts of the potential energy landscape at a much smaller computational cost.
Here, we will present a set of software packages that we developed over the course of the last years that allow one to rapidly construct and sample models that map vibrational (hiPhive [1]) and/or compositional degrees of freedom (icet [2]).
These tools can be integrated in a homogeneous Python based workflow, which combines structure generation, electronic structure calculations, model construction, thermodynamic sampling, and subsequent analysis.
Example applications include the prediction of alloy phase diagrams, temperature driven phase transitions, thermal conductivities, and ordering phenomena.

[1] F. Eriksson, E. Fransson, and P. Erhart, Adv. Theory Simul. 1800184 (2019)
[2] M. Ångqvist, W. A. Muñoz, J. M. Rahm, E. Fransson, C. Durniak, P. Rozyczko, T. H. Rod, and P. Erhart, Adv. Theory Simul. 1900015 (2019)

Propelled partly by the Materials Genome Initiative, and partly by the algorithmic developments and the resounding successes of data-driven efforts in other domains, machine learning (ML) strategies are beginning to take shape within several subfields of materials science[1-4]. One area of enormous importance is atomic-level materials phenomena, which spans numerous fields from electronic structure theory to mechanical failure, and has been dominated by either quantum mechanics (QM) based methods [4-7]—which are time-intensive, but accurate and versatile—or semi-empirical/ classical methods[8,9]—which are fast but are significantly limited in veracity, versatility and transferability. ML methods have the potential to bridge the chasm between the two extremes and can combine the best of both worlds. We have created a platform for the rapid prediction of properties such as potential energy, atomic forces, stresses, charge density, and the electronic density of states. Our ML-models are trained on accurate QM reference data, and can reproduce the QM results with the same level of accuracy but several orders of magnitude faster [10-15]. The ML models can also be progressively improved in quality by periodically (or on-demand) exposing them to fresh QM data in regions of poor performance, a feature that currently is either impossible or daunting with modern empirical/classical methods. Our ML models may thus be used to perform large-scale and/or long-time simulations of important materials phenomena previously beyond the reaches of QM based methods. Here, we demonstrate the power and versatility of this new platform in correctly capturing electronic, thermodynamic, mechanical, and diffusive properties for a variety of systems, with the hope of ushering in a new era of atomic-level understanding of materials.

References:

[1] Chan, H.; Cherukara, M.; Narayanan, B.; Loeffler, T.; Benmore, C.; Gray, S.; Sankaranarayanan, S. Nat. Comm. 2019, 10, A397
[2] Behler, J.; Parrinello, M. Phys. Rev. Lett. 2007, 98, 146401
[3] Botu, V.; Ramprasad, R. Int. J. Quant. Chem. 2015, 115, 1074–1083
[4] Kolb, B.; Lentz, L.; Koplak, A. Scientific Reports, 2017, 7, A1192
[5] Jones, R. O. Rev. Mod. Phys. 2015, 87, 897-923
[6] Hohenburg, P.; Kohn, W. Phys. Rev. 1964, 136, B864-871
[7] Kohn, L., W; Sham, K. Phys. Rev. 1965 , 140, A1133-A1138
[8] Daw, M. S.; Baskes, M. I. Phys. Rev. B 1984 , 29 , 6443–6453
[9] Hrennikoff, A. IABSE, 1968
[10] Botu, V.;Batra, R.;Chapman, J.;Ramprasad, R. Jour. Phys. Chem. C, 2017, 121, 511–522
[11] Chandrasekaran, A.; Kamal, D.; Batra, R.; Kim, C.; Chen, L.; Ramprasad, R. npj Computational Materials, 2019, 5, 22
[12] Batra, R.; Chandrasekaran, A.; Chapman, J.; Kim, C.; Chen, L.; Huan, T.; Ramprasad, R. Jour. Phys. Chem. C, 2019
[13] Huan, T.; Batra, R.; Chapman, J.; Kim, C.; Chandrasekaran, A.; Ramprasad, R. (manuscript under review)
[14] Chapman, J.; Batra, R.; Ramprasad, R. (manuscript in preparation)

Dislocations, which are line defects in solid, are a dominant factor of plastic deformation of solid. For advanced material design, it is very important to clarify the details of dislocation mechanics and dynamics. The energetics of dislocation mainly consists of the following two-part: atomic dislocation core and long-range elastic field. The dislocation core structure and behavior determine the mobility of individual dislocations and the condition of reactions between them. The elastic interaction determines the structure of mesoscopic aggregate of dislocations, such as dislocation cell walls. This structure is directly connected macroscopic mechanical behavior of crystals, such as work hardening.
To investigate the long-range elastic interaction, discrete dislocation dynamics (DD) modeling is a powerful tool, because it can treat micrometer scale in a realistic time scale. DD modeling is based on dislocation elastic theory. So, DD modeling requires the mobility and condition of the reaction of dislocations as input parameters. To investigate the dislocation core structure, the atomistic modeling is very suitable, because dislocation core structure and behavior are severely affected by the discreteness of atomic spacing. High accuracy density functional theory (DFT) calculations have established a leading position in atomistic modeling in decades. Due to the high-cost calculation of DFT, the number of atoms in DFT calculation is limited to a few hundreds, practically.
One solution to link the DD modeling and DFT atomic modeling is an atomic potential based on artificial neural network (A) framework [1]. By the universal approximation theorem, the ANN potential can compute any function. Therefore, with enough DFT training data, ANN potential becomes sophisticated replica potential of DFT calculation. In this work, we construct the atomic ANN potential for iron (Fe).
For DFT training date, we prepare the 5173 atomic structure of Fe calculated by Quantum ESPRESSO package [4]. For crystal structure data, we select the FCC, HCP, BCC and simple cubic structure. For defect data, we select the vacancy, divacancy, self-interstitial (T-site, O-site, [100], [110], [111]-dumbbell) and (100), (110), (111) and (112) surfaces model. To predict the exact dislocation core structure, we also prepare generalized stacking fault energy (GSF) surface on (110) and (112). The Chebyshev descriptor developed by Artrith et.el. is used for the input descriptor of the atomic environment [2,3]. For neural network training, we use the aenet package by Artrith et el [3].
By using ANN potential, We investigate the a/2<111> screw dislocation core structure in BCC iron, where a is lattice constant. The screw dislocation core predicted by constructed ANN potential has nondegenerate compact structure. We also evaluate the Peierls barrier of a/2<111> screw dislocation by using the nudged elastic band (NEB) method. The Peierls barrier is estimated as 36.6 meV in units of length of Burgers vector. These results are a good agreement with previous DFT calculation [5].
References:[1]J. Behler and M. Parrinello, Phys. Rev. Lett., 98 (2007) 146401. [2]N. Artrith, A. Urban, and G. Ceder, Phys. Rev. B 96 (2017) 014112. [3]N. Artrith and A. Urban, Comput. Mater. Sci. 114 (2016) 135. [4]P. Giannozzi, et al., J. Phys.:Condens. Matter., 21 (2009) 395502. [5]M. Itakura, H. Kaburaki, and M. Yamaguchi, Acta Mater. 60 (2012) 3698.

Light interstitial atoms such as hydrogen, helium, oxygen, carbon and nitrogen are known to have profound effects on the material properties experimentally. Ab initio and atomistic methods are employed to address their equilibrium segregation, effect on thermally activated processes, as well as non-equilibrium effects. A conceptual framework based on Gibbs adsorption isotherm, but generalized to the saddle state, with activation strain-volume/entropy/excess, is developed to interpret the simulation and experimental results.

The approach consisting in coupling spin and molecular dynamics is an increasingly popular
numerical method allowing to investigate the reciprocal influence of lattice vibrations and magne-
tization dynamics. Through the recently added SPIN package, LAMMPS now allows to perform
such simulations at unprecedented accuracy and scale. A brief overview of this implementation
and of some of its recent improvements will be first presented.
Although this implementation greatly improved the available algorithms for this methodology,
the lack of magneto-elastic potentials remains a strong limitation to its development. Indeed, the
most commonly used approach consists in overlaying the magnetic interactions on top of a purely
mechanical potential (typically, a magnetic exchange interaction on top of an EAM potential), all
those interactions being parametrized separately. As the magnetic interactions generates rather
small forces compared to the mechanical potential, this approach proved to be sufficient to obtain
numerous first order results.
However, it remains a rather crude approximation and cannot be used to recover fine exper-
imental results. In this presentation, a new methodology based on high-fidelity machine-learned
inter-atomic potential approaches and enabling to circumvent this difficulty will be presented.

Interactions between heterogeneous systems are challenging to model with high fidelity at the atomic scale across significant length scales. This presentation describes recent developments of third-generation charge optimized many-body (COMB3) potentials to enable metal-carbon interactions in classical molecular dynamics simulations. The potentials are applied to investigate the interaction of carbon materials, including metal-carbide-derived-carbons (CDCs) with titanium and aluminum, and the interactions of graphene and buckyballs with copper and aluminum surfaces. The results provide new insights into the bonding character at these interfaces and the role of defects, bond-angle, and temperature on interatomic interactions.

First principles molecular dynamics efficiently calculates the internal energy of metals, both in liquid and solid form. We show that the liquid state correlation functions that are a byproduct of the simulation allow a similarly efficient calculation of absolute entropy based on a simple variant of the classical "S2" method. The method is quite accurate for nearly free eletron metals, as demonstrated by comparison with experimental data for elemental liquid metals. For liquid alloys, hybrid Monte Carlo/molecular dynamics is utilized to enhance ensemble sampling. We combine the energy and entropy to calculate absolute free energies in liquid metals in order to explain the occurence of phase separation in liquid Li-Na compared with the eutectic in Na-K.

Chemically complex materials, namely materials comprised of multiple chemical components, are playing an increasingly important role in material applications including high entropy alloys (HEAs) for structural materials and complex mixed ion perovskite materials for solar cells as well as optoelectronic devices. The microstructures of these materials, for example, chemical short-range orders or phase segregations are critical for the material performance; however, characterization of microstructures of these chemically complex materials poses grand challenges to both experimentalists and theorists. From computer simulation perspective, exhaustive sampling of microstructures is computationally prohibitive for density functional theory (DFT) calculations owing to both system size and computation time limitations. Herein, we demonstrate that by utilizing neural network energy predictors trained from numerous images from DFT calculations, a 10⁵ computational speedup can be obtained by implementing the neural network predictor into the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), while retaining high fidelity to DFT calculations. Hence, the neural network energy predictor allows us to perform exhaustive sampling of systems of thousands of atoms of CoCrFeNi HEAs and FA_xMA_1-xPbI_yBr_3-y perovskite materials using Monte Carlo in combination with molecular dynamics simulations, thereby providing atomistic insights into microstructures and chemical short-range orders of these chemically complex materials with DFT accuracy.

New dielectric polymer materials with high energy density while maintaining low loss at high applied electric fields are needed for increased energy and power density required in modern electronics. We propose a polarizable reactive molecular dynamics-based framework that offers near-quantum accuracy for estimating dielectric properties of various polymer systems within a small fraction of quantum-calculation time. We use this framework to demonstrate the essential role of morphological complexity that governs dielectric constant of industrially relevant model polymer systems. Our framework supports high dielectric performance of PDTC-HK511 and TDI-EDR148 polymer systems, which have been identified through rational co-design strategy based on advanced computational screening procedures. The proposed framework can augment screening procedures for identifying polymer dielectrics that are significantly better than the current state of the art.

Gaussian process regression was previously demonstrated to provide quantitative Bayesian estimates of prediction uncertainty, enabling on-the-fly active learning of force fields in molecular dynamics of rare events¹. However, the time complexity in this method of prediction grows linearly with the training data size. We build on earlier attempt to exploit the structure of the N-body kernel function, to decompose the prediction expression into contributions that can be interpolated in order to approximate atomic forces, without loss of accuracy compared to the full regression². Crucially, this approach no longer depends on training size and accelerates the Gaussian process prediction by several order of magnitudes. We present a method of mapping (1) force fields from the underlying kernel, and (2) predictive variance (uncertainty), combined with dimension reduction approaches. This method can be combined into an end-to-end on-the-fly Bayesian active learning workflow, accelerating prediction of both the force values and variances, with ab-initio accuracy and rigorous Bayesian uncertainty quantification. We demonstrate this workflow on large-scale molecular dynamics simulations of properties of 2D stanene and related structures, including adatom migration and melting transition.

[1] Jonathan Vandermause, Steven B Torrisi, Simon Batzner, Alexie M Kolpak, and Boris Kozinsky. Submitted. “On-the-Fly Bayesian Active Learning of Interpretable Force-Fields for Atomistic Rare Events.” arXiv preprint arXiv:1904.02042

[2] Glielmo, Aldo, Claudio Zeni, and Alessandro De Vita. "Efficient nonparametric n-body force fields from machine learning." Physical Review B 97.18 (2018): 184307.

Ability to compute reaction rates from free energy is the key to understand many chemical and physical biology problems, for example, surface catalytic reactions and protein folding. Enhanced sampling methods, such as metadynamics and umbrella sampling, can accelerate slow transitions and improve their statistical estimates by biasing molecular dynamics in a low-dimensional reaction coordinate / collective variable (CV) space. The choice of these CVs coordinates determines the efficiency and accuracy of sampling. However, choosing good CVs can be challenging, requiring chemical intuition, experience and many trial-and-error tests. Here, we proposed to use a machine learning dimensionality reduction approach to automate the design of collective variables.
A multi-task machine learning algorithm is employed to learn collective variables from transition path ensembles. In this algorithm, one neural network is trained to map atomic configurations to a lower dimensional latent space, and two additional neural networks are trained to map the latent space to potential energies and metastable state labels. A specific training procedure optimizes the latent space to preserve the information of potential energies and the progression of reaction pathways. The resulting latent space can be used as CVs for metadynamics and umbrella sampling in a wide variety of reactive systems. This approach has been tested on model systems that have two meta-stable states with one reaction pathway, including a Muller-Brown potential model and alanine dipeptide.

We present results from the study of a capacitor made from graphene in contact with water, a simple model of an electrochemical cell. Electrochemistry can be very challenging to simulate at the atomic scale for three reasons: many non-equivalent atoms can be required; the system is at finite temperature with many atoms being mobile; chemical reactions can occur with electrons entering and leaving the simulation cell. From this we can see that the essential ingredients of any successful simulation are: a reliable description of bonding that permits transfer of electrons; an efficient description of bonding suitable for molecular dynamics; a low cost scheme for electron open boundaries. For the description of bonding we use Tight Binding with the Hairy Probe method for open boundaries [1,2]. We see how the water molecules orient themselves under the influence of an applied bias, and how this modifies the capacitance.

[1] “Efficient simulations with electronic open boundaries”, A.P. Horsfield, M. Boleininger, R. D'Agosta, V. Iyer, A. Thong, T. N. Todorov, and C. White, Phys. Rev. B 94 075118 (2016)
[2] “Efficient electron open boundaries for simulating electrochemical cells”, Zauchner, M. G., Horsfield, A. P., and Todorov, T. N., Physical Review B, 97 045116 (2018)

Most of the previous works on the prediction of molecular crystal structures have been based on the (static) lattice energy minimization (LEM), which neglects thermal and zero-point molecular vibrations. This leads to the neglection of thermal effects such as thermal expansion and phase transition. Another more serious drawback of LEM is that it results in too many polymorphs, which can be considered as a kind of artifacts, compared to experiments. Here we propose a new approach for the prediction of molecular crystal structures at finite temperatures by taking into account thermal and zero-point phonons. Our approach, free energy minimization (FEM), minimizes the Gibbs free energy which contains static lattice energy, phonon energy, and entropy-related term. Since the Gibbs free energy can be defined in equilibrium states only, the search space of FEM is not the whole set of the Cartesian product of lattice parameters and atom positions as in LEM but is its subset with constraints that all atoms are at equilibrium positions. In other words, only lattice parameters are regarded as independent variables, and atom positions are dependent variables determined by the constraints. As a proof of concept, we demonstrate the FEM simulation of crystal structures of solid argon. The lattice constants by LEM was found 1.6% smaller than the experimental value at 10 K. On the other hand, the lattice constants by FEM exhibited smaller mismatch, 0.5%, with experiments and also reproduced the thermal expansion quite well. In addition, FEM showed that the minimum point of the Gibbs free energy disappears above 87 K, which is consistent with the experimental boiling point of argon, 87.30 K. The reproduction of the unstability of argon crystal above the boiling point implies that FEM can properly eliminate polymorphs which are unstable at specified temperatures.

Lithium-air batteries are an active area of research because of their potential to have a much higher energy density than traditional lithium-ion batteries. However, they are not yet commercially viable due to poor efficiency, high charging voltages, and low cycle lifetimes. Experimental studies of Li-air batteries with aprotic solvents have shown that the O₂ reduction starts when superoxide (O₂^-) forms in solvent and reacts with Li⁺ to form lithium superoxide (Li⁺-O₂^-). Solid Li₂O₂ then forms as the final discharge product on the cathode. Recent experimental work has suggested that the choice of solvent and the presence of any lithium salts in the system may have a large impact on how the discharge product forms at the cathode. We therefore modelled the clustering of lithium salt molecules in solvent systems without LiO₂ present by performing explicit solvent calculations using both classical molecular dynamics and ab initio molecular dynamics simulations. Properties such as the vibrational density of states, coordination number, and structure factor were calculated for these simulations, and the results from the classical and ab initio simulations were compared both with each other and with experimental data. Large differences in properties such as coordination number were observed for classical versus ab initio simulations, which was not expected based on the widespread use of classical molecular dynamics for modeling solvent behavior. A summary of our results comparing and benchmarking computational approaches for properties such as vibrational density of states and coordination numbers for pure solvent and solvent-salt systems will be presented.

Biomineralization is a widespread process in living organisms through which hierarchical mineral assemblies are formed under moderate conditions of temperature and pH. To this end Silaffin, a naturally occurring protein known to mediate the silica biomineralization process in diatoms, has been widely studied. In particular, R5 - the 19-residue segment of Silaffin, was found to precipitate silica in vitro, and has been extensively employed as a model peptide in many experiments and simulations to elucidate the driving forces of biomineralization. For R5 to facilitate silica or titania precipitation in solution, there must be sufficient negative charge in the system to compensate for the relatively high positive charge on the peptide. In vitro, this has been accomplished with phosphate buffer to mimic post translational modification in the form of serine phosphorylation known to occur in vivo. Sprenger et al. used molecular dynamics (MD) to study the atomic mechanism underlying the phenomenon of protein-mediated biomineralization, focusing on the role of pH and phosphorylation. The interaction of wild type R5, along with two mutants (local phosphorylation and global phosphorylation) to silica was investigated, and it was found that phosphorylation decreased the binding free energy at acidic and neutral pH. While the trends match experimental results of R5 binding on silica, peptide-surface binding is not the primary indicator of biomineralization efficiency. The key role of phosphate suggests peptide-peptide interactions play an important role in the nascent phase of biomineralization and in this work, we study R5 peptide dimerization, to shed light on inter-peptide interactions that enable the uniform templating of silica. To match closely to experimental systems, we study R5 dimers under five distinct conditions. These are (1) wild type R5 (2) R5 with phosphate ions (3) R5 with phosphorylation on Ser1 (4) R5 with phosphorylation on Ser14 (5) globally phosphorylated R5. Enhanced sampling using Parallel Bias Metadynamics with Partitioned Families (PBMetaD-PF) was employed to ensure convergence of the free energy profiles. Results indicate that R5 dimerizes in the presence of phosphorylated serine, or with phosphate ion in solution; and impacts of the number/location of phosphoserine in peptide is also observed. Dimer conformations are analyzed using K-means clustering algorithm. This work provides molecular insight on peptide aggregation in biomineralization, and could potentially inform rational design of protein and nanomaterials.

We present an experimental/simulation integrated study examining crack propagation in nanoporous gold. We utilized a simulation technique that combines continuum fracture mechanics with molecular dynamics. We applied this to a large atomistic digital sample of nanoporous gold to implement mode I crack propagation. Crack propagation tests were also conducted on an experimental sample, prepared via chemical dealloying and observed via in-situ TEM microscopy. We observed cracks in both samples propagating by the same mechanisms of sequential individual ligament failure. A series of nanowire computational deformation tests were also conducted to understand individual ligament behavior, and how this influences the overall sample fracture. This iterative direct experiment/simulation comparison provides new insight into the failure response of nanoporous gold.

Vacancy diffusion under hydrogen (H) environments is important to the mechanical performance of structural materials, such as in hydrogen embrittlement and hydrogen blistering. The vacancy diffusion is believed to be retarded by trapping H-atoms forming H-vacancy complexes. The diffusivity of the H-vacancy complex can be estimated using long-time molecular dynamics simulations or static energetics calculations provided that the thermodynamic equilibrium distribution of H-atoms around the vacancy has been achieved. From the static energetics calculations based on the harmonic transition state theory, we found that the vacancy diffusion coefficient in Cu-H system decreases exponentially as the H-concentration increases by considering the apparent probabilities of various H-vacancy complexes at its ground state and the energy barrier of vacancy diffusion calculated from an embedded-atom-method potential. However, this H-concentration dependency of vacancy diffusivity is denied by the direct molecular dynamics simulations, which shows that H-atoms strongly enhance the diffusion of the vacancy with increasing H-concentration. It should be noted that the number of trapped H-atoms along the diffusion pathway was kept constant in the static energetics calculations. From the viewpoints of statistical sense, the varied number of trapped H-atoms along the diffusion pathway should be allowed. Therefore, we used the potential-of-mean-force method to study the vacancy diffusion under H-environments so that the thermodynamic equilibrium distribution of H around the vacancy can be achieved along the diffusion pathway. We found that the vacancy can trap more number of H-atoms at its saddle-point than that at its ground state, which can be characterized by an activation excess. According to the generalized Gibbs isotherm, the higher chemical potential of H-atoms significantly decreases the energy barrier of the vacancy diffusion. Thus, the vacancy diffusivity is enhanced by H-atoms, which is in consistent with the results of the direct molecular dynamics simulations. This effect of H is believed to be general in other thermal activation processes that can attract more H-atoms at its saddle-point.

WC-Co cemented carbides combine superb hardness with high toughness making them ideal for usage in cutting applications and in wear resistance tools. Thin films of cubic stacking have been observed experimentally at the phase boundary between WC and Co when doping with Ti and V. These films could potentially inhibit grain growth and are thus important to understand as grain size play a crucial role for the mechanical properties of the material. Experimental results have also shown that these types of films are present in undoped WC-Co.

Here we construct, using ab-initio calculations and modeling, an interface phase diagram for thin films in undoped and Ti doped WC-Co . We employ cluster-expansions to model mixed metal layers and carbon vacancies. Monte Carlo simulations are used to sample the configurational entropy. We include vibrational contribution to the free energy by using force-constant regression in order to extract the harmonic free energy. Configurational and vibrational entropy are all shown to be important factors in order to understand this stabilization of cubic thin films in WC-Co. The calculated thermodynamic properties of the cubic thin films are compared to experimental studies and show good agreement.

Alloy cluster expansions (CEs) provide an accurate and computationally efficient mapping of the potential energy surface of multi-component systems that enables comprehensive sampling of the many-dimensional configuration space.
Here, the integrated cluster expansion toolkit (icet), a flexible, extensible, and computationally efficient software package, is introduced for the construction and sampling of CEs.
icet is largely written in Python for easy integration in comprehensive workflows, including first-principles calculations for the generation of reference data and machine learning libraries for training and validation.
The package enables training using a variety of linear regression algorithms with and without regularization, Bayesian regression, feature selection, and cross-validation.
It also provides complementary functionality for structure enumeration and mapping as well as data management and analysis.
Potential applications are illustrated by several examples, including (1) studying chemical ordering and associated properties in a series of intermetallic clathrates as a function of composition and temperature and (2) by predicting the phase diagrams of bulk and surface alloys.

[1] M. Ångqvist, W. A. Muñoz, J. M. Rahm, E. Fransson, C. Durniak, P. Rozyczko, T. H. Rod, and P. Erhart, Adv. Theory Simul. 1900015 (2019)

Density functional theory is the workhorse of theoretical materials investigations. Due to the
shortcoming of (semi-)local exchange correlation potentials, hybrid functionals have been established for practical calculations to describe surfaces, molecular adsorption, and defects. These functionals operate by mixing between semi-local and Hartree-Fock exchange semi-emprically. However, their parameters have to be optimized for every material separately. To treat materials with a more physics driven approach and without the need of parameter optimization is possible with many-body approaches like GW, but at an immense increase in computational costs and without the access to total energies and hence geometry optimization.

We propose a novel exchange correlation potential for semiconductor materials, that is based on
physical properties of the underlying microscopic screening. We demonstrate that it reprocuduces
the low temperature band gap of several materials. Moreover, on the example of defects in semicon- ductors, it respects the required linearity condition of the total energy with the fractional occupation number, as expressed by the generalized Koopman’s theorem. It is shown, that all technologically relevant Gallium-based semiconductors can be treated with a common choice of the functional. We also can treat the imfamous case of ZnO with the same functional. As the parameters of the calcu- lation can in principle be determined from ab initio calculations, our approach can be seen as a non-emprical approximation.

Presently available exchange-correlation functionals in density functional theory (DFT) suffer from a number of errors, including many electron self-interaction error or delocalization error. Eliminating delocalization error (i.e., deviation from piecewise linearity) has been shown to improve property estimation, including relative energies and especially properties directly tied to the improved orbital energies such as the band gap. However, all common approaches to mediating this delocalization error (e.g., global or range separated hybrid functionals or DFT+U) do so at the cost of increasing fractional spin error i.e. the static correlation error, which worsens other essential properties (e.g., bond dissociation energies and activation energies), thereby presenting a challenge for studying correlated materials where both these errors are simultaneously prevalent. We develop a general approach to avoid this trade-off by returning to the exact flat-plane constraint that arises from the union of the fractional spin and piecewise linearity constraints. We examine the shapes of semi-local DFT errors with respect to the flat plane condition and use these deviations to construct new, approximate, functional forms that can recover the flat-plane condition. Our judiciously-modified DFT (jmDFT) approach for constructing few-parameter, low-order corrections to conventional xc-functionals adds no computational overhead to semi-local DFT while recovering the flat-plane constraint. We describe how our approach relates to commonly used Hubbard corrections (i.e., DFT+U and DFT+U+J), thereby emphasizing their relation to improving DFT+U. Using this observation, we outline expressions that can be used to extract jmDFT parameters in a first-principles manner and study how such expressions compare with those employed in the ab initio computation of Hubbard parameters. Finally, we show how to employ this low-cost, non-empirical framework to the study of transition-metal based correlated materials.

The fast-growing field of regenerative medicine is closely related to the advances in tissue engineering and material sciences. Biomaterials have a prominent role as an architectural framework, similar to the innate Extracellular Matrix (ECM), which can re-establish normal functions of deteriorated tissues. In this research, computational methods were used to study the relationship in both molecular and cellular levels of mesenchymal stem cells (MSCs) differentiated into neuronal lineage in electrospun polycaprolactone-graphene (PCL-GRA) scaffold. Specifically, we used molecular dynamics (MD) simulation to investigate the interactions of the scaffold with actin protein and Cellular Potts Model (CPM) to study cellular behaviors such as cell attachment, differentiation, and proliferation. For a better understanding of neuron cell mobility on PCL-GRA scaffold, we investigated the interactions of F-actin protein, which controls the mobility and influences the dendritic spines for neuronal cells. For computational simplicity, we used the actin monomer; globular actin (G-actin). For the entire duration of 100 ns MD simulation, the Van der Waals (VDW) energy between F-actin and PCL-GRA scaffold shows negative values, which represent strong attractive forces between our scaffold and G-actin. We considered G-actin residues within 5 Å from PCL-GRA surface as “interacted residues.” Hydrogen bonds established between LYS 326 residue of G-actin and polycaprolactone chains, which also suggest non-bonding electrostatic interaction between two molecules. For the cellular level interaction, we are using CPM, a computational lattice-based model which uses the stochastic Monte Carlo’s probabilistic model by minimizing the system’s total energy. MSCs extend filopodia and spread along the strain orientation, where PCL fibers are stiffest. The mobility and adhesion of MSCs are controlled by cell-to-cell forces and cell-to-scaffold forces which are incorporated in our CPM model. Our CPM simulation results show that in randomly aligned PCL fibers, cells tend to travel up gradients toward the areas where the graphene concentration is higher. Addition of graphene to PCL scaffold improves cell attachment and viability as compared to only PCL scaffold. Our results indicate that besides the intrinsic factors assigned to the stem cells, the extrinsic parameters of the PCL-GRA interaction have a significant role in the creation of focal adhesion points, cell motility, and cellular attachment.

Catalytic chemical vapor deposition (CCVD) is a widely used method to synthesize carbon nanotubes (CNTs). To control the property (e.g. length, diameter and chirality) of CNT, it is essential to understand the formation mechanism of CNT. Molecular dynamics (MD) simulation is a useful method to observe the formation process. However, MD simulation can only reach the time scale of nanoseconds and is not enough to directly observe the growth of CNTs at low temperature where catalysts remain solid. Therefore, we adopt a recently proposed acceleration method called collective variable-driven hyperdynamics (CVHD) method [1]. In CVHD method, firstly a collective variable (CV) is set and the CV is biased gradually in the same way as metadynamics [2]. Unlike metadynamics, CVHD reset the CV and the bias on it once the system undergo a transition. In this way, CVHD method can reproduce correct dynamics of the system. However, CVHD has a problem that the magnitude of acceleration decreases as the number of target reaction increases. This problem will be fatal in the simulation of CNT because CNT formation process contains migrations of multiple carbon atoms. In this study, we solve this problem by executing CVHD method in parallel (parallel CVHD hereinafter). Specifically, multiple CVs are set according to the number of carbons present in the system, and bias potentials are gradually added to each CVs. When a carbon atoms undergo a transition, the bias potential of the corresponding carbon is reset based on the CVHD method.
In order to confirm the correctness of parallel CVHD, activation energy of carbon atom migration is compared with not accelerated MD. The activation energies are derived from the Arrhenius plot of migration frequency of two carbon atoms in bulk iron. The activation energy gained with parallel CVHD was consistent with not accelerated MD. Subsequently, in order to observe the migration and segregation process of carbon atoms, which is important in the formation of CNT, 20 carbon atoms are randomly put in a slab iron of 1100 atoms and parallel CVHD is performed at 600K. In the 1 ns calculations of not accelerated MD, almost no migration of carbon was observed. Whereas in the parallel CVHD, an acceleration rate of about 700 times is realized, and it was possible to observe atom migration to the surface.

[1] K.M. Bal, E.C. Neyts, J. Chem. Theory Comput. 11 (2015) 4545.
[2] A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. 99 (2002) 12562.

We present a Langevin based stochastic model with spatial correlations and local
environment dependence capable of capturing both the electron-phonon interaction
as well as electronic stopping in classical molecular dynamics (MD). The model is
parametrised by fitting to data from smaller scale first principles simulations, such
as real-time time-dependent density functional theory (TDDFT), resulting with a
solution that does not have tunable additional parameters. The framework is readily
extendable to multicomponent systems and is implemented as an extension to widely
used LAMMPS code. We show the application of our model in the case of Ni crystal,
where we illustrate the non-equilibrium dynamics of phonon modes and a NiFeCr
solid solution. In the latter case the electronic stopping regime as well as thermal
equilibration is illustrated by collision cascade MD simulations at different initial
conditions (PKA energies). Our model is able to excite and damp phonon modes
with different wavevectors with rates in good agreement with quantum theory, which
has been lacking so far in classical simulations studying non-equilibrium dynamics.
Also, in collision damage simulations the model transitions smoothly from electronic
stopping mode, where the coupling is strong, to the electron-phonon regime, weak
coupling, as the projectile slows down and does not explore high electronic-density
regions. In summary, the developed model is able to describe the non-equilibrium as
well as equilibrium dynamics of a system due to electron-ion interaction and could be
further applied to various phenomena such as, laser-matter interaction, compression
shockwave dynamics, and radiation damage.
This work performed under the auspices of the U.S. Department of Energy
by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344
and funded by the Energy Dissipation to Defect Evolution Center, EDDE, an En-
ergy Frontier Research Center of the U.S. Department of Energy (Award Number
2014ORNL1026).

Kinetic Monte Carlo techniques are extensively used to simulate the evolution of defects in materials and to study irradiation effects in particular. All these methods are based on a master equation involving a transition rate matrix. Transition rates are derived from an energetic model in which defects and atoms may either occupy substitute and interstitial sites of a crystalline lattice or exhibit off-lattice in a continuous space. Unfortunately, Kinetic Monte Carlo methods become inefficient when the transition rate matrix associated with the master equation exhibits a broad spectrum of frequencies. Vacancies perform a huge number of transitions between a few atomic configurations connected to each other by small energy barriers. These connected configurations form trapping basins from which the typical escape time for a particle is much higher than the typical time to cross the small barriers inside the basins. The system remains stuck in metastable thermodynamic states. Acceleration techniques based on the theory of absorbing Markov chains are currently being used in Kinetic Monte Carlo simulations to overcome kinetic trapping in low energy basins. This is achieved by drawing escaping events to distant locations from the exact first-passage and no-passage distributions. These two distributions can be formally expressed through the eigenvalue decomposition of the transition rate matrix. Assuming reversible diffusion processes, a property usually satisfied by defects in metals and alloys in and out of equilibrium, we show that the involved linear and eigenvalue problems to be solved can be symmetrized. Furthermore, the slowest modes, those associated with the smallest eigenvalues, are observed to contribute the most to the first-passage and no-passage distributions. As a result, Krylov subspace projection techniques implementing reverse iterations provide us with an efficient tool for accurately approximating the two desired distributions at a relatively low cost. We discuss the convergence, scalability and range of applicability of the approach. We demonstrate its efficiency by computing sink strengths and transition currents for the emission and absorption of vacancies from and to cavities in Aluminium at low defect concentrations.

We present a molecular-dynamics (MD) informed phase field simulation approach to study the microstructure evolution during solidification of Ni-Ti alloys. In the conventional phase-field model developed by Echebarria et al. (EFKP) [1] the kinetic coefficient for solidification of dilute binary alloys is assumed to be zero. However, in the newly-developed phase-field model by Pinomaa and Provatas [2], based on the Crystal Growth Model (CGM), it has been shown that consideration of the kinetic coefficient may play an important role for capturing the effect of solute trapping. Using MD simulations, we calculate the anisotropic kinetic coefficient and crystal-melt interfacial free energy by employing the free solidification and capillary fluctuation method, respectively. Then, using the MD-calculated crystal-melt interfacial properties, we parameterize both the EFKP and Pinomaa- Provatas phase field models and investigate how the partition coefficient changes with the solidification velocity. Work supported by the U.S. National Science Foundation (NSF OIA-1541079).

[1] B. Echebarria, R. Folch, A. Karma, M. Plapp, Quantitative phase-field model of alloy solidification, Physical Review E, 70 (2004) 061604.
[2] T. Pinomaa, N. Provatas, Quantitative phase field modeling of solute trapping and continuous growth kinetics in quasi-rapid solidification, Acta Materialia, 168 (2019) 167-177

The interactions between solute atoms and crystalline defects such as vacancies, dislocations, and grain boundaries are essential in determining physical, chemical and mechanical properties of alloys. Here we present a general correlation between two electronic factors and the solute-defect interaction energies in binary alloys of body-centered-cubic (bcc) refractory metals (such as W and Ta) with transition-metal substitutional solutes. One electronic factor is the bimodality of the d-orbital local density of states for a matrix atom at the substitutional site, and the other is related to the hybridization strength between the valance sp- and d-bands for the same matrix atom. Remarkably, the correlation is independent of the types of defects and the locations of substitutional sites, following a linear relation for a particular pair of solute-matrix elements. The correlation plus a residual-correction model can provide a novel and quantitative guidance to predict the solute-defect interactions in alloys based on electronic structures.

Reverse osmosis membranes are increasingly being used for water desalination and purification. These membranes consist of a very thin (100-200 nm thick) selective layer, often made of crosslinked aromatic polyamide (PA), on top of a porous support. Understanding the mechanism of diffusion across the selective layer is key to improve the performance of current membranes. So far, classical molecular dynamics has been very useful to explain the hydration, water diffusion, and salt rejection mechanisms in plain PA membranes. Recently, nanocomposite membranes containing carbon nanotubes (CNTs) mixed with PA have shown high permeation and salt rejection and good chlorine resistance. However, the mechanism of diffusion is still not completely understood and in fact, there is no agreement on the how CNTs increase the permeability of these membranes. While some researchers suggest that water flows within the embedded nanotubes, other researchers suggest they are impermeable. Unfortunately, the characteristic dimensions of CNTs (diameter around 10-12 nm, and several hundred nm length) requires using large simulation boxes that forbids the use of molecular dynamics. In this work, we used a multiscale approach by combining classical molecular dynamics and a simplified kinetic Monte Carlo approach to study the water diffusion across CNT-PA nanocomposite membranes. Briefly, we built a relatively small cell of PA and PA containing a single CNT. By using classical molecular dynamics, we calculated the free-energy-potential field in both cells. We implemented later a simplified kinetic Monte Carlo algorithm to simulate the diffusion of water molecules across 200 nm thick CNT-PA nanocomposite membranes and we were able to study the effect of CNT's concentrations on water permeation across the membranes.

Deterministic simulations of the rate equations governing cluster dynamics in materials are limited by the number of equations to integrate. Stochastic simulations are limited by the high frequency of certain events. We propose a coupling method combining deterministic and stochastic approaches. It allows handling different time scale phenomena for cluster dynamics. This method, based on a splitting of the dynamics, is generic and we highlight two different hybrid deterministic/stochastic methods. These coupling schemes are highly parallelizable and specifically designed to treat large size cluster problems. The proof of concept is made on a simple model of vacancy clustering under thermal ageing.

The efficient extraction of force constants (FCs) is crucial for the analysis of many thermodynamic materials properties. Approaches based on the systematic enumeration of finite differences scale poorly with system size and can rarely extend beyond third order when input data is obtained from first-principles calculations. Methods based on parameter fitting in the spirit of interatomic potentials, on the other hand, can extract FC parameters from semi-random configurations of high information density and advanced regularized regression methods can recover physical solutions from a limited amount of data. Here, the hiphive Python package, that enables the construction of force constant models up to arbitrary order is presented. hiphive exploits crystal symmetries to reduce the number of free parameters and then employs advanced machine learning algorithms to extract the force constants. Depending on the problem at hand, both over and underdetermined systems are handled efficiently. The FCs can be subsequently analyzed directly and or be used to carry out, for example, molecular dynamics simulations. The utility of this approach is demonstrated via several examples including ideal and defective monolayers of MoS2 as well as bulk nickel.
[1] Advanced Theory and Simulations 2019, 1800184

Optical characterization of carbonaceous materials (e.g., Laser Induced Incandescence, LII) is important in fire detection and manufacturing of various grades of carbon black and even their impact on climate largely depends on light absorption. For example, the contribution of carbonaceous materials to global warming (~25 % of the total anthropogenic contribution) has currently the highest uncertainty among all contributions (CO₂, CH₄, etc.), as their optical properties depend on their precise composition quantified by C/H and morphology. Understanding the light absorption of carbonaceous particles is essential to their characterization and selective sensing by fire detectors and can improve the accuracy of climate models.
Here, the impact of composition on the light absorption from carbonaceous particles is investigated by coupling mesoscale Discrete Element Modeling (DEM) [1] with Discrete Dipole Approximation (DDA) during their surface growth and agglomeration [2]. The Mass Absorption Cross-section, MAC, of these agglomerates is estimated by DDA and validated against atomistic point dipole interactions and mesoscale DDA simulations [2]. Using a refractive index, RI, for mature soot with large C/H > 4.5 [3] yields constant MAC and absorption function E overestimating the light absorption up to 75 % of nascent soot with C/H < 4.5. The RI is interpolated between those of nascent and mature soot for wavelength, λ = 1064 nm to account for quantum confinement and evolving number of clustered sp²-bonded rings. This affects the optical band gap, E_g, resulting in excellent agreement of the DEM-derived evolutions of MAC and E with the corresponding LII measurements in methane [4] and ethylene premixed flames [5]. The E_g of nascent soot decreases during agglomeration, increasing MAC and E by 65 %. The good agreement between DEM and LII data confirms that carbonaceous particle dynamics by surface growth and agglomeration strongly correlate with particle composition and RI that are essential for quantifying carbonaceous particle light absorption.

References:

[1] Kelesidis, G. A., Goudeli, E., & Pratsinis, S. E. (2017). Proc. Combust. Inst., 36, 29-50.
[2] Kelesidis, G. A., & Pratsinis, S. E. (2019). Proc. Combust. Inst., 37, 1177-1184.
[3] Yon, J., Bescond, A., & Liu, F. (2015). J. Quant. Spectrosc. Radiat. Transf., 162, 197-206.
[4] Bejaoui, S., Batut, S., Therssen, E., Lamoureux, N., Desgroux, P., & Liu, F.S. (2015). Appl. Phys. B, 118, 449-469.
[5] Olofsson, N.E., Simonsson, J., Torok, S., Bladh, H., & Bengtsson, P.E. (2015). Appl. Phys. B, 119, 669-683.

Mesoporous materials have attracted considerable attention in the past for their applications in catalysis, sorption, separation, hydrogen storage, oil & gas extraction, etc. and as host matrices for optically and electronically active components. In the recent past, mesoporous siliceous materials (MSM’s) have also shown tremendous potential in drug delivery technologies. Their superior thermophysical, mechanical, and chemical properties such as variable pore size, pore morphology, and surface functionalizability create avenues for optimized interactions between drug molecules and carriers. A widely adopted method of synthesis of these structures is evaporation-induced self-assembly (EISA). In EISA, control of pore morphology and liquid crystalline mesophases is accomplished by controlling the environmental conditions such as temperature, orientation of micelles, pH, and also the solvent evaporation rate and the block copolymer used as precursor. An important concern is that the packing of individual mesoporous solids in to a container can result in internal pores between individual mesoporous structures of larger sizes (>50 nm). Hence a monolithic solid is desired and has to be synthesized directly from the precursors. Synthesis of monoliths is highly challenging as it results in the formation of cracks. Prior studies have shown that cracks can be partially alleviated by controlling solvent evaporation rates, supercritical CO₂ drying, or by using liquid paraffin in the synthesis mixture.
While the usual practice is to experiment with the block copolymers, surfactants, rate of evaporation, and the time involved in polymerization, proper guidelines for the rational synthesis of mesoporous siliceous monoliths is not available yet. Hence, as the first stride to address this issue, here, we present a molecular dynamics (MD) simulation routine that can serve as a guideline to the synthesis of mesoporous structures. Combining non-reactive and reactive MD, we simulate the non-reactive self-assembly process followed by chemical events leading to the evolution of porous structures. We implement a 3-step process using a water-solvated pluronic-based tri-block copolymer having both hydrophilic and hydrophobic groups in it. In the first step, we equilibrate the system with concentration > critical micelle concentration and temperature higher than the critical micelle temperature allowing for the non-reactive self-assembly of the surfactant molecules. OPLSAA force field [1] is used to drive this co-operative assembly. Secondly, we introduce silicic acid as the silica precursor for the micelle system and equilibrate using ReaxFF reactive force field [2] to enable chemical reactions leading to the formation of skeletal Si-O-Si bonds. To condense the silica precursor to the organic template, we use the accelerated bond-boosting scheme [3] within the ReaxFF molecular dynamics framework, which essentially bypasses the transition state complexes with high activation energy requirements, thereby achieving laboratory timescales in MD simulations. Following this, as the final step, by preferentially evaporating water solvent, MSM structures are obtained. We study the structural and mechanical properties of MSM’s to gain insights on molecular structure, localized heterogeneities, crack formation, etc., thereby elucidating the process of synthesis of mesoporous siliceous materials.

References
[1] Jorgensen et al., J. Am. Chem. Soc., 1996, 118, 11225-11236
[2] van Duin et al. JPC-A 2001, 105, 9396
[3] A. Vashisth et al., JPCA 2018, 122, 6633-6642.

Hydrogen is a promising clean fuel. A key bottleneck to its widespread utilization is the susceptibility of materials used in hydrogen distribution infrastructure (typically steels) to embrittlement and failure. This necessitates the use of permeation barriers to prevent hydrogen ingress into the steel. Al₂O₃ is a popular hydrogen barrier, due to its low hydrogen permeability. Our recent study on bulk Al₂O₃[1] using an ab initio thermodynamics approach has shown that the hydrogen solubility and diffusivity in Al₂O₃ can be further reduced by dopant engineering. Specifically, doping with Si, Ti, Fe and Cr reduces the solubility of free protons to negligible amounts, and traps protons at aluminum vacancy sites with high binding energy (~3 eV), in turn, reducing diffusivity.
While the intrinsic hydrogen permeation resistance of bulk Al₂O₃ is attractive, it is more lucrative to use a multilayer structure of alternating Al₂O₃ and Al layers as the coating. Such a structure is advantageous on multiple fronts: (a) it has high toughness and fatigue resistance, preventing cracking of Al₂O₃, thus avoiding fast diffusion pathways for hydrogen ingress to the steel, (b) it is self-healing, and (c) space charge engineering of the high interfacial area can help reduce hydrogen permeation further. The latter is the aim of this study.
Density functional theory was used to obtain the formation energies of various native, hydrogen and dopant (Mg, Si, Ti, Fe, Cr, B, C, N) defects in the bulk and defect segregation energies to the Al₂O₃|Al interface core. Using these energies as input, a continuum model was developed, which self-consistently solved Poisson’s equation to obtain the defect concentrations at the interface. Results show that in the undoped case, free hydrogen interstitials are depleted in the space charge region and accumulate at the core. This hydrogen buildup can lead to blistering of the interface, thus, deteriorating coating performance. Dopants were then utilized to prevent hydrogen accumulation at the core and ensure trapping of hydrogen in the space-charge region. To identify the impact of the interface on hydrogen diffusivity, effect of dopants on hydrogen migration barriers across the interface was studied using nudged elastic band method. The work of adhesion of doped interfaces was also calculated, which has implications for the effect of dopants on the spall tendency of the interface.
The robust coupling of ab initio methods and continuum modeling in this study helps take advantage of the accuracy of first principles calculations while overcoming the spatial limitation of such methods. It is emphasized that the fundamental results of this study can be extended to other systems where the metal/oxide interface affects device efficiency, such as the electrode/electrolyte interface in resistive switching devices and the metal/dielectric interface in metal-oxide-semiconductor devices.

[1] V. Somjit, B. Yildiz, Doping α-Al2O3 to reduce its hydrogen permeability: Thermodynamic assessment of hydrogen defects and solubility from first principles, Acta materialia 169 (2019) 172-183.

Chiral (Dzyaloshinskii-Moriya) magnetic interactions give rise to a wide range of complex magnetic phases, from helical and cycloidal spin textures to topological phases such as skyrmions. Due to the sensitivity of these interactions to precise atomistic environments and their symmetry, even mild strain can profoundly change the character and relative stability of the resulting magnetic phases. We first use phenomenological arguments to chart out how strain may alter chiral magnetic interactions, identifying the strain fields necessary to stabilize various types of magnetic phases in materials of any symmetry. In order to test these predictions, we construct strain-coupled magnetic phase diagrams from first-principles, by relying on density functional theory, cluster expansions, Hamiltonian Monte Carlo, and phase field techniques. This scale-bridging model reveals that strain fields readily accessible by a wide range of experimental techniques can be used to precisely control magnetic phase behavior and the emergence of topological phases in common high-symmetry ferromagnetic materials.

Electronic, optical and electrical properties of the point defects in cubic unit cell ZrO₂ were calculated by using the density functional theory (DFT) calculation based on full potential linear augmented plane wave (FP-LAPW). The exchange-correlation potential is treated by the Gradient Generalization Approximation (GGA) and the recently proposed by Tran and Blaha the modified Becke-Johnson potential (TB-mBJ) which successfully corrects the band-gap problem found with GGA for a wide range of materials. The energetics of point defects is studied after optimization of positions. Our results indicate that the Zr_O antisite has the lowest formation energy and could thus be the acceptor defect responsible for the p-type conductivity of undoped cubic ZrO₂. However, based on our electronic structure obtained using the TB-mBJ method we have calculated various optical properties, including the complex index of refraction, absorption and transmittance coefficients. Furthermore, the electrical properties of each point defects have been investigated using the Boltzmann transport theory.

Neuromorphic computing paradigm holds much promise for applications that require big data analysis and artificial intelligence. The key components in neuromorphic circuits are devices that exhibit a history-dependent conductivity modulation, which could efficiently represent synaptic weights in artificial neural networks. However, an accurate and fast modulation of the conductance over a wide dynamic range is proved to be challenging. Recent experimental studies showed that a gradual change in electronic conductivity of WO₃ can be achieved by hydrogen doping, giving rise to the ability to have a large number of conduction states. However, the relationship between hydrogen concentration and electric conductance is not linear and WO₃ could undergo phase transitions upon hydrogenation. In this work, atomistic density functional theory calculations are used to establish the correlation between H content and H_xWO₃ conductivity and determine the rate-limiting steps of hydrogen intercalation process. This allows us to determine the most relevant hydrogen concentration range and guide experimental studies.

Electronic transport has been traditionally studied by means of the reciprocal lattice method, which has successfully explained many physical properties of crystalline solids through the electronic band theory, such as the optical absorption and the thirty-orders-of-magnitude difference in the electrical resistivity between a good conductor and a typical insulator [1]. However, current electronic devices frequently contain multiple non-periodic structural interfaces and their presence avoids the use of reciprocal space. An alternative way to investigate the electron and phonon transport at atomic scale in macroscopic thermoelectric devices, which allow a direct conversion between thermal and electrical energies without harmful pollution, could be by using the real space renormalization method [2]. This method combined with the convolution theorem has the advantage of being numerically accurate and computationally efficient, i.e., able to address 10²⁴ non-periodically located atoms without introducing additional approximations [3]. It is well known that the metals are poor thermoelectric materials. In contrast, segmented metallic nanowires have a band structure by design, which combines with a properly placed chemical potential by applying a gate voltage could archive a good thermoelectric power [4,5]. In this work, we study thermoelectric properties of periodic and quasiperiodically segmented macroscopic-length nanowires with and without branches by means of the Kubo-Greenwood formula, in which tight-binding and Born models are respectively used for the calculation of electric and lattice thermal conductivities [6-8]. The results show a substantial growth of the thermoelectric figure of merit (ZT) induced by the long-range quasiperiodic disorder, because it diminishes the thermal conduction of long wavelength acoustic phonons being such phonons usually not altered by local defects neither impurities [9].

This work has been partially supported by UNAM-IN106317, UNAM-IN115519 and CONACyT-252943. Computations were performed at Miztli of DGTIC, UNAM.

[1] C. Kittel, Introduction to Solid State Physics, 8^th Ed. (Wiley, New York, 2005).
[2] V. Sanchez, L. A. Perez, R. Oviedo-Roa and C. Wang, Phys. Rev. B 64, 174205 (2001).
[3] V. Sanchez and C. Wang, Phys. Rev. B 70, 144207 (2004).
[4] Y. Tian, M.R. Sakr, J.M. Kinder, D. Liang, M.J. MacDonald, R.L.J. Qiu, H.-J. Gao, and X.P.A. Gao, Nano Lett. 12, 6492 (2012).
[5] S.C. Andrews, M.A. Fardy, M.C. Moore, S. Aloni, M. Zhang, V. Radmilovicbd and P. Yang, Chem. Sci. 2, 706 (2011).
[6] C. Wang, F. Salazar and V. Sanchez, Nano Lett. 8, 4205 (2008).
[7] J. E. Gonzalez, V. Sanchez and C. Wang, J. Electron. Mater. 46, 2724 (2017).
[8] J. E. Gonzalez, V. Sanchez and C. Wang, MRS Communications 8, 248 (2018).
[9] F. Sanchez, C. Amador-Bedolla, V. Sanchez and C. Wang, J. Electron. Mater. (2019) doi: 10.1007/s11664-019-07298-0.

Spontaneous self-organization of solid binding peptides on single layer atomic materials (SAP/SLAM) offers a great potential in exploiting these systems for bioelectronic and biosensing applications. The molecular interaction at the bio/nano soft interface between the peptide and the solid material is dependent on the conformations of the peptides and the surrounding conditions such as pH and temperature that affect the molecular structures, conformations, binding, and self-assembly in a myriad of ways. Although molecular dynamics has been used to study protein folding problems, here we demonstrate a closed-loop data pipeline molecular-dynamics implementation in SAP/SLAM systems. We perform molecular dynamics simulation on a wild type peptide, obtained from directed evolution experiments, that self-assembles on graphene at neutral pH and room temperature. Several mutants of the WT peptide are used to investigate whether they display spontaneous self-assembly under different pH and temperatures. Principal components analysis is paired with gaussian mixture model clustering to obtain full-peptide conformation propensities and automated analysis is done to obtain information about the conformations at different processing conditions on graphene and in aqueous solution. We show that while the WT peptide assembles on graphene at neutral pH and room temperature, the charge-neutral mutant assembles under a wide range of pH (3.0 < pH <9.0) and has a highly similar conformations over the range of pH conditions. Similarly, the tyrosine to tryptophan mutant assembles at higher temperatures with a clear change in conformations from that of low temperature conditions. Knowing conformations of peptides and how the mutants affect conformation at different conditions is crucial in understanding the molecular footprint of the peptide on the surface. The knowledge of soft bio/nano interface structure will facilitate DFT calculations downstream to understand effect of self-assembly on the band structure changes happening in the SAP/SLAM systems. These studies eventually will lead to the predictive atomistic algorithms design of unique bio/nano systems defined by genetically selected peptides and single atomic layers materials towards biosensors, bioelectronics and biomolecular fuel cells. As part of the Materials Genome Initiative, the research is supported by NSF-DMREF program through the grant DMR-1629071.

Direct methane-to-methanol conversion, which enables transformation of gas-phase methane to liquid methanol under mild conditions, is considered a dream reaction. At present the reaction is often realised via a three-step cycle; an activation phase, a reaction phase, and finally an extraction phase. Copper ion-exchanged zeolites have been suggested as promising catalysts for this reaction, although the nature of the actives sites is still under debate. Here we examine, using ab-initio molecular dynamics and first-principles thermodynamics, oxidation state and coordination of Cu-dimers inside the chabazite structure Cu-SSZ-13 under relevant experimental conditions. At room temperature, the formation energy of the Cu-dimers shows a multitude of Cu₂(H_xO_y) clusters being exergonic. However, at the higher temperatures of the relevant reaction conditions only Cu₂O and Cu₂(OH) remain as thermodynamically stable structures. In the three step cycle, these are thus identified as the relevant structures for the activation and extraction phase, respectively.

Rosette Nanotubes (RNTs) are tubular soft materials self-assembled from guanine-cytosine (GΛC) hybrid molecules, which can be covalently functionalized for use in various applications. They have been shown to be effective assemblies in regenerative medicine, drug display and delivery, and catalysis.

Previous spectroscopic studies suggest they have a structure which is formed by hexameric rings, maintained by self-complementary hydrogen bonds, that are further stacked and supported by π-π interactions, which form the ring-stacked RNTs. While this mode of association maximizes the hydrogen bonding interactions and results in efficient π-π stacking, it is also possible to envision that the GΛC modules assume a helical organization defining a tubular core. Imaging the nanotube using scanning tunneling microscopy (STM) proved inconclusive in differentiating the two configurations. Furthermore, comparing the calculated chemical shifts of the nitrogen atoms of the GΛC monomers of both the stacked-tube (ST) and helical tube (HT) configurations showed no significant difference between them.

We further investigated this possibility by using the lysine-functionalized RNT (K1-RNT) and applying multi-scale molecular modeling methods, which include Monte-Carlo conformational search, molecular dynamics, and the statistical mechanical theory of molecular liquids, 3D-RISM theory. We considered three structures of the K1-RNT: ST, left-handed HT (LHT), and right-handed HT (RHT). Our results suggest that the formation of ST, LHT, and RHT K1-RNTs in water are favorable and are enthalpically driven. Moreover, 3D-RISM analysis suggests that the RHT conformation is more probable than the ST and LHT RNTs and this is due to a more favorable solute-solvent interaction energy. Using the 3D distribution of solvent sites around the RNTs (from 3D-RISM calculations), we were also able to determine the solvation structure and estimate the free energy of binding of each water molecule in the RNT channel. The results also corroborate the higher probability of the RHT conformation occurring. Studies are being done to verify these findings.

The chemistry of graphene is relatively well-known at relatively low temperatures, when the ambient gas is artificially limited to atomic or molecular oxygen in ultra-low pressure conditions for Electron Microscopy studies and epoxy is typically the only adsorbed surface species. The same applies at regular ambient conditions when air consists only of molecular species such as oxygen, nitrogen or water vapor, which have negligible small degrees of dissociation. Chemisorption of these molecular species is effectively prevented by fairly large activation energies, which severely limits their chemisorption concentrations and therefore play no role in the graphene surface chemistry.

However, carbon materials are extensively used in high-temperature conditions such as Thermal Protection Systems for Hypersonic flight applications. Carbon oxidation plays in these applications a critical role because it actually gasifies carbon in the form of gaseous CO and CO₂ therefore removing the protective layer. However, the chemistry of the graphene surface is not well known in these high temperature environments, which are characterized by large concentrations of highly reactive radicals such as hydroxyl or atomic oxygen, nitrogen and hydrogen. These species chemisorb on the graphene surface with very low or even no energy barrier at all, so they can reach high surface concentrations, competing with oxygen for carbon sites and therefore effectively reducing the oxygen coverage. However, the energetics and kinetics of these species on the graphene surface is almost completely unexplored.

Using DFT calculations, we have systematically studied the kinetics and dynamics of surface processes such as chemisorption, desorption and surface hopping on the graphene surface of atomic hydrogen, nitrogen, oxygen and hydroxyl, which typically populate the graphene surface in these high-temperature environments. We have also studied the interaction and reactions among them. Our ultimate goal is to develop a detailed kinetic mechanism including surface species and reactions among them in much the same way as those already available in combustion of hydrocarbon fuels. We have found for instance that oxygen can chemisorb in an on-top-of-carbon configuration, in addition to the well-known epoxy form. We have studied as well interaction effects among near neighbors, showing actually that clustering is energetically favored, an effect already observed in the case of atomic chemisorbed hydrogen.

On the other hand, we have studied the oxidation dynamics. Compared to the dynamics of these surface species, oxidation occurs in times scales several orders of magnitude slower, which make oxidation events very infrequent. In order to overcome this rare-event limitation, we have used kinetic Monte-Carlo simulations to study the gasification dynamics. In particular, we have been able to simulate and explain the well-known pitting, so characteristic of graphene and graphite oxidation. We show in particular how, in the more relevant conditions for applications high-temperature environments, the presence of all those other surface species modifies the oxidation and pit growth rates with respect to the corresponding values found in the highly controlled conditions of Electronic Microscopy studies.

Improved understanding of aqueous solutions at charged graphitic interfaces is critical for designing new carbon-based materials for energy storage and water desalination. However, many mechanistic details of the interfacial interactions remain unclear, including how interfacial structure and response are dictated by intrinsic properties of solvated ions under applied voltage. This knowledge gap is partly related to the challenge in simulating aqueous interfaces in the presence of a bias potential. Here, we combine a newly developed hybrid first-principles/continuum simulations with electrochemical measurements to investigate adsorption of several alkali-metal cations at the interface with charged graphene and within graphene slit-pores. We confirm that adsorption energy increases with ionic radius, while being highly dependent on pore size under confinement conditions. In addition, in contrast with conventional electrochemical models, we find that interfacial charge transfer contributes non-negligibly to this interaction and can be further enhanced by confinement. Overall, the measured interfacial capacitance trends result from a complex interplay between voltage, confinement, and specific ion effects--including ion hydration and degree of charge transfer. The results provide broad implications for optimizing electrochemical interfaces for energy storage and water desalination.

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

Performance of dielectric polymers under high electric field is limited by the electrical breakdown, which is commonly understood as an avalanche of processes such as carrier multiplication and defect generation, triggered by field-accelerated hot electrons and holes. We model the hot-carriers transport in dielectric polymer, Polyethylene, with excited-state quantum molecular dynamics simulations in presence of electric field, which reveal multiple microscopic processes induced by hot electrons and holes under an electric field ranging between 100 and 1200 MV/m. We found that electronic-excitation energy is rapidly dissipated due to strong electron-phonon scattering. This, in turn, creates other electron-hole pairs to cause carrier multiplication. The key chemical damage occurs due to localization of holes of the surface of slabm which weaken carbon-carbon bonds on the surface. We will also discuss the effects of field strength and pre-existing surface defects (e.g. –COOH and –OH) on the field-induced hot-carrier dynamics and chemistry. Such quantitative information can be incorporated into first principles-informed, predictive modeling of dielectric breakdown.

This work was supported by the Office of Naval Research through a Multi-University Research Initiative (MURI) grant N00014-17-1-2656. The computations were carried out at the Center for High Performance Computing of the University of Southern California

The mesoscale connects the molecular and atomic to the bulk. This size regime has posed a conundrum for modeling, as it is outside the realm of density functional theory (DFT) and molecular dynamics (MD) may not be appropriate. We use a chemical graph theory approach and compute the adjacency matrix for the clusters. This leads to identification of the coordination of each site. We demonstrate mesoscale relationships for the magic formulas of clusters, the radial distribution function (RDF) and nanocatalysis. The number of bonds, atoms, and coordination numbers exhibit magic number character versus the number of shells, n, as the size of the clusters increases. Next, we use gold nanoclusters as model systems to generate their corresponding RDFs. The integrated intensity of the first peak in the RDF is related to magic formulas for the number of bonds (neighbors) and the number of atoms, as a function of n. Each unique cluster type and size has a structural fingerprint in its corresponding RDF. Lastly, the coordination style method is key to studying nanocatalysis and creating a mesoscale model free of arbitrary parameters for sizes ~ 3 nm < D < 100 nm. We apply DFT, informational statistical mechanics, and thermodynamics to complete the model. This provides a theoretical scheme to determine the size and shape dependence of the entropy and enthalpy of nanocluster-adsorbate systems. We present examples from the literature, where others model cluster properties such as melting temperatures, cohesive energy, magnetic moments, Curie temperature, and the melting entropy and enthalpy using coordination numbers. We summarize the potential of these techniques to reach any mesoscale size, using magic formulas.

Molecular dynamics (MD) simulations provide a general technique to study thermal vibrations/phonons and their contributions to thermal transport. However, despite the widespread usage of MD to study phonons, direct comparisons between experiments and simulations are often associated with low fidelity, due to a lack of accurate interatomic potentials (IPs). This issue therefore imposes a major barrier to utilizing MD for studying thermal transport and deriving new physically meaningful insights. Towards solving this problem, we present a new approach to making EIPs specifically optimized for accurately describing phonons and thermal vibrations. The approach enables nearly exact reproduction of ab initio phonon dispersion relations (i.e., < 1%) error), along with accurate forces and thermal conductivity (i.e., < 10% error), with low computational expense like that of traditional EIPs.

A significant challenge in the development of functional nanomaterials is understanding the growth of colloidal metal nanocrystals. One bottleneck is understanding the evolution of nanocrystal seeds from nuclei. Seeds play a critical role in nanostructure formation, as the shape of the seed often determines the final nanocrystal shape. The interplay between structural transitions, deposition, and etching as seeds grow from nuclei leads to a myriad of metastable nanocrystal morphologies. To understand this interplay and its ramifications for seed structure, we conducted replica-exchange molecular dynamics simulations to quantify how the structures Ag and Cu seeds evolve in vacuum, as well as in a solution environment, with increasing size and as a function of temperature. Using common-neighbor analysis, we create barcodes to describe the distribution of nanocrystal structures as a function of size and temperature. We obtain energy barriers for the transformations of selected seeds from one shape to another. Several critical branching sizes emerge, at which intervention by solution processing may influence the final nanoocrystal shape. Knowledge from these studies will enable the development of strategic synthesis protocols to achieve selective nanostructure yields.

Phase transformation microstructures have an abrupt change in lattice geometry that significantly affects material properties. For example, a cubic to tetragonal lattice transformation induces polarization in ferroelectric materials, and a similar lattice transformation during charging/discharging of battery materials induces large volume changes. Current material models decribe phase transformation phenomena at different length and time scales. On the one hand, continuum methods such as phase-field models describe microstructural evolution by typically homogenizing the crystallographic texture of the material, and do not model individual lattice transformations. On the other hand, atomistic models such as the first principle methods and molecular dynamics models provide atomistic insights into lattice transformation, but are not applicable to the length scale and time scales for microstructural evolution. In my talk, I will introduce a novel theoretical framework that couples a Cahn-Hilliard model (continuum scale) with a phase-field crystal model (atomistic scale) to describe phase transformation microstructures in crystalline materials. I will show how this modeling technique can be applied to engineer crystallographic features of phase transformation materials (such as battery materials, ferroelectrics, ferromagnets) to enhance their material properties and lifespans.

Nanocrystalline alloys can be designed to be stable by intentional grain boundary segregation: by using alloying elements that will segregate to the boundaries and lower the grain boundary energy, the driving force for grain growth can be attenuated or even eliminated. The main thermodynamic data necessary for predicting stability in nanocrystalline alloys is the segregation strength of the alloying element, quantified by the segregation energy. Segregation energy is usually calculated precisely only for specific high-symmetry boundaries that are not representative of the spectrum of states in a true polycrystal, or estimated in a semi-empirical manner for some “average” grain boundary assumed typical of a polycrystal. In this work, we outline the thermodynamic and computational framework to address equilibrium grain boundary segregation in full polycrystalline grain boundary networks, by developing a variety of atomistic tools. Specifically, we aim to understand the full spectrum of segregation energies, as well as the role of solute-solute interactions at the boundary. Our efforts to apply these methods across many alloy systems will also be addressed.

The interactions between solute atoms and crystalline defects such as dislocations and grain boundaries play an essential role in determining physical, chemical and mechanical properties of solid-solution alloys. In recent years, the ability to predict solute segregation at high symmetry grain boundaries from first principles have been widely studied. However, previous algorithms have mainly focused on the simple grain boundary structures for dilute solute cases due to the costly computation power needed by density functional theory (DFT). Here, we present a general atomistic approach to optimize the structures and simulate solute segregation trends of grain boundaries in multiple component systems by the combination of a highly efficient genetic algorithm and the grand canonical ensemble, in which components are not restricted to dilute or stoichiometric cases. Different chemical potential can be used as input for creating different reservoirs for the grain boundary phases. In our study, thousands-atom grain boundary systems will be investigated by well-established empirical potentials (MEAM or EAM potentials) for Mg-based alloy systems, like Mg-Y or Mg-Zn, which are potential candidates for lightweight structural components as a result of their low density and high specific strength. Because of the complicated potential energy landscapes (PEL) coming from both geometric and occupational freedom, we will either average a good amount of small configurations by Boltzmann statistics based on their energy distributions for patterned segregation systems or use a large supercell to study cluster segregation systems across the grain boundaries. Final structures will then be used to investigate the effect of solute on mechanical behaviors of grain boundary systems.

We will present ultra-large scale Molecular Dynamics simulations of aluminum single crystals subjected to uniaxial tension. We observe that appearance (or not) of three-stage hardening depends on the initial crystallographic orientation of the straining axis and results from crystal rotation. In its turn, crystal rotation is a direct consequence of co-axiality forced on the specimen by the testing machine, the view widely accepted in the phenomenological Crystal Plasticity community since the classical studies of Schmid predating dislocations. Remarkably, stress-strain behaviors, slip system activity and crystal rotations observed in our high-rate MD simulations are in an exact qualitative agreement with quasistatic experiments suggesting that physics of crystal plasticity scales, i.e. dislocation mechanisms defining crystal plasticity response remain the same over twelve decades of straining rates.

When the dynamics of a high dimensional system is characterized by long periods in metastable energy basins with rare interbasin transition, a rigorous coarse grained representation is a Markov chain with suitable transition rates. If well constructed, such models can provide a valuable, robust connection between atomistic and mesoscale material models. However, with finite resources sampling is always incomplete, which if uncontrolled can have catastrophic consequences on the model validity. Furthermore, commonly used harmonic rate theories can fail at surprisingly low homologous temperatures. I will talk about recent work that uses a robust Poisson-Bayes estimator of sampling completeness to autonomously manage and optimize massively parallel sampling of defective crystal energy landscapes[1] and a novel path-based scheme to evaluate anharmonic contributions to transition state theory predictions of reaction rates[2].

[1] TD Swinburne and D Perez, Self-optimized construction of transition rate matrices from accelerated atomistic simulations with Bayesian uncertainty quantification, Physical Review Materials 2018
[2] TD Swinburne and MC Marinica, Unsupervised calculation of free energy barriers in large crystalline systems, Physical Review Letters 2018

A variety of experiments on submonolayer metal heteroepitaxial growth have demonstrated that the combined effects of strain and bonding can lead to interesting shape transitions and/or defect formation, although the corresponding mechanisms and/or defect structure are not well understood. In order to better understand the underlying mechanisms as well as separate out the effects of strain from “chemical” effects, we have carried out temperature-accelerated dynamics simulations [1] of a simple model system corresponding to the submonolayer growth of Cu on a biaxially strained Cu(100) substrate at 200 K. In the case of 2% compressive strain we find - as was previously found for Cu/Ni(100) growth [2] – that multi-atom “pop-out” events which create open step-edges are enhanced, leading to the formation of blobby islands with a mixture of open and closed step-edges. In contrast, for larger (4%) compressive strain, the presence of submonolayer islands is found to promote the formation of substrate vacancies which eventually (once they come close to each other near an island) leads to the formation of stacking faults in both the substrate and island. The kinetics of this transition is characterized by a large activation barrier along with a prefactor which is tens of orders of magnitude larger than is typical of atomic processes in fcc metals. A theoretical estimate for the prefactor based on an estimate of the change in entropy in going from the initial state to the transition state is found to be in good agreement with our simulation results. In contrast, the prefactor for the reverse transition is much smaller, leading to a negligible rate at 200 K. Finally, we consider the case of 8% tensile strain for which the dominant monomer diffusion mechanism corresponds to exchange with the underlying substrate. In this case a shape transition is observed with increasing island size from isotropic to strongly anisotropic islands. A theoretical estimate of the critical size for the shape transition based on recent theoretical results [3] will be presented which is in good agreement with our simulations.
[1] M.R. Sorensen and A.F. Voter, J. Chem. Phys. 112, 9599 (2000).
[2] Y. Shim and J.G. Amar, Phys. Rev. Lett. 108, 076102 (2012).
[3] J.G. Amar, Y. Shim, and R.T. Deck, Surf. Sci. 616, 120 (2013).

We present a simple and accurate computational technique to determine the frequency prefactor in harmonic transition state theory without necessitating full phonon density of states (DOS) calculations. The atoms in the system are partitioned into an “active region,” where the kinetic process takes place, and an “environment” surrounding the active region. It is shown that the prefactor can be obtained by a partial phonon DOS calculation of the active region with a simple correction term accounting for the environment, under reasonable assumptions regarding atomic interactions. This model enables the use of plane wave density function theory calculations for diffusion description in large systems where the full Hessian calculation becomes intractable. Convergence with respect to the size of the active region is investigated for different systems and different diffusion processes, as well as the reduction in computational costs when compared to full phonon DOS calculation. Additionally, we introduce an open source implementation of the algorithm that can be added as an extension to Large-scale Atomic/Molecular Massively Parallel Simulator software.

Extended Lagrangian Born-Oppenheimer molecular dynamics enables stable simulations without relying on an expensive self-consistent field optimization prior to the force evaluations. The formulation is different from extended Lagrangian Car-Parrinello molecular dynamics and samples the potential energy surface with an accuracy that is of 4th-order in the integration time step instead of 2nd-order. No fictitious mass parameters need to be tuned and the integration time steps is limited only by the integration time step of regular Born-Oppenheimer molecular dynamics. However, the formalism includes a particular metric kernel that needs to be approximated. For most systems a simple scaled delta function is sufficient, but for more challenging materials a more accurate approximation is needed. I will present some recently developed low-rank kernel approximations using Krylov subspace expansions in combination with preconditioners that allow stable (SCF-free) simulations even of highly challenging systems including chemical reactions.

This talk will present the results of a recent phase field crystal modelling study of electromigration. The newly formulated model successfully captures some well-known electromigration induced phenomenologies like the Blech effect and mean time to failure, as well as the failure mechanisms in interconnects at atomic scale through formation of voids at triple junctions and surface steps motion. We also discuss extensions and future directions of coupling this approach with first principles calculations.

Rough energy landscapes appear in a variety of applications including disordered media and soft matter. In this work, we examine challenges to sampling from Boltzmann distributions associated with rough energy landscapes. Here, the roughness will correspond to highly oscillatory, but bounded, perturbations of a fundamentally smooth landscape. Through a combination of numerical experiments and asymptotic analysis, we demonstrate that the performance of Metropolis Adjusted Langevin Algorithm can be severely attenuated as the roughness increases. In contrast, we prove, rigorously, that Random Walk Metropolis is insensitive to such roughness. We also formulate two alternative sampling strategies that incorporate large scale features of the energy landscape, while resisting the impact of roughness; these also outperform Random Walk Metropolis. Numerical experiments on these landscapes are presented that confirm our predictions. Open analysis questions and numerical challenges are also highlighted.

A multiscale modelling strategy based on molecular statics and molecular dynamics simulations in combination with transition state theory has been developed to determine the flow stress of Al-Cu alloy containing Guinier-Preston zones as a function of temperature. The flow stress is assumed to depend on two contributions. The athermal contribution is given by the Taylor model and only depends on the elastic constants of the Al alloy, the dislocation density and the Burgers vector. The thermal contribution can be calculated following transition state theory from the obstacle strength and the free energy barrier.

Molecular statics simulations in combination with thermal annealing was used to determine the obstacle strength (the stress necessary to overcome the obstacle in the absence of thermal energy) for edge dislocations interacting with Guinier-Preston zones in two different orientations (0 and 60 degrees). It was found that the introduction of thermal annealing after each strain increment was critical to avoid the locking of the simulations in local minima which leads to an overestimation of the obstacle strength and to unrealistic interaction mechanisms between the dislocation and the Guinier-Preston zone in both orientations. The rate at which dislocations overcome the Guinier-Preston zone was determined as a function of the applied strain and temperature from molecular dynamics simulations. This information was used to the determine the thermodynamic quantities that control this process, namely the activation energy, the activation entropy and the activation volume for both dislocation/ Guinier-Preston zone orientations.

The predictions of the model were compared with experimental data in the literature as well as new data obtained by means of micropillar compression tests in an Al - 4 wt. % Cu alloy naturally aged at ambient temperature to obtain a homogenous distribution of Guinier-Preston zones. They were in good agreement at ambient temperature although modelling results overestimated the flow stress at low temperature. Overall, this talk presents a coherent multiscale methodology to determine the effect of Guinier-Preston zones in the flow stress of Al alloys by means of atomistic simulations.

Glass is used in a wide variety of fields, from nuclear waste encapsulation to biomedicine. Computer simulation is often used to understand the connection between glass composition, structure and properties. Existing modelling typically uses Molecular Dynamics approaches to simulate the cooling of glass and these simulations can only simulate cooling rates down to about 10¹⁰K/s. This is several orders of magnitude faster than experimental cooling rates even for the most rapidly quenched glasses.

The cooling rate of the glass is known to affect the structure of the glass including the connectivity of the 3D structure and how ions are incorporated into the glass. The gap between experimental and modelling cooling rates thus make it difficult to compare results and to gain insight into the role of glass composition on its properties.

In this work we will discuss using a parSplice approach to model the cooling of glasses. We will discuss the differences observed in the glass structure when lower cooling rates are explored through this approach and the potential for the application of this methodology to other systems. The efficiency of the methodology and the limits this places on its applicability will also be discussed.

We derive and analyze a novel approach for modeling and computing the mechanical relaxation of incommensurate 2D heterostructures. Our approach parametrizes the relaxation pattern by the compact local configuration space rather than real space, thus bypassing the need for the standard supercell approximation and giving a true aperiodic atomistic configuration. Our model extends the computationally accessible regime of weakly coupled bilayers with similar orientations or lattice spacing, for example materials with a small relative twist where the widely studied large-scale moire patterns arise. Our model uses a generalized stacking fault energy for interlayer interactions and makes possible the simulation of the relaxation of multi-layer heterostructures for which a planar moire pattern does not exist.

We will present recent works which aim at understanding the mathematical foundations of first order kinetics and the harmonic transition state theory. The cornerstone of the analysis is the notion of quasi stationary distribution. This analysis sheds a new light on the accelerated molecular dynamics algorithms which have been proposed by Arthur Voter and co-workers in the late nineties. It also opens the route to generalizations of the original algorithms, in particular of the Parallel Replica algorithm.
References:
- G. Di Gesù, T. Lelièvre, D. Le Peutrec and B. Nectoux, Sharp asymptotics of the first exit point density, Annals of PDE, 5(1), (2019).
- C. Le Bris, T. Lelièvre, M. Luskin and D. Perez, A mathematical formalization of the parallel replica dynamics, Monte Carlo Methods and Applications, 18(2), 119–146, (2012).
- T. Lelièvre, Mathematical foundations of Accelerated Molecular Dynamics methods, In: Andreoni W., Yip S. (eds) Handbook of Materials Modeling. Springer, Cham, (2018).

One of the most stringent limitation of conventional molecular dynamics is the extremely short timescales that are amenable to direct simulation, even when massively-parallel ressources are invested in the calculation. The Parallel Replica Dynamics (ParRep) method, introduced by A. Voter, unlocks the potential of parallel architectures by allowing for an effective parallelization of the problem in the time-domain. However, in order to benefit from this extension of timescales, practitioners of ParRep need first address two issues: i) the identification of metastable sets on the energy landscape of the system, and ii) the accurate estimation of the time required to relax to the quasi-stationary distribution in each set. Both these problems can be readily addressed on smooth energy landscapes, but become extremely difficult for complex systems, where considerable prior knowledge often needs to be injected in the construction of the metastable states; a process which, even then, remains error prone.

We show how both these problems can be tackled simultaneously using data analysis techniques derived from the stochastic Koopman operator formalism. In this new method, states are locally defined and characterized on-the-fly by analyzing the trajectories produced as part of the ParRep simulation itself. To do so, the algorithm requires only the specification of a (potentially large) number of atomistic descriptors in terms of which the states will be implicitly defined. We demonstrate the new approach through long-time simulations of biomolecules and of soft defects in metals.

Metallic nanoclusters are functional materials with many applications owing to their unique physical and chemical properties, which are sensitively controlled by their shapes and structures. An in-depth understanding of their morphology stability is therefore of crucial importance. It has been well documented by transmission electron microscopy (TEM) studies that metallic nanoclusters can interconvert between different isomers. However, the relevant mechanisms remain elusive because the timescales of such shape fluctuations are too short to be resolved experimentally and yet too long for conventional atomistic simulations. By employing a recently introduced Accelerated Molecular Dynamics method, Parallel Trajectory Splicing, we present simulations that reached timescales of milliseconds and thus provide a clear description of the dynamic process of the experimentally observed shape fluctuation in metallic nanoclusters. We observe transformations back and forth between face-centered-cubic (fcc) and structures with five-fold symmetry (decahedron or icosahedron). These transitions occur following either by a partial-dislocation-mediated twinning mechanism or by a surface-reconstruction driven process. The identified pathway is in remarkable agreement with the existing microscopy results and serves as further evidence that shape fluctuation can occur directly through thermal activation, without involving melting or other external factors.

A range of specialized Molecular Dynamics (MD) methods have been developed in order to overcome the challenge of reaching longer timescales in systems that evolve through sequences of rare events. In this talk, we consider Parallel Trajectory Splicing (ParSplice) which works by generating large number of trajectory segments in parallel in such a way that they can later be assembled into a single statistically correct state-to-state trajectory, enabling parallel speedups up to the number of parallel workers. In practice, the ability for ParSplice to scale significantly improves when it is possible to predict where the trajectory will be found in the future. With this insight in mind, we develop a maximum likelihood transition model that is updated on the fly and make use of an uncertainty-driven estimator to approximate the optimal distribution of trajectory segments to be generated next.

This work was supported by the U.S. Department of Energy Office of Science Graduate Student Research (SCGSR) program, administered by the Oak Ridge Institute for Science and Education under DE SC0014664 and by the Exascale Computing Project (17-SC-20-SC).

The mechanical behavior of materials operating under high temperatures is strongly influenced by creep mechanisms such as dislocation climb which are controlled by the diffusion of vacancies. However, atomistic simulations of these mechanisms have traditionally been impractical due to the long time scales required. To overcome these time scale challenges, we use Parallel Trajectory Splicing (ParSplice), an accelerated molecular dynamics method, to simulate dislocation climb in nickel. We focus on modeling the activity of a vacancy near a jog on an edge dislocation in order to observe vacancy pipe diffusion and vacancy absorption at the jog. Using ParSplice, we are able to achieve a total simulation time of over 4 μs for a system containing approximately 16,000 atoms, during which over 2000 vacancy diffusion and absorption events occurred. Based on these results, we estimate average rates for pipe diffusion and vacancy absorption into the jog; results are in good agreement with molecular static calculations.

Capturing plasticity at realistic dislocation densities with high configurational complexity requires a continuum-level discrete dislocation dynamics (DDD) description. However, many features controlling dislocation motion are inherently atomistic, such as the interaction of dislocations with solutes, precipitates, cracks, and interfaces. Here, we review a recent multiscale concurrent coupled atomistics/discrete-dislocation model in 3d that couples the DDD code Paradis and the MD code LAMMPS [1]. The need for precise matching of continuum dislocation properties to atomistic dislocation properties is demonstrated, and issues in making this connection are also discussed. A full coupling of atomistic and continuum domains requires solution of an elasticity problem in the continuum domain. Standard methods such as FE are computationally infeasible. We thus present a Greens function implementation where H-matrices and a Newton-Krylov solution scheme are employed to achieve computational efficiency with controllable accuracy [2].
[1] G. Anciaux et al., “The coupled atomistic/discrete dislocation method in 3d. Part I: Concept and algorithms”, J. Mech. Phys. Solids 118, p. 152-171 (2018); M. Hodapp et al., “Coupled atomistic/discrete dislocation method in 3d Part II: Validation of the method”, J. Mech. Phys. Solids 119, p. 1-19 (2018); J. Cho et al., “The coupled atomistic/discrete dislocation method in 3d. Part III: Dynamics of hybrid dislocations”, J. Mech. Phys. Solids 118, p. 1-14 (2018).

[2] M. Hodapp et al., Computer Methods in Applied Mechanics and Engineering 348, 1039-1075 (2019).

Novel 2D materials have unusual properties, many of which are coupled to their large scale structural properties. Microscopic modeling of structure and transport is a formidable challenge due to a wide span of length and time scales involved. I will review recent progress in structural multiscale modeling of 2D materials and thin heteroepitaxial overlayers [1], and graphene and h-BN [2,3] in particular, based on the Phase Field Crystal (PFC) model combined with Molecular Dynamics and Quantum Density Functional Theory. The PFC model allows one to reach diffusive time scales at the atomic scale, which facilitates quantitative characterisation of domain walls, dislocations, grain boundaries, and strain-driven self-organisation up to micron length scales. I will discuss how PFC generated multigrain systems can be further used for Molecular Dynamics simulations of thermal conduction in selected 2D materials [3,4].

[1] K. R. Elder et al., Phys. Rev. Lett. 108, 226102 (2012); Phys. Rev. B 88, 075423 (2013); J. Chem. Phys. 144, 174703 (2016).
[2] P. Hirvonen et al., Phys. Rev. B 94, 035414 (2016).
[3] H. Dong et al., Chem. Phys. Phys. Chem. 20, 4263 (2018).
[4] Z. Fan et al., Phys. Rev. B 95, 144309 (2017); Nano Lett. 7b172 (2017); Phys. Rev. B 99, 064308 (2019); K. Azizi et al., Carbon 125, 384 (2017); K. Xu et al., Mater. Sci. and Eng. 26, 8 (2018); Phys. Rev. B 99, 054303 (2019).

*This work has been supported in part by the Academy of Finland QTF Centre of Excellence program (project 312298).

For n-type transparent conducting oxide (TCO) materials such as SnO2, In2O3 and ZnO, native defects play a key role in electronic conductivity. Depending on their electronic structure, energetics and geometries, defects can act as donors, resulting in intrinsic n-type conductivity, or can compensate extrinsic donors such as Sn in In2O3. Predictive modelling of the properties of defects in such systems requires a detailed description of the dielectric response of the host material, which can be difficult to obtain using standard supercell techniques. Here, we employ the hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster method, a multi-region approach that allows us to model defects at the true dilute limit, with polarisation effects described in an accurate and consistent manner. Moreover, we develop techniques to analyse the energetic balance between electrons bound to donors in diffuse and compact states, a difficult problem regardless of the model employed. We benchmark our results where possible and find good agreement with experiment for a variety of defect-related properties.

In long time scale simulations based on the adaptive kinetic Monte Carlo (AKMC) method (sometimes referred to as on-the-fly KMC, or self-learning KMC) the essential task is to find all relevant, low lying first order saddle points on the energy ridge surrounding a given initial state minimum without preconceived notion of the possible final states. In current implementations of the method, several random initial displacements from the minimum are used as starting points for climbs up the energy surface to sample the saddle points. While schemes have been developed to improve the sampling and reduce the occurrence of reconvergence on known saddle points [1], a significant concern are saddle points that are hard to reach for example when the basin of attraction does not reach the boundary of the convex region around the local minimum, and the waisted computational effort when the same saddle point is found repeatedly. A systematic approach where a scalar quantity based on semi-local information about the energy surface is used to move from a local minimum along a linear path with bifurcation points so as to map out all first order sadddle points will be presented [2]. Also, the use of Gaussian process regression in this context to reduce the number of energy and atomic force evaluations will be discussed, analogous to what has already been reported for nudged elastic band calculations [3].

[1] M.P. Gutiérrez, C. Argáez and H. Jónsson, Journal of Chemical Theory and Computation 13, 125 (2017).

[2] V. Ásgeirsson, K. Blöndal, A. Laio and H. Jónsson (in preparation).

[3] O-P. Koistinen, V. Ásgeirsson, A. Vehtari and H. Jónsson, J. Chem. Phys. 147, 152720 (2017).

Aluminum is well-known for its exceptional resistance to high-temperature corrosion via oxidation. This is attributed to the formation of a scale, i.e. a thin oxide film on the surface. This protective layer acts as a barrier to diffusion and prevents further deterioration. However, despite tremendous efforts, the mechanisms that underlie the slow diffusion through the scale are poorly understood. Unfortunately, this limits the development of highly corrosion resistant materials that are required, for instance, in the high-temperature environments encountered in the petrol industry.

In this work, we lay the foundation for further computational research on diffusivity of defects by computing formation energies of point defects in aluminum oxide from first-principles. We specifically explore the influence of exchange and correlation by comparing standard density functional theory to hybrid-functional results and we also study the influence of spin polarization. This allows us to confirm previous studies showing that a peroxide split interstitial defect is energetically favorable to the traditionally considered, quasi-atomic oxygen interstitial defect at the octahedral site. Interestingly, we also find such a displaced ground-state configuration for aluminum interstitials. In order to efficiently sample possible starting atomic geometries for these as well as extrinsic interstitials, we explore using Bayesian analysis. In our work, we also investigate the viability of octahedral configurations for extrinsic interstitials and we examine the impact of defect charge state on the ground-state geometry. These so far overlooked lower-energy geometries possibly point at the possibility of a more general mechanism in this material for interstitial formation and can, hence, influence defect migration paths and predictions for diffusion properties.

Transport properties of point defects are crucial physical quantities governing the microstructural evolution of metals and alloys. Characterizing the kinetics of defects in presence of energetic and entropic traps is notoriously difficult due to the emergence of well separated time scales in the simulations. We herein describe fast first-passage algorithms based on the theory of absorbing Markov chains assuming that defects undergo reversible diffusion. We show that the absorbing transition rate matrix can be transformed into a symmetric definite-positive matrix and that this similarity transformations considerably simplifies the computational tasks. The efficiency of the approach is demonstrated in simulations of diffusion and aggregation of Copper and Manganese clusters in $\alpha$-iron. Substantial accelerations are achieved by directly sampling the important rare events contributing to atomic transport.

Modeling physical systems as stochastic processes helps capture natural fluctuations and improves long-time accuracy. Yet, the physical processes which may be confidently formulated in terms of an easily computable stochastic process are limited to those which are driven exclusively by either rare events or Gaussian white noise. Additionally, the stochastic integration schemes used in the development and implementation of these algorithms have been traditionally based on solving the Itô and Stratonovich integrals, which do not preserve the chain rule for differentiation or the Newton-Leibniz integration rule when applied to Lévy processes driven by combined rare events and continuous-in-time Brownian motion, e.g. liquid-solid phase transformations. In this presentation we explore the consequences of solving the Marcus canonical integral, which preserves the aforementioned rules when applied to Lévy processes, to stochastic simulations of material phase transformations.

The recent development of machine learning interatomic models (IMs) provides a promising avenue to extend the scope of first principles methods to materials problems on larger length and time scales. Machine learning IMs interpolate a training set of density functional theory (DFT) calculations and can provide close to DFT accuracy at a fraction of the computational cost. A challenge is the inherent lack of transferability of machine learning IMs; their predictions "far" from the training set are suspect. Using 2D materials as a test system, we address this limitation in two ways: (1) we develop a hybrid IM that combines the known physics of long-range dispersion with a neural network (NN) to model short-range interactions for multilayer graphene systems; (2) we apply an NN dropout technique to assess the uncertainty in predictions of the IM to prevent its use outside its range of validity. These techniques along with tools for automated training have been implemented within a general-purpose fitting framework called KLIFF (KIM-Based Learning-Integrated Fitting Framework). KLIFF is an open source Python package compatible with the Knowledgebase of Interatomic Models (KIM) application programming interface (API) standard, allowing IMs fit with it and distributed through https://openkim.org to be used by a wide variety of simulation codes compatible with the KIM API.

Atomistic level investigations are required for the discovery of novel material properties which, in turn, results in advancements in technology such as batteries, capacitors, etc. Most of the atomistic-scale modeling methods that have been developed to conduct these investigations are restricted to small systems due to associated computational costs. On the other hand, it is possible to simulate large atomistic systems using ReaxFF, which is a very reliable and transferable reactive force field concept used by a broad academic and industrial user community. In ReaxFF, inter and intra-atomic interactions are calculated by functionals that involve system-specific parameters. These parameters must be optimized for each materials system, a task requiring exploration of multi-dimensional parameter landscape to obtain a high level of accuracy. However, the current ReaxFF optimization method is a sequential search algorithm which is limited in terms of exploring a wider parameter space. Moreover, due to the consecutive nature of the method, the optimization algorithm is non-parallelizable. In this respect, this study develops an initial design-based machine learning (ML) framework to enable a fast and high-quality force field development to improve the ReaxFF optimization algorithm.

Initially, a Latin Hypercube Design (LHD) algorithm is used to explore the parameter landscape extensively. The LHD passes the information about explored regions to the ML method for model training. The ML method is composed of regression and classification models. A deep learning model finds the minimum discrepancy regions, while a classification model classifies the physically-feasible/unfeasible regions in parameter landscape. This ML method eliminates unfeasible regions, which originate from the unphysical atomistic interactions, and constructs a more comprehensive understanding of a physically meaningful parameter landscape. In addition, the developed framework is easily parallelizable; thus, it will significantly reduce the force field optimization time while increasing the quality by benefiting from ML. The performance of the framework will be evaluated by its application to systems with available force field parameter sets including 2D materials and alloys.

Kinetic Monte Carlo (KMC) is the method of choice for studying diffusion in solid materials. The method is often limited in accuracy by the quality of the migration energy barriers used to describe the system of interest. Here is KMC's tradeoff between precision and computational cost - a more detailed description of the energy barriers will typically require either heavier precalculations, or resolving barriers on-the-fly.

Like in the realm of potential energy functions, machine learning (ML) offers a way to bridge the computational gap in this problem also. We have expanded a KMC model for the Cu surface, earlier developed in our group, by employing ML for the purpose of avoiding some of the barrier calculations. Our most recent pursuit is for an active learning scheme, which as a first stepping stone requires the ability to estimate uncertainties of the predicted barrier values. For some ML methods, such as Gaussian process regression, this property comes out-of-the-box, while for others, like artificial neural networks, we have had to resort to more specialized tricks. In this talk, we present the progress we have made in estimating the uncertainty of various regression models, as well as the developing framework around the active learning scheme — when to trust the regressor during the KMC simulation, and when should the oracle be queried for more precise data.

Cross-linked Polyethylene (XLPE) has been recognized as an outstanding insulator for high-voltage power cables due to its favorable structural integrity at high temperature, lower moisture sensitivity, chemical resistance, and lower rates of failure due to aging. Dicumyl Peroxide (DCP) is the main PE crosslinking agent, however, DCP-initiated crosslinking leads to the formation of different byproducts such as cumyl alcohol, acetophenone, methylstyrene, water and methane that remain in the final product and may have adverse effects on cable function and its long-term properties. Previous studies stated that the byproducts such as cumyl alcohol (CA), acetophenone (AP), α-methylstyrene and antioxidant additives can cause a space charge formation in the final product and affect electrical properties of the XLPE. Currently, our knowledge is limited about the effects of the byproducts on the properties of the final XLPE product. By understanding on how the crosslinking byproducts change the XLPE properties, improvements may be formulated relative to the conventional XLPE cable. In order to design and optimize the XLPE cables, it is crucial to obtain detailed, atomistic-scale insight of the XLPE chemistry since each and every byproduct in the XLPE affect differently to the polymer’s electrical properties, thus the effects of each of the byproducts should be investigated. Quantum Mechanic (QM) methods are considered the most accurate technique among atomistic simulation methods, however, due to high computational costs, these first-principles techniques can only be viable for relatively small length and short time scales events. The ReaxFF reactive force field method can model chemical reactions based on the bond order concept without the full expense of QM methods, and it has been successfully applied to a wide range of systems such as polymeric fibers, covalent, metallic, and metal oxide/hydride/carbide materials. e-ReaxFF method is an extension of the ReaxFF method with description of an explicit electron-like or hole-like particle. In this study, we used e-ReaxFF technique to investigate the role of XLPE byproducts, and processing variables on the dielectric strength of the PE. To achieve this, we developed an e-ReaxFF reactive force field verified against DFT data obtained for XLPE chemistry. Then, using this force field, we performed MD simulations to find the effects of different byproducts, temperature, density, the ratio of byproducts to PE on the dielectric strength of the PE.

Molecular dynamics (MD) simulations with Reactive force field (ReaxFF) provide valuable insight due to accurate description of chemistry of the material. This comes with a computational cost which limits size and time scale of the simulation. Non-reactive models, such as OPLS, may provide accurate description of the dynamics of the system without chemical reactions. Typically OPLS runs fifty times faster compared to ReaxFF. We developed a hybrid framework which utilize ReaxFF for chemical reactions while using non-reactive force field to accelerate MD simulation. In this framework, chemically active regions such as surfaces, radical sites, described with ReaxFF and rest of the system described with OPLS. A tracking tool utilized to define reactive and non-reactive regions during the simulation. ReaxFF uses a unique parameter set for each atom, bond, angle and dihedral types, however OPLS parameters depends on local bond topology. Since ReaxFF allows chemical reactions, bond topology of the system changes during the simulation. We used a Machine Learning based model to predict OPLS parameters in order to dynamically transfer atoms from reactive to non-reactive regions. As a test case, we consider cross-linking of poly ethylene with dicumyl peroxide with large scale hybrid simulations. We achieved significant improve in the computational efficiency and were able to simulate realistic size systems.