Modeling of nanostructure evolution in solids requires comprehensive data on the properties of intrinsic point defects, foreign atoms and defect clusters. Since most processes occur at elevated temperatures not only the energetics of these species in the ground state but also their temperature-dependent free energies must be known. These data can be used to obtain improved, temperature-dependent input parameters for atomistic or object kinetic Monte Carlo simulations and rate theory.
The first-principles calculation of contributions of phonon and electron excitations to free formation, binding, and migration energies is illustrated in the case of bcc-Fe. First of all, the ground state properties of the defects are determined under constant volume (CV) as well as zero pressure (ZP) conditions, and relations between the results of both kinds of calculations are discussed. Second, vibrational and electronic contributions to defect free energies are calculated using the equilibrium atomic positions determined in the ground state for the CV and the ZP case. Additionally, the quasi-harmonic approach is applied to ZP-based data in order to obtain results closest to the experimental conditions at elevated temperatures. However, in most cases considered this leads only to minor modifications. In contrast to ground state energetics the CV- and ZP-based defect free energies do not become equal with increasing supercell size. A simple transformation is found between the CV- and ZP-based frequencies and between the corresponding free energies. Finally, self-diffusion via the vacancy mechanism is investigated. The ratio of the respective CV- and ZP-based results for the vacancy diffusivity is nearly equal to the reciprocal of that for the equilibrium concentration. This behavior leads to almost identical CV- and ZP-based values for the self-diffusion coefficient. Obviously, this agreement is accidental and cannot be generalized to other cases.The consideration of the temperature dependence of the magnetization yields self-diffusion data in very good agreement with experiments.

First principles quantum mechanics based modeling schemes, such as density functional theory (DFT), are powerful tools to study the static and dynamical evolution of processes (e.g., chemical reactions, phase transformations, transport). Nevertheless, such methods have several practical drawbacks. For instance, owing to the computational expense, a typical DFT-based simulation can only handle system sizes of at most 100s of atoms or span timescales of the order of just picoseconds (with a maximum reachable timescale of about a nanosecond). Here, we show that DFT-based simulations can be significantly accelerated, by using machine-learning methods to make rapid high-fidelity predictions, based on past data or knowledge.
For many atomistic simulations, e.g. molecular dynamics, geometry optimization, identifying reaction barriers, materials properties, etc., determining the force on an atom is key. If we are able to rapidly estimate (with acceptable accuracy) the atomic forces of a new configuration, given past similar configurations (determined, say, using DFT), then we can significantly speed up the simulations. The ‘training&’ underlying this force prediction capability requires a critical amount of prior information or results, adequate descriptors (or fingerprints) that uniquely represent our configurations, and a suitable measure of (dis)similarity between configurations. We show here that such training, prediction, and consequently, the acceleration of DFT-based simulations, are indeed possible. Illustrations of this new development are made for the case of the self-diffusion of vacancies in bulk Al, the diffusion of an Al adatom on an Al surface, finding optimal geometries of bulk Al structures with more than 100s of atoms, identifying reaction barriers for simple elementary processes, and lastly, determining thermodynamic properties of bulk Al. Such a framework lays the foundations, for understanding the behavior of chemical reactions or materials at increased time and length scales, all at first principles accuracy, a capability that is well beyond the current realm of first principles methods.

Theoretical calculations have been proven to be very useful tool in understanding the fundamental principles of the structural and electronic properties of energy relevant materials. Here we discuss the fundamental ab-initio calculation schemes used and overview examples of how these wave-based and density functional theory (DFT) methods relates to experiments.[1,2,3] We also briefly present hybrid methods enabling time-resolved surface reconstructions of semi-conductor nanoparticles in solution up to 5 ns. Some fundamental problems in applying the time-independent Schrodinger or Kohn-Sham equations to material properties are exemplified where the time-dependent Dirac equation instead is preferred. In the view of the success of quantum field theory and the idealized environment approximation (IEA) with one or few particle excitations, we present a high-fidelity numerical solution to the time-dependent Dirac equation [4] and show that this can be calculated at a lower computational cost than the time-dependent Schrodinger equation utilizing the hyperbolicity of the Dirac equation. We finally present time-resolved calculations of Klein tunneling in 3D and spin selective reflectance [4] at both idealized potential barriers and situations approaching material systems.
[1] Lissau j.S., Nauroozi D., Santoni M.-P., Edvinsson T., Ott S., Gardner J. M., Morandeira A.
What Limits Photon Upconversion on Mesoporous Thin Films Sensitized by Solution-Phase Absorbers? J. Phys. Chem. C 2015, 119, 4550minus;4564
[2] Park B., SM Jain S.M., Zhang X., Hagfeldt A., Boschloo B., Edvinsson T.
Resonance Raman and Excitation Energy Dependent Charge Transfer Mechanism in Halide-Substituted Hybrid Perovskite Solar Cells, ACS Nano, 2015, 9, 2088-2101
[3] Unger, E. J., Edvinsson, T., Roy-Mayhew, J. D., Rensmo, H, Hagfeldt, A., Johansson E. M. J. and Boschloo, G. , Excitation energy dependent charge separation at holetransporting dye/TiO2 hetero-interface, J. Phys. Chem. C , 2012, 116, 21148minus;21156
[4] Almquist, M., Mattsson, K.; Edvinsson, T.
High-fidelity numerical solution of the time-dependent Dirac equation
J of Comput. Phys, 262, 2014, 86-103

In this talk, multiscale simulation will briefly be introduced as an enabling technology that bridges the gap between scales. A specific application of multiscale simulation is the prediction of macroscopic behavior from first principles. A more impactful avenue of research is how one could use multiscale modeling in reverse engineering for predicting new materials. We will demonstrate how descriptor-based modeling can enable such a search of novel materials and assess this framework with experiments. We will demonstrate this methodology for the ammonia decomposition on bimetallic catalysts using first principles kinetic Monte Carlo (KMC) simulations. We will address challenges rooted at efficient parameterization of the KMC engine and the lack of efficient tools for carrying out sensitivity analysis of complex reactions on surfaces. We will show that the microstructure of materials plays a profound role and poses interesting optimization problems.

Alloying allows the properties of a material to be tuned to a specific application, and advanced materials for a wide variety of applications are often alloys. Radiation-induced point defects can modify and degrade material properties. In an alloy, defect energies are sensitive to the occupations of nearby atomic sites and thus vary with location in the alloy, which leads to a distribution of defect properties. When radiation-induced defects diffuse from their initially non-equilibrium locations, this distribution becomes time-dependent. Furthermore, the defects can become trapped in energetically favorable regions of the alloy leading to a diffusion rate that slows dramatically with time. Density Functional Theory (DFT) allows the accurate determination of ground state and transition state energies for a defect in a particular local environment in the alloy but requires thousands of processing hours for each such calculation. Kinetic Monte-Carlo (KMC) can be used to model defect diffusion and the changing distribution of defect properties but requires energy evaluations for millions or billions of local environments. We have used the Cluster Expansion (CE) formalism to “glue” together these seemingly incompatible methods in order to model defect diffusion in alloys. In the CE approach, the occupation of each alloy site is represented by an Ising-like variable, and products of these variables are used to expand quantities of interest. Once a CE is fit to a training set of DFT energies, it allows very rapid evaluation of the energy for an arbitrary configuration, while maintaining the accuracy of the underlying DFT calculations. These energy evaluations are then used to drive our KMC simulations. We will demonstrate the application of our DFT/MC/KMC approach to model thermal and carrier-induced diffusion of intrinsic point defects in InGaAs.
Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under Contract DE-AC04-94AL85000.

Graphene, and its recently discovered layered “cousin” materials such as h-BN, MoS₂, silicene, topological insulators, and phosphorus, have been occupying the central stage of modern materials science. Each important member in this big family distinguishes itself with its unique properties, and different structural combinations of any two members can further enrich their functionality for broader technological significance. Aside from their uniqueness, these materials also share easily identifiable commonalities, especially the predominant role of van der Waals forces in their interlayer or interface bonding. Therefore, an in-depth understanding of each individual member could be greatly beneficial in our developing a more complete view of the whole materials family. Here we use multiscale modeling approaches combining first-principles calculations and nonequilibrium rate equation analysis to elucidate the atomistic growth mechanisms of graphene and related materials on metal substrates. Our studies (a) reveal why Cu substrates are highly preferred to catalyze mass production of quality graphene, (b) identify a novel kinetic pathway towards effective suppression of undesirable grain boundaries in graphene growth on Cu(111) via chemical vapor deposition, (c) propose the use of aromatic molecules as carbon sources for graphene growth at dramatically reduced growth temperature, (d) identify carbon dimers as the dominant feeding species in the epitaxial growth of graphene on both Cu(111) and Cu(100) substrates, and (e) exploit the conditions at which second layer of graphene can be grown underneath the first layer. We finally contrast the different growth behaviors of graphene and h-BN, and the main findings will be discussed in close comparison with experiments.

Far more challenging than either fully lithiated or delithiated states is the intermediate state of delithiation, in which two kinds of solid solution order may coexist in cathode oxides. One is Li/vacancy order and the other one is charge-order of transition metal cations in mixing valences. The former may introduce phase transition such as O3-O1 stacking transition during delithiation of LiCoO₂ cathode compound. The latter, which is coupled with the former, may impact the electronic and ionic transportation.
From the point view of Pauling's rules of electrostatic valence and ionic crystal as a periodic piling of coordination polyhedra, cathode compound can be regarded as a network of neutral polyhedron unit (NPU), in which the charge of anion at vertex is fractioned by coordination number to compensate the charge of cation at center. Based on counting of pairwise interaction between NPUs, we proposed a new method to calculate Madelung constant, expressed in matrix form.
What Madelung matrix describes is theCoulomb interaction betweenpolyhedral building-blocks, relating to Pauling&’s third rule about polyhedral combination. By investigation of Madelung matrices for eutectic HCP and CCP lattices, we reveal the physical rule indicating the electrostatic stability of edge-sharing between tetrahedra and providing theoretical validation of the structure stability of multi-lithium compounds, which violates Pauling's third rule. Interestingly, we found the O3-O1 phase transition is electrostatically spontaneous, driving by the inter-layer octahedral interaction.
This new method is adeptto analyze order configurations of any specific charge distribution on cation sites, includingvariable stoichiometry (intermediate phases during deintercalation of cathodes), as well as variable polymorphisms (such as A₂BSiO₄ silicates), can be extended to spin-polarized systems, e.g., AFM (anti-ferromagnetic) order. As an example, we identified 25 classes of charge-order configurations within Fd-3m lattice, predicting local charge-order configurations coexisting rather than global charge-order in LiMn₂O₄.

Evaluating mechanical properties of lightweight metal alloys for fuel-efficient vehicles, investigating liquid lithium films for fusion reactor walls, and studying charge-discharge cycles of next generation Li-ion battery anodes are three projects that may appear to have nothing in common other than that they could be part of a larger energy research portfolio. However, all three exploit a quantum mechanics method - orbital-free density functional theory (OFDFT) - that directly evaluates electron distributions instead of wavefunctions. This technique is orders of magnitude faster than standard Kohn-Sham DFT because it scales quasilinearly with a small prefactor. As such it can be used to study many thousands of atoms with quantum mechanics, or to perform longer time scale ab initio molecular dynamics on smaller samples. Consequently, OFDFT is able to explicitly study, e.g., plasticity in metals and liquid metal dynamics. I will give a brief history of our work in this field and then present our recent advances in OFDFT methods and applications that now furnish accurate treatment of semiconductors and transition metals, extending the reach of OFDFT nearly to the full periodic table. Most recently, our code PROFESS has been enabled to run on petascale platforms. The workhorse routine of PROFESS is the Fast Fourier Transform (FFT), which in its typical form is limited to operate across 1-2 thousand cores. We have extended Wang&’s small box FFT idea for the Hartree kernel (X.-W. Jiang, S.-S. Li and L.-W. Wang, Comp. Phys. Comm., 184, 2693 (2013)) to work for all terms, even the nonlocal kinetic energy density functional. The code now scales across 100,000 cores, opening the door to first principles simulations of many millions of atoms.

Accurate and detailed electronic structures are prerequisites for our understanding of materials properties. Ideally, one just needs to solve the Schroedinger equation which has been introduced for over 80 years. Unfortunately, the many-body nature of the Schroedinger equation makes itself extremely difficult to solve. Theories of various levels of accuracy exist in the literature. Very accurate methods, such as the configuration interaction method, often have a prohibitive computational cost that scales exponentially with system sizes. Efficient methods, such as the Kohn-Sham density functional theory, often have large errors that are difficult to determine and control. All these difficulties severely limit the predictive power of materials simulations. A promising way to overcome this obstacle is to perform quantum mechanics embedding simulations, in which the key region in materials is solved by highly accurate methods, with the unimportant regions to be solved by fast, physically correct methods. In this seminar, I will present our recent advances in quantum mechanics embedding theory.[1,2] I will show how to perform quantum mechanics embedding simulations of materials in a seamless and first-principle way. I will also discuss its application to two long-term puzzles that are related to surface catalysis and metal corrosion: (a) the adsorption of carbon monoxide on copper surface [2] and (b) the counterintuitive oxidation process of aluminum surface. [3]
[1] C. Huang and E.A. Carter, J. Chem. Phys., 135, 194104 (2011).
[2] C. Huang, M. Pavone, and E. A. Carter, J. Chem. Phys., 134, 154110 (2011).
[3] F. Libisch, C. Huang, P. Liao, M. Pavone, and E.A. Carter, Phys. Rev. Lett., 109, 198303 (2012).

We present our recent development and applications of computational techniques for local energies and local stresses based on density-functional theory (DFT) [1-3]. In conventional plane-wave DFT methods, total energies and stress tensors are given as integrated or averaged quantities, and local distributions of energies or stresses are not obtained. Previously, the schemes to calculate energy density and stress density were proposed [4], while the inherent gauge-dependent problem associated with the choice of symmetric or asymmetric forms of the kinetic terms prevented practical applications. Within the PAW-GGA framework, we have developed the computational techniques to obtain local energy and local stress as unique physical quantities via integrating the energy and stress densities inside local atomic or layer regions to satisfy the gauge-independent conditions [1, 2]. In this talk, we present our recent applications [2, 3] to i) grain boundaries (GBs) in Fe, Al and Cu, ii) impurity segregation to GBs, iii) tensile tests of GBs, iv) local bulk moduli of Fe-Si alloys, and v) Fe/TiC interfaces. About GBs, we have observed the presence of tighter sites with smaller atomic volumes and looser sites with larger atomic volumes. The former sites show compressed stresses and lower local energies, while the latter sites show tensile stresses and higher local energies. These local features have correlation with local magnetic moments in Fe GBs. For impurity segregation, the local-energy decomposition is quite effective to analyze the segregation mechanism. For Si and Mg impurities in Fe, Al and Cu GBs, we have observed that the tighter and looser sites show quite different segregation energies and segregation mechanisms, depending on impurity and metal species. For the first-principles tensile tests of metallic GBs, the variations of local energies and local stresses of all the atoms during the tensile deformation and failure provide valuable information on the mechanical properties of GBs. By using local stress changes for small volumetric compression and tension of supercells, we can obtain local bulk modulus of each atom or atomic group in Fe-Si alloys, which are useful to understand the origin of peculiar bulk-modulus changes of Fe-Si alloys depending Si concentration. For Fe/TiC coherent interfaces, the local-energy and local-stress distributions are useful to understand the origins of interface adhesion and interface misfit stresses. Finally, we discuss future applications of our schemes and the physical meaning of local energy and local stress.
[1] Y. Shiihara et al., Phys. Rev. B 81, 075441 (2010); ibid.87, 125430 (2013).
[2] H. Wang et al., J. Phys.: Condens. Matter 25, 305006 (2013); submitted (2015).
[3] S. Kr. Bhattacharya et al., J. Phys.: Condens. Matter 25, 135004 (2013); ibid.26, 355005 (2014); submitted (2015).
[4] N. Chetty and R.M. Martin, Phys. Rev. B 45, 6074 (1992); A. Filippetti and V. Fiorentini, ibid.61, 8433 (2000).

A subspace formulation of time-dependent density functional theory (TDDFT) is proposed for large-scale calculations based on density functional perturbation theory. The formulation is implemented in conjunction with projector augmented-wave method and plane-wave basis set. A key bottleneck of conventional TDDFT method is circumvented by projecting the time-dependent Kohn-Sham eigenvalue equations from a full Hilbert space to a substantially reduced sub-Hilbert space. As a result, both excitation energies and ionic forces can be calculated accurately within the reduced subspace. The method is validated for several model systems and exhibits the similar accuracy as the conventional TDDFT method but at a computational cost of the ground state calculation. The Born-Oppenheimer Molecular Dynamics can be successfully performed for excited states in C₆₀ and T₁₂ molecules, opening doors for many applications involving excited state dynamics.

We are developing a modern framework for quantum based molecular dynamics simulations that combines some of the best features of regular Born-Oppenheimer and Car-Parrinello molecular dynamics. The new framework is based on an extended Lagrangian formulation of Born-Oppenheimer molecular dynamics that allows, for the first time, efficient energy conserving Born-Oppenheimer molecular dynamics simulations with a computational complexity that scales only linearly with the system size. Often only a single diagonalization per time step is required. The efficiency and accuracy of the new dynamics can be understood from a backward analysis. Instead of integrating an underlying exact dynamics with approximate forces, exact forces that do not rely on the fulfillment of the Hellmann-Feynman theorem are used to integrate the equations of motion for an approximate shadow Hamiltonian. In this way properties such as the total energy can be controlled. This geometric approach to integration is widely used in classical molecular dynamics, e.g. in the velocity Verlet algorithm. Our new framework allows this geometric technique to be applied also to self-consistent field theory. Extended Lagrangian Born-Oppenheimer molecular dynamics represents a general approach that can be applied to a broad variety of quantum based molecular dynamics applications.

In this work we introduce RESCU, a powerful MATLAB-based Kohn-Sham density functional theory (KS-DFT) solver. We demonstrate that RESCU can compute the electronic structure properties of systems comprising many thousands of atoms using modest computer resources, e.g. 16 to 256 cores. For example, RESCU converged the total energy and electronic density of a 5,832 Si atoms supercell in 5.5 hour, a 8,788 Al atoms supercell in 23.9 hours, a 1,372 Cu atoms supercell in 9.1 hours and a small DNA molecule submerged in 1,713 water molecules for a total 5,399 atoms in 9.6 hours using 256 Xeon E5-2670 cores. Its computational efficiency is achieved from exploiting four routes. First, a vector space spanning the occupied Kohn-Sham subspace is built using the Chebyshev filtering technique proposed by Zhou et al.[1] Second, we use a hybrid numerical atomic orbital (NAO) plus real space grid scheme to generate efficiently an initial subspace which is crucially needed in the Chebyshev filtering paradigm. Third, we relax the computational and memory requirements associated with the Kohn-Sham orbitals by discretizing them on a coarser grid than the potentials and non-local projectors as advocated by Hirose and Ono[2]. Finally, by judiciously analyzing and optimizing various parts of the procedure in RESCU, we delay the O(N³) scaling to larger N, and our tests show that RESCU scales consistently as O(N^2.3) from a few hundred atoms to more than 8,000 atoms when using a real space grid discretization. The scaling is better or comparable in a NAO basis up to the 14,000 atoms level. Additionally, we developed a partial Rayleigh-Ritz algorithm to achieve efficiency gains for systems comprising more than 10,000 electrons. Our results suggest that the RESCU method has reached a milestone of solving thousands of atoms by KS-DFT on a modest computer cluster.
[1] Y. Zhou, et al., Phys. Rev. E 74, 066704 (2006)
[2] T. Ono, K. Hirose, Phys. Rev. Lett. 82, 5016 (1999)

In this talk, I will present our new development in the divide-and-conquer linear scaling three dimensional fragment (LS3DF) method. The LS3DF divides a large system into fragments, then patches the quantities (charge density, and kinetic energy) of the fragments to yield the values of the global system. LS3DF has been used to simulate systems with tens of thousands of atoms, e.g., a Moire&’s patterned double layers of 2D materials and electrostatic fluctuation in the hybrid perovskite system. Recently, we have implemented the LS3DF method with the GPU machine, and achieved significant speed up of the calculations. This will allow efficient atomic relaxation and molecule dynamics simulations. In another topic, I will also present a new method to accelerate the atomic relaxation using density functional theory (DFT). This is by on-the-fly fitting of the classical force field, and uses it as a guide to accelerate the DFT atomic relaxation. 3 to 6 times speedups can often be achieved using this method.

The implementation of a combined quantum mechanics and molecular mechanics (QM/MM) scheme, employing Kohn-Sham density functional theory (KS-DFT) and a single-center multipole expansion (SCME) model [1], is described. Both the electronic density of the QM part as well as the induced dipoles and quadrupoles of the MM part are solved for self-consistently under their mutual influence. The effect of the multipoles on the electronic density is represented as an additional external potential term which enters the conventional KS-DFT calculation, whereas the resulting electronic density is included in the SCME models as additional electric field terms which perturb the dipoles and multipoles.
SCME explicitly includes the instantaneous dipole as well as quadrupole polarizability of the MM part and hence presents an ideal model for solvents surrounding e.g. charge transfer reactions between reactant molecules, or reactant and surface. The QM part is the grid-based projector augmented wave code GPAW [2] which describes molecules and condensed matter systems on equal footing with a real-space grid representation of the wave-functions, and is massively parallelizable.
Example use includes QM water and MM water molecular dynamics simulations as well as QM metal surface and MM water interfaces. Although the present model has only been applied to simulations with the SCME water potential the model is general to any combination of QM and MM. The SCME does however require several parameters which need to be determined for any given material.
[1] K.T. Wikfeldt, E. R. Batista, F. D. Vila, and H. Joacute;nsson, PCCP 15(39), 16542 (2013)
[2] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, PRB 71, 035109 (2005)

Many important materials processes take place on time scales that vastly exceed the few microseconds accessible to molecular dynamics simulation. Typically, this long-time dynamical evolution is characterized by a succession of thermally activated infrequent events involving defects in the material. In the accelerated molecular dynamics (AMD) methodology, known characteristics of infrequent-event systems are exploited to make reactive events take place more frequently, in a dynamically correct way. The key feature is that the trajectory itself is allowed to find its own way out of each state, so that no prior assumptions need to be made about the transition states or available reaction paths. These methods have proven powerful for a range of processes, including metallic surface diffusion and growth, radiation damage annealing processes, and nanotube and nanowire dynamics. In this talk, I will discuss some recent advances that extend the power and range of applicability of the AMD methods. We recently formulated a local version of the hyperdynamics method that gives constant boost as the system size is increased, in contrast to standard hyperdynamics, for which the boost decays towards unity as the system size is increased. This opens the door for performing massively parallel hyperdynamics on arbitrarily large systems. While it was shown that this method gives accurate results for a homogeneous system, i.e., one in which all atoms are equivalent, test calculations indicate it is surprisingly accurate on inhomogeneous systems as well. I will discuss our recent advances in understanding the generality of this local hyperdymanics approach -- whether it is in fact valid, or can be made valid, for arbitrarily inhomogeneous systems. I will also discuss a novel approach for parallelizing the temperature accelerated dynamics (TAD) method, in which we speculatively accept every transition that is discovered at high temperature, spawning it as a new TAD run on its own processor, allowing it to proceed in parallel until/unless it is found to be an unneeded branch. In this approach, the correct low-temperature trajectory is generated as fast as the system can find the correct transitions at high temperature. Finally, time permitting, I will discuss our recent study of how close one can come to the ideal limit of hyperdynamics, when the bias potential is designed to "fill the basin" as completely as possible. This work has been supported in part by United States Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineering Division.

Although Molecular Dynamics (MD) is unrivaled at predicting the dynamical evolution of interacting atoms, the traditional method is often limited to time-scales that are much too short to capture critical phenomena of interest. Accelerated Molecular Dynamics (AMD) methods address this issue by accurately projecting the MD trajectory onto an accelerated state-to-state evolution. Of these methods, Temperature Accelerated Dynamics (TAD) hastens the projected evolution of a system by exploring possible transitions from each official state at an elevated temperature. To date, implementations of the TAD procedure have yet to leverage distributed computing to parallelize the exploration of multiple concurrent states. Here we present the Speculatively Parallel TAD (SpecTAD) method, which allows for the parallel scaling of serial TAD performance by executing a dynamically generated tree of speculative states on distinct computational cores. We introduce a few variations of the SpecTAD algorithm, and demonstrate its performance using specific materials applications.

Exascale computing presents a challenge for the scientific community as new algorithms must be developed to take full advantage of the new computing paradigm. Atomistic simulation methods that offer full fidelity to the underlying potential, i.e., molecular dynamics (MD) and parallel replica dynamics, fail to use the whole machine speedup, leaving a region in time and sample size space that is unattainable with current algorithms. We present an extension of the parallel replica dynamics algorithm [A. F. Voter, Phys. Rev. B 57, R13985
(1998)] by combining it with the synchronous sublattice approach of Shim and Amar [Y. Shim and J. G. Amar, Phys. Rev. B 71, 125432 (2005)], thereby exploiting event locality to improve the algorithm scalability. This algorithm is based on a domain decomposition in which events happen independently in different regions in the sample. We develop an analytical expression for the speedup given by this sublattice parallel replica dynamics algorithm and compare it with parallel MD and traditional parallel replica dynamics. We demonstrate how this algorithm, which introduces a slight additional approximation of event locality, enables the study of physical systems unreachable with traditional methodologies and promises to better utilize the resources of current high performance and future exascale computers.

In this invited talk I will discuss recent attempts at overcoming the timescale challenges of atomic-resolution simulations: (a) strain-boost hyperdynamics [Phys. Rev. B 82 (2010) 184114] for simulating primarily displacive events, (b) diffusive molecular dynamics (DMD) [PRB 84 (2011) 054103; PRB 86 (2012) 014115] for microstructural evolution driven by repetitive diffusion events and coupled displacive-diffusive processes. Ties to recent experimental observations of stress-driven glass crystallization [PNAS 110 (2013) 19725] and other diffusive-displacive transitions [Nature Mater. 13 (2014) 1007] will be discussed.

Classical molecular dynamics (MD) method is a powerful tool for simulating materials on the atomic scale, which can be used to compute a wide range of equilibrium and dynamical properties, such as elastic moduli, diffusivity, free energy barrier and reaction rate. However, classical MD neglects quantum mechanical nuclear zero-point energy and the nuclear tunneling effects. For systems containing light elements at low temperatures, quantum corrections are necessary in order to obtain a quantitative correct description of the system.
While solving the full time-dependent Shrodinger equation is usually limited to small systems, there are several approximate methods that have been developed for dealing with those complex, many-body systems. Among them, the ring-polymer molecular dynamics (RPMD) method is a popular semi-classical approach to include quantum nuclear effect in a classical MD simulation. RPMD is based on the path-integral formalism of statistical mechanics, which yields real-time MD trajectories that preserve the exact quantum Boltzmann distribution. It has been demonstrated in several systems that RPMD method does very well when used to calculated thermal reaction rates.
A long-standing problem is that conventional MD is limited to relatively short simulation timescales, which leaves many important phenomena out of reach. These so-called activated processes are typically characterized by long periods of uneventful vibrational dynamics punctuated by rapid transitions. While very powerful, RPMD suffers from an even more severe timescale problem as trajectories are generated not on a single system, but on a ring polymer where each bead is a complete replica of the system.
In classical MD simulation, accelerated molecular dynamics (AMD) methods address this challenge by concentrating on the sequence of transitions and shortening the inter-transition intervals. AMD methods require no apriori knowledge of the possible transitions, but to modify the dynamics so that transitions occur faster. In the present work, RPMD was coupled with Parallel Replica Dynamics (ParRep) method to accelerate the simulation of systems containing light atoms at low-temperatures. An example system of the diffusion of helium atom in bulk Fe, W, and Fe-Cr alloy at various temperatures will be discussed.

We performed parallel replica dynamics (PRD) simulations of cation diffusion in pyrochlores, a class of complex oxides that have been proposed for nuclear energy applications such as inert matrices for advanced nuclear fuel and nuclear waste encapsulation. In these applications, it is critical to understand the behavior of the material under irradiation, which notably induces disorder in the form of cation antisites. The investigation of cation transport in pyrochlores is challenging due to the two time scales involved: the slow diffusion of the cations prohibits the use of classical molecular dynamics (MD), while the concurrent fast diffusion on the oxygen sublattice proves hugely inefficient for techniques such as kinetic monte carlo (KMC) or temperature accelerated dynamics (TAD). In contrast, PRD allows one to exploit this separation of timescales in order to define states in terms of cation degrees of freedom alone. Using this approach, the effect of the chemistry on cation diffusion is investigated by looking at two compositions of pyrochlore (Gd₂Zr₂O₇ and Gd₂Ti₂O₇), known experimentally to have contrasting behaviors. The effect of the level of disorder is also discussed.

Dislocation nucleation from interfacial defect, such as grain boundary (GB), surface, dominates the mechanical response of nano-materials, which may have limited number of dislocation source inside. Molecular dynamics (MD) methods have been widely used to reveal atomistic details and energetics of the dislocation nucleation events. However, the regular MD simulation is limited to a timescale of nanoseconds and high strain-rates conditions (sim;10⁸ s^-1). In order to extend the timescale of MD simulation, some accelerated MD methods have been developed, such as hyperdynamics MD, bond-boost MD, strain-boost MD and adaptive-boost MD. Dislocation nucleation from grain boundary has been systematically studied using regular MD, but the nucleation process is unclear at a long time-scale. In the present research, we studied the dislocation nucleation process from a Σ9{221} symmetric tilt GB in Cu using the adaptive-boost MD method and an embedded-atom potential. The nucleation frequency, activation free energy, activation enthalpy, activation entropy and activation volume were given at a range of uniaxial tensile stresses and temperatures. The time acceleration factor is up to 10¹¹ (corresponding to a nucleation frequency of 1s^-1) at 300K, which is in agreement with the dislocation nucleation frequency in experiments. The activation entropy decreases with increasing tensile stress and is independent on temperature. Comparing to the activation energy at 0K given by nudged elastic band method, a dramatic decrease in activation free energy is observed with increasing temperature because of the entropic effects. The activation free energy firstly decreases smoothly with increasing tensile stress, where only one dislocation nucleates from the GB. Then, when the tensile stress approaches a critical stress, the regular MD can be used in this condition and the activation energy decreases rapidly induced by multi-dislocations collective nucleation from the GB. This result shows the different nucleation processes given by the accelerated MD with a long timescale and the regular MD with a nanoseconds timescale.

Fundamental insight into solidification in metals can significantly contribute to improve our understanding of their behavior during processing and under operating conditions. In the last decades, advances in computer simulation techniques made it possible to obtain relevant atomistic insight into the fundamental processes of solidification. Molecular dynamics and Monte Carlo simulations can be used to provide information concerning the dynamics of the interface during solidification, but the modeling of the initial nucleation remains challenging due to the extended timescales of the process associated to high free energy barriers. Nowadays, advanced computational methods like transition path sampling (TPS) have enabled the investigation of nucleation on the atomistic level. In this work, we employ TPS to investigate the nucleation during solidification in nickel, a material of high technological relevance. We initially focus on homogeneous nucleation in elemental nickel as a function of the undercooling. Here, a comparison of the temperature dependence of the free energy barriers to the predictions of classical nucleation theory (CNT) is discussed. As a second step towards more complex materials, we extend our study by including additional alloying elements and small Ni-clusters as seeds during heterogeneous nucleation. The transition state ensemble obtained from our TPS simulations provides atomistic insight into the structure and size of critical nuclei for different nucleation mechanisms (homogeneous and heterogeneous nucleation at defects), as well as nucleation barriers and rate constants. Furthermore, optimal candidates for reaction coordinates are identified based on local structural parameters and the nucleus shape using TPS and maximum likelihood analysis of the committor function. Such results provide relevant information to validate and improve existing thermodynamic models describing nucleation.

Chemical reactions between fast oxygen atoms and silicon dioxide reconstructed surfaces are important for vehicle re-entry into the atmosphere [1]. Therefore, there is interest in understanding the implications of Eley Rideal oxygen-silica reactions for mechanisms of energy transfer - especially surface-catalyzed recombination - at surfaces subjected to high-energy impacts. To enable simulations of oxygen-silica collisions, our goal is to design an empirical potential, based upon an integration of quantum-mechanical electronic structure data into a classical framework. In our scheme, we separate the potential energy into a bulk term and a defect term. The former is treated using a Tersoff-type empirical bond-order potential [2]. The latter consists of two- and three-body contributions and is described by the environment-dependent REBO method, which incorporates functional dependencies of the bond geometries and energetics on the numbers of each kind of nearest neighbor attached to each atom of the pair or trio [3]. The fitting strategy was designed to reproduce both the structure and energetics of surface defects of the types that were previously uncovered and described using first-principles calculations [4]. The chemical reactions between incoming high-energy oxygen atoms and defects at a defected α-quartz surface will be discussed in this talk alongside the details of the empirical interatomic potential scheme. The authors are grateful to Tom Schwartzentruber for helpful discussions. This work was supported in part by the Air Force Office of Scientific Research under grant no. FA9550-12-1-0486.
[1] I. Cozmuta, 39th AIAA Thermophys. Conf., AIAA-2007minus;4046 (2007).
[2] J. Tersoff, Phys. Rev. B 37, 6991 (1988); Phys. Rev. B 39, 5566 (1989).
[3] A. Kutana and K. P. Giapis, J. Chem. Phys. 128, 234706 (2008).
[4] R. Meana-Pañeda, Y. Paukku, K. Duanmu, P. Norman, T. E. Schwartzentruber, and D. G. Truhlar, J. Phys. Chem. C 119, 9287 (2015).

Large-scale atomistic simulations using classical molecular dynamics/molecular mechanics (MD/MM) methods enable modeling of materials, biological systems, and chemical reactions at surfaces/interfaces that are not easily accessible with experiments. The predictive capability and performance of MD/MC methods greatly rely on the empirical force field (EFF) used, which describes all interatomic interactions. The EFF parameters are generally fitted using least-squares local minima search algorithms to represent interatomic interactions and material properties obtained from quantum-mechanical based methods such as density functional theory (DFT), as well as available experimental data. Evolutionary algorithms such as genetic algorithm (GA) enable finding the global optimum of parameter space, even for EFF with very complex functional forms. However, systematic studies of the efficiency and accuracy of different optimization techniques for the determination of EFF parameters are still lacking. In the present work, we test the efficiency and effectiveness of a variety of global and local optimization schemes in the determination of EFF parameters. We first generated a training data set using DFT for Ir-O and Zr-N systems, which included pertinent structural and thermodynamics properties for a variety of structures. The EFF functional forms were chosen to be Morse for the Ir-O system and modified embedded-atom method (MEAM) for the Zr-N system. The optimization of the EFF parameters against the training data was carried out using (1) multi-start local optimization with common techniques including Simplex and Levenberg-Marquardt, and (2) single-objective GA. Using the random search as a base line, we compared the effectiveness and efficiency of the algorithms in terms of reaching the lowest error in representing DFT data, and number of function evaluations. We also investigated the use of multi-objective GA (MOGA) to remove the ambiguities in selection of weighting scheme in defining the fitness function. We analyzed accuracy and performance of MOGA in finding the Pareto-optimal solutions with respect to population size and number of objectives. Overall, this study provides a systematic approach for selecting a suitable optimization method to determine the optimal EFF parameters for use in classical simulations, which in turn enable more accurate prediction of material properties and chemical phenomena.

Organic monolayers covalently anchored to silicon surfaces represent a topic of fundamental and applied interest because such layers may be used as a thin dielectric, as a protection barrier or as a primer layers for use in microelectronics. The functionalization of silicon surface is also of interest in the development of chemical or biochemical sensors.
A key issue for such applications is the chemical and thermal stability of the monolayers. In previous works we investigated the chemical stability of hydrogenated and alkylated Si(111) surfaces (1, 2) towards O₂ and H₂O oxidizing species using Density Functional Theory.
In this work have we have employed reactive molecular dynamics simulations using the ReaxFF force field in order to understand the mechanisms of thermal decomposition of minus;CH₃, minus;CH₂CH₃, minus;CH₂CH₂CH₃, minus;CH₂(CH₂)₂CH₃, minus;CH₂(CH₂)₃CH₃ and minus;CH₂(CH₂)₈CH₃ grafted to the Si(111) surface. A maximum theoretical coverage of around 69% is predicted (3) for alkane monolayers on Si(111), implying that hydrogenated silicon atoms still remain on the surface. The only exception is the minus;CH₃ monolayer which has every atop Si atom is bound to a methyl group. Surface SiH groups play a key role in the stability of the monoalyers.
The thermal decomposition of the monolayers was investigated in the temperature range from 500 to 1500 K. The first step in the decomposition mechanism involves the breakage of surface SiH bonds and the diffusion of hydrogen atoms into the bulk, leaving surface silyl radicals on the surface. These radicals first abstract hydrogen atoms from the CH₂ group of beta carbon atoms which produces the desorption of alkene molecules. As the reaction proceeds, the surface coverage of organic molecules decreases and this allows the molecules to lie down on the surface as a consequence of the interaction with surface sylil radicals. This initiates the dehydrogenation of all the methylene groups. At the highest temperatures a full dehydrogenation occurs giving rise to the formation of silicon carbide on the surface.
The full coverage minus;CH₃ monolayer as an exceptional thermal stability due to the absence of SiH groups. For this monolayer, the dehydrogenation of the methyl group occurs by a H transfer to the Si atom bound to the C atom which has a much higher energy barrier than the transfer of an H atom from a methylene group to an adjacent sylil radical as is the case for the long chain monolayers.
For the most relevant elementary reaction steps identified in the MD simulations, Nudged Elastic Band calculations based on DFT energies were performed in order to obtain the energy profile along the reaction path from which activation energy barriers are calculated.
(1) F. A. Soria, E. M. Patrito, P. Paredes-Olivera. J. Phys. Chem. C, 2012, 116, 24607.
(2) F. A. Soria, P. Paredes-Olivera, E. M. Patrito J. Phys. Chem. C.2015, 119, 284.
(3) L. Scheres, B. Rijksen, M. Giesbers, H. Zuilhof Langmuir, 2011, 27, 972.

Metal thin films grown from the vapor phase on an insulating substrate is characterized by the nucleation and growth of three-dimensional islands and results in the formation of a polycrystalline thin film. The morphological evolution of the film microstructure is governed by kinetic limitations, which are a result of the dynamical interplay between the times-scales or different atomistic mechanisms, e.g., the adatom diffusion rate in different local environments and the vapor arrival rate. To understand how to influence the growth process, these mechanisms can be favorably simulated by the kinetic Monte Carlo (kMC) method. In order to construct a realistic simulation model of the growth process, detailed understanding of which atomistic mechanisms to include is required. However, existing literature on the mechanisms governing 3D nucleation and growth is somewhat divided. The thin film community argues that 3D island formation occurs by top-layer nucleation facilitated by limited downwards interlayer mass transport. On the other hand, the catalysis community advocates the process of upwards diffusion as the dominant mechanism of growth in the out-of-plane direction.
In this work, we attempt to bridge the divide between these communities by developing a physical model that encompasses aspects from both theories, within the framework of the kMC method and applied to FCC metals growing with a (111) surface parallel to the substrate. The goal of our model is to be able to replicate experimental results from both epitaxial as well as metal-on-insulator thin film growth cases. Here, selected results are presented for the temperature dependence of island shapes in Ag/Ag(111) homoepitaxy, as well as the annealing of metal islands on insulating substrates for different initial geometries.

The fracture toughness degradation of Fe based alloy used for nuclear pressure vessels is associated with the formation of Cu rich precipitates and matrix damages caused by neutron irradiation since these impurities become obstacles to moving dislocations. A Cu atom has a place in the Fe matrix as a substitutional solid solution, so it needs vacancies to freely migrate in the matrix. The vacancy mechanism is considered to be most probable mechanism for the precipitate formation in irradiated materials. The detailed atomistic process of Cu diffusion via vacancy mechanism has been investigated using molecular dynamics and kinetic Monte Carlo method in previous studies. However, their special or time scale is limited because of their own applicability. We have developed self-evolved kinetic Monte Carlo (SEAKMC) code [1] for Fe-Cu binary system with the EAM potential developed by Ackland et al. [2] to study the Cu diffusion via vacancy mechanism for a long time scale in the order of ms. SEAKMC is one of the on-the-fly kinetic Monte Carlo methods, which evaluates activation energies around defects in the simulation box at each time step using EAM potential by themselves. The activation energy variations caused by strain fields around defects or solute atoms are naturally considered in it. We applied the code to investigate the behavior of Cu-vacancy complex, diffusion coefficient of Cu and the mechanism of Cu-rich precipitate formation. The simulation result showed that the each life time of a Cu-vacancy complex is so short, but Cu- rich cluster is gradually formed over time.
[1] H. Xu, Y.N. Osetsky, R.E. Stoller, J. Phys.: Condens. Matter 24 (2012) 375402.
[2] G. J. Ackland, D. J. Bacon, A. F. Calder, and T. Harry, Philos. Mag. A 75 (1997) 713.

A novel process, the formation and growth of metallic Ag nanowires (NWs) and nanoparticles (NPs) on α-Ag₂WO₄ upon electron beam irradiation has been extensively investigated by different experimental techniques including transmission electron microscopy (TEM), field emission scanning electron microscopy (FE-SEM), energy dispersive spectroscopy (EDS) characterization, among others [1-3]. This is route widely employed to produce efficient photocatalysts, ozone sensors and bactericides.
The nucleation and formation of metallic Ag on α-Ag₂WO₄ initiate when Ag atoms diffuse from the interior material to the surface. This work uses DFT calculations and ab initio molecular dynamics simulations to investigate the geometrical and electronic structure of the most favorable surfaces: (100) and (001). This semiconductor exhibits a complex structure, which can be understood as an arrangement of AgO_x (x = 2, 4, 6, and 7) clusters used as constituent building blocks. Therefore, the relaxation process upon cleaving the surfaces from the bulk can induce significant rearrangements. Our calculations supply an atomistic approach to the local geometry and the electronic structure of the surfaces exposed to the electron beam irradiation. We also show that there are specific Ag atoms in the sub-surface positions that are prone to undergo diffusion processes with very low energy barrier (< 0.1 eV). Furthermore, our results point out that the injection of electrons decreases the activation barrier for this diffusion step [4]. This work combines experimental results and computer simulations to give theoretical insights in the NW and NC growth process.
References
[1] E. Longo, L. S. Cavalcante, D. P. Volanti, A. F. Gouveia, V. M. Longo, J. A. Varela, M. O. Orlandi and J. Andres, Sci.Rep., 3, 1676 (2013).
[2] J. Andrés, L. Gracia, P. Gonzalez-Navarrete, V. M. Longo, W. Avansi Jr., D. P. Volanti, M.M. Ferrer, P. S. Lemos, F. A. L. Porta, A. C. Hernandes and E. Longo, Sci. Rep. 5, 5391 (2014).
[3] E. Longo, D. P. Volanti, V. M. Longo, L. Gracia, I. C. Nogueira, M. A. P. Almeida, A. N. Pinheiro, M. M. Ferrer, L. S. Cavalcante and J. Andrés, J. Phys. Chem. C, 118, 1229 (2014).
[4] W. da Silva Pereira, J. Andrés, L. Gracia, M. A. San-Miguel, E. Z. da Silva, E. Longo and V. M. Longo, Phys. Chem. Chem .Phys. 17, 5352 (2015).

Many chemically complex phenomena are currently still beyond the reach of first principles molecular dynamics techniques. This is not an accuracy issue: the problem mostly arises because the necessary model system sizes are too large, and/or the required simulation times are too long. In many situations using classical simulation techniques is not a viable alternative, as suitably general and accurate “reactive” force fields are not available, nor are fitting databases a priori guaranteed to contain the information necessary to describe all the chemical processes encountered along the dynamics. Standard “QM/MM” techniques combining quantum and classical zones in a single calculation are also not devoid of difficulties, especially for applications on processes involving sustained mass transport into and out of the (e.g., fast moving) QM zone.
In this talk I will argue that this situation makes it necessary to use molecular dynamics techniques capable of incorporating accurate QM information generated at run time during the simulations [1]. It also creates a novel market for databases coupled with specially-tuned Machine Learning force fields which minimise the computational workload by allowing QM subroutine calls only when “chemically novel” configurations are encountered along the system&’s trajectory. I will present one such “Learn On the Fly” scheme, effectively unifying First-Principles Molecular Dynamics and Machine Learning into a single, information efficient simulation scheme capable of learning/predicting atomic forces through Bayesian inference [2]. Strategies for dealing with large, dynamically evolving QM zones in simulations running on high-end parallel platforms will also be discussed [3,4].
[1] A. Gleizer, G. Peralta, J. R. Kermode, A. De Vita and D. Sherman, Phys. Rev. Lett., 112, 115501 (2014)
[2] Z. Li, J. R. Kermode and A. De Vita, Phys. Rev. Lett., 114, 096405 (2015)
[3] Cf., e.g., the US-DOE INCITE on SiO2 ML-Fracture Project https://www.alcf.anl.gov/projects/sio2-fracture-chemomechanics-machine-learning-hybrid-qmmm-scheme
[4] M. Caccin, Z. Li, J. R. Kermode and A. De Vita, Int. J. of Quantum Chemistry 2015, DOI: 10.1002/qua.24952

The lattice dynamics (LD) method forms the foundation for our understanding of the properties of crystalline systems. In it, an assumption of small lattice vibrations permits the use of a harmonic approximation, in which the energy is expanded to second order in the atom displacements. The resulting Hamiltonian can be solved in closed form for the dynamic and thermodynamic behaviors, using either quantum or classical mechanics. Volume-dependent thermal effects can be evaluated with the quasi-harmonic extension of this treatment.
At low temperature or high pressure, LD works very well. However, atomic vibrations grow with increasing temperature and decreasing density, and the harmonic approximation begins to fail. Consequently, there are many conditions of interest for which LD is inadequate. In such situations, the most reliable alternative is molecular simulation.
Molecular simulation as normally practiced for the evaluation of anharmonic thermal properties does not exploit the harmonic character of the crystal to improve calculation of averages. Thus, even at conditions where LD provides an excellent description, simulation essentially “starts from scratch” in evaluating the properties, making no use of the LD characterization. Consequently, the precision of the properties computed by molecular simulation is severely compromised, inasmuch as the stochastic averaging must contend with fluctuations contributed by the harmonic component of the intermolecular potential. There is a clear inefficiency in computing stochastically a large contribution that is already known analytically.
We present a set of new methods to compute properties of crystalline phases by molecular simulation. The methods are highly efficient to the extent that the modeled system is harmonic, but at the same time they do not rely on the harmonic character to produce correct results. Rather they yield data that focus strongly on characterizing the anharmonic contributions directly, so that precise results can be obtained without noise introduced from sampling harmonic effects. The methods are applied to several model systems, including the Lennard-Jones model, an extended Finnis-Sinclair model of bcc tungsten, clathrate hydrates modeled with the TIP4P potential, and fcc aluminum modeled using density functional theory. The approaches are shown to provide data of precision equal to that given by conventional methods, while requiring many orders of magnitude less computation time.

We have extended the Nested Sampling algorithm[1][2][3][4] to simulate materials under periodic boundary conditions and constant pressure. We show how this methodology can, in turn, be used to determine phase diagrams of materials directly from the potential energy of the configurations generated during the Nested Sampling. This process is both efficient and can be done in a highly automated fashion. Thus, the only inputs required to generate a phase diagram are the composition of the material and the desired pressure and temperature ranges. In particular, solid-solid phase transitions are recovered without any a priori knowledge about the structure of the solid phases. We have applied this technique to the Lennard-Jones system, aluminium, and the NiTi shape memory alloy.
References
[1] J. Skilling, AIP Conf. Proc. 735, 395-405 (2004)
[2] J. Skilling, Bayesian Anal. 1, 833-859 (2006)
[3] L.B. Pártay, A.P. Bartoacute;k and G. Csányi, J. Phys. Chem. B 11, 10502-10512 (2010)
[4] R.J.N. Baldock, L.B. Pártay, A.P. Bartoacute;k, M.C. Payne and G. Csányi, submitted (2015)

In systems containing atoms with different charges and electronegativities, charge transfer and the electrostatic interactions between the resulting charged atoms are of critical importance. In molecular dynamics simulations charge equilibration (QEq) approach based on the electronegativity equalization principle allows interatomic charge transfer and assigns new charge on atoms based on their surrounding environment. Charge equilibration can be solved self consistently or with an extended Lagrangian scheme. We introduce an advanced extended Lagrangian charge equilibration solver that is based on the fast inertial relaxation engine (FIRE) algorithm (qeq/fire). It is a damped dynamics method with additional velocity modifications and adaptive time steps. Results show that charge equilibration via the qeq/fire algorithm converges significantly faster than the standard extended Lagrangian scheme.
Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy National Nuclear Security Administration under contract DE-AC04-94AL8500

Despite the availability of large multi-processor computer architectures, physically relevant timescales are often out of reach in molecular dynamics (MD) simulations, due to the sequential nature of this approach. We have developed a framework that parallelizes an MD trajectory by dividing it up in the time domain where separate processes are started from predicted accessible states. The problem of parallelization is then rephrased as a problem of prediction. Being able to predict a sufficient set of basins of trajectory directly translates into being able to splice the long-time simulation in the time domain. Using machine learning tools that were built for accurate predictions and decision-making on related topics, we show how certain long-time simulations of rare-event dynamics can be efficiently parallelized even for systems that have a heterogeneous distribution of transition rates between basins.

“...everything that living things do can be understood in terms of the jigglings and wigglings of atoms.”
-Richard P. Feyman
With the advent of supercomputers molecular dynamics (MD) has become routine to investigate physical phenomena, with reasonable accuracy. However, it has two major shortcomings namely, size and time limitation. Lately, efficient parallelization of computers has become a major tool to overcome the size limitation obstacle. During the past decade many attempts have been made to create new algorithms to speed up atomistic simulations. One of such new methods is Diffusive Molecular Dynamics (DMD), to achieve longer time scale this algorithm treats displacive (jiggling) and diffusive (wiggling) mechanisms differently. Here we introduce a new standalone open source MD/DMD package, MIT Atomistic Parallel Package (MAPP). While equipped with most of popular atomistic potentials, MAPP has several parallelization schemes, minimization schemes, and time integration schemes with variable time steps. Several physical examples will be discussed and compared to experiments and other popular modeling methods.

A key application of molecular dynamics (MD) simulations of crystal defects concerns the calculation of defect-phonon couplings, namely the thermal drag force that gives rise to the viscous dynamical law used in all coarse grained crystal plasticity simulations. Rather than simulate long defect trajectories in large boxes, we present a method to directly extract the defect force autocorrelation from much smaller boxes, and relate this through Zwanzig's projection technique to the frictional force acting on crystal defects. The defect force, position and velocity is extracted in a novel manner, by defining a projection operator for a crystal defect parametrised from static calculations[1]. We demonstrate the efficiency our approach for small mobile defects in Iron and Tungsten, extracting accurate diffusion constants from the defect force autocorrelation over very short simulation times far smaller than those required to evaluate the diffusivity from trajectory analysis.
We have also used the projection operator to derive a clear analytic form for the defect-phonon interaction, resolving a long standing failure of phonon scattering theory to predict the anomalous temperature independent drag force seen in many MD simulations of nanoscale defects. We show that an anomalous temperature independent coupling arises as vibrations orthogonal to defect motion produce a defect force even to linear order in the vibrational amplitude, which is expressly forbidden in phonon scattering calculations.
[1] TD Swinburne, SL Dudarev, AP Sutton, Phys. Rev. Lett. 113, 215501. http://dx.doi.org/10.1103/PhysRevLett.113.215501

The design of next-generation functional materials and the quest to gain insights into dynamical processes at the atomic scale entail efficient and accurate atomistic simulations. The success of these simulations hinges on the accuracy, robustness, and transferability of the force fields (FFs) employed to describe interatomic interactions. One aspect of developing these FFs is to determine the set of variable parameters in their functional forms that closely reproduce available experimental or quantum mechanical data. In general, accurate FFs possess complex functional forms with numerous parameters, and call for sophisticated fitting strategies. Genetic algorithms (GA) coupled with local optimization such as the Simplex method provide an efficient way to scan the parameter space for such an optimization problem. Here, using technologically relevant Co-C and Zr-N binary systems as examples, we demonstrate the applicability of GA for force field fitting using single and multiple objectives for two different functional forms: (a) Tersoff-like bond order (Co-C), and (b) Modified embedded atom method (Zr-N). These newly developed FFs accurately describe the structural, elastic, thermodynamic and defect properties of the bulk condensed phases in excellent agreement with density functional theory calculations and experiments. We employ these FFs to (a) investigate self-assembly of Co₅₀₀ clusters and C₆₀ fullerenes into porous metal-organic hybrid structures and the morphology of obtained hetero-structures using large-scale molecular dynamics, and (b) atomic-scale pathways of defect migration in cubic ZrN under radiation conditions. These representative examples are used to highlight the effectiveness of evolutionary strategies in bridging the atomistic and electronic length scales.

A new interatomic potential model that can simultaneously describe metallic, covalent and ionic bonding nature has beed developed by combining the second nearest-neighbor modified embedded-atom method (2NN MEAM) and charge equilibration (Qeq) method. Paying a special attention to the removal of known problems found in the original Qeq model, a mathematical form for the atomic energy is newly developed and carefully selected computational techniques are adapted for energy minimization, summation of Coulomb interaction and charge representation. The model is applied to the Ti-O and Si-O binary systems selected as representative oxide systems for a metallic and a covalent element. The reliability of the present 2NN MEAM+Qeq potential is evaluated by calculating fundamental physical properties of a wide range of titanium and silicon oxides, and comparing with experimental data, DFT calculations and other calculations based on (semi-) empirical potential models.

In chemistry, condensed matter physics and materials science it is often the case that only the initial and final states of a transition are known. The task is to use computer simulations to determine the minimum energy path (MEP) which shows the mechanism of the transition (i.e. the evolution of the system during the transition) and gives an estimate of the activation energy. Knowing the mechanism can be helpful in tuning the system to favor or inhibit a given transition. The nudged elastic band (NEB) method is often used for relaxing an initial path towards an MEP, usually the closest to the initial path. This can be sufficient for simple systems. However, for more complex systems, a large number of MEPs can be present between the given initial and final states and the task is to find which MEP involves the lowest increase in energy and thereby the smallest activation energy.
We first discuss the efficiency of different versions of the NEB method. The improved tangent NEB (IT-NEB) [1] has been shown to exhibit better convergence than the original NEB (O-NEB) [2]. However, for relaxing long paths with a limited number of images, we show that the O-NEB is more efficient than the IT-NEB. We then introduce a modified NEB method which retains the advantages of both the O-NEB and IT-NEB and is efficient for relaxing both short and long reaction paths.
We then present a procedure for carrying out global optimization of reaction paths. The procedure is based on dividing up the path by generating intermediate configurations from a heredity transformation which consists in mixing other configurations. Paths between consecutive configurations are relaxed toward MEPs using the modified NEB method.
The nucleation of a dislocation typically involves displacements of thousands of atoms, and
finding the optimal MEP is a challenging task. Ge films deposited on a Si(001) substrate are used in many applications in electronics and photonics. The lattice constant difference between Ge and Si induces a large strain in the film and leads to the formation of dislocations. Understanding the mechanisms for the formation of the dislocations can help grow high quality films. By carrying out a global optimization, we examine and review various mechanisms for the formation of dislocations in the Ge/Si(001) system [3].
[1] G. Henkelman and H. Joacute;nsson, J. Chem. Phys. 113, 9978 (2000)
[2] H. Joacute;nsson, G. Mills, and K. W. Jacobsen, Nudged Elastic Band Method for Finding Mini-
mum Energy Paths of Transitions, in Classical and Quantum Dynamics in Condensed Phase
Simulations, edited by B. J. Berne, G. Ciccotti, and D. F. Coker (World Scientific,1998)
[3] E. Maras, A Stukowski, T. Ala-Nissila and H. Joacute;nsson, (submitted 2015).

We will report on developments in adaptive kinetic Monte Carlo (KMC) approaches for modeling problems in materials science using empirical potentials. We will present a study of the influence of the prefactor in the Arrhenius equation for the long time scale motion of defects, simulated by the adaptive KMC technique. We will show that calculated prefactors vary widely between different defects, thus it is important to determine them accurately within the simulations. Even calculating the prefactor to one significant figure accuracy made a great impact during the simulations. The results were verified by reproducing many events using a combination of Molecular Dynamics and Temperature-Accelerated Dynamics simulations.
We will also present an approach to building superbasins on the fly in an adaptive KMC framework to massively accelerate simulations in the presence of low barriers. This will be shown in the context of the use of automatic defect identification and classification, using graph theory to both classify defects and graph isomorphisms to maximize the reuse of barriers. We will illustrate how this has been utilized in both surface growth problems and the modelling of the evolution of radiation damage to extend timescales to be experimentally relevant.
Examples of the growth of both ZnO and CdTe surfaces will be used to illustrate the power of these methods in identifying very complicated growth processes that are critical in these materials. The modeling of radiation damage in Fe, FeCr and NiCr systems will be used to illustrate the ability of these methods to lengthen timescales to be industrially relevant and to model the effect of multiple radiation damage events at realistic rates in these systems.

The microstructure of metal films formed via deposition processes has significant effects on their performance and longevity for applications in microelectronics and as conductive coatings. The ability to predict, and ultimately control, metal film microstructure is of significant importance and requires a fundamental understanding of the deposition process used for its formation. In order to address this, a three-dimensional on-lattice kinetic Monte Carlo (KMC) method was developed for atomistic simulations of single crystal and polycrystalline metal electrodeposition. The method uses the semi-empirical multi-body embedded-atom method (EAM) potential and a multiple lattice approach. The resulting computational method, KMC-EAM, enables highly descriptive simulations of metal deposition processes that capture both atomistic and experimentally relevant scales.
As a sample application of the KMC-EAM method, copper electrodeposition is studied under two different sets of conditions: (i) kinetically controlled electrodeposition onto single crystal under galvanostatic direct-current conditions and (ii) polycrystalline copper under potentiostatic direct-current conditions. Four types of surface processes are considered during the electrodeposition process: deposition, dissolution, surface diffusion and grain boundary diffusion. Simulation results were found to agree well with findings from experimental studies that the evolution of the root-mean-squared roughness of the deposit during the early stages of deposition follows a power law relationship t^0.5. These results support the use of the KMC-EAM method to bridge the atomistic-to-lab scale gap for simulation of metal deposition processes.

Using a combination of molecular dynamics (MD) and adaptive kinetic Monte Carlo (KMC) it is possible to investigate dose effects in radiation cascade studies. When a large energy is imparted to an atom such as from an energetic neutron within a material or by accelerating an ion beam and injection across a surface, the initial motion can be well represented by classical MD. The system then evolves more slowly and by using an adaptive KMC approach longer time scales can be accessed.
Examples are given using this hybrid approach to model sputtering from gold surfaces and cascades generated in bcc Fe. It is shown that for low energy sputtering of metal surfaces, a large amount of recrystallisation occurs between successive particle impacts. It is further shown how certain preferred facets can form if the surface is not low index.
A second example shows how He bubbles can aggregate in bcc Fe and how the size of the bubbles becomes limited by the large energy barriers that the He atoms have to overcome when the strain field around the bubble becomes too large. The predicted bubble sizes are in close agreement with experiment.

Kinetic Monte Carlo (KMC) methods have been a successful technique for accelerating the time scales and system sizes possible in atomistic simulations. However, a requirement for its success is a priori knowledge of all relevant reaction pathways. This can be difficult for many systems with complex chemistry, such as high-temperature shock-compressed materials or heat shields undergoing phenolic pyrolysis, which can consist of hundreds of molecular species and thousands of distinct reactions. In this work, we develop a method for statistically learning a reduced set of elementary reactions and corresponding rate coefficients from molecular dynamics (MD) simulations that can be used to build an accurate and efficient KMC model.
We consider a MD simulation of shock-compressed liquid methane using the ReaxFF potential, and use bond length criterion to identify almost 200 species and more than 1000 elementary reactions over the course of 216 picoseconds. Rate coefficients are then estimated for each elementary reaction based on their observed frequency. We show that a KMC simulation on this set of reactions and rate coefficients reproduces the molecular concentrations of the dominant species in the MD simulation to acceptable levels of accuracy. The rate coefficients can then be further refined by solving a convex relaxation of the maximum likelihood estimation problem on a Langevin model for chemical reaction networks. We then add an L1-regularization term to reduce the number of reactions in the model while maintaining an acceptable level of accuracy in reproducing the MD.
We find that we can reduce the original chemical reaction network by more than 75% of reactions while successfully reproducing key features of the MD simulation. This method has the potential to make the construction of complex KMC models a routine and computationally efficient task independent of chemical intuition. The reduced model also ensures better performance and decreases the number of expensive quantum chemical computations needed to refine rate coefficients using transition state theory.

Time evolution of an atomistic model can be often represented as a sequence of discrete random incremental changes in the system's configuration. We have found an efficient solution to the notorious computational problem of kinetic trapping in which simulated random trajectories repeatedly visit the same uninteresting configurations in a low energy basin without exploring the rest of the configuration space. Nature's brute-force solution to kinetic trapping is to wait for the trajectory to eventually escape the trap after revisiting it over and over. Termed Kinetic Path Sampling (KPS), our method modifies the trajectory to avoid revisiting any trapping configurations even once and to predict where the system will escape the trap. Then, based on this knowledge, KPS reconstructs an entire stochastic trajectory ending in the predicted escape. Unlike other existing methods, KPS is exact and correctly predicts where and when the trajectory will escape without any prior knowledge of the trap. Being essentially immune to kinetic trapping, the method's accuracy and efficiency is demonstrated in simulations of anomalous diffusion on a random energy landscape and of the kinetics of phase separation in a lattice model of FeCu binary alloy [1].
[1] M. Athenes and V. V. Bulatov, Path Factorization Approach to Stochastic Simulations, Phys. Rev. Lett. 113, 230601 (2014).

Direct simulation with Molecular Dynamics (MD) is an extremely powerful tool to shed light on the atomistic behaviors that control many of the properties of materials. The predictive power of MD however comes at a steep computational price, limiting the system sizes and simulation times that can be achieved in practice. While the size limitation can be efficiently addressed with massively parallel implementations of MD based on spatial decomposition strategies, the same approach usually cannot extend the timescales much beyond microseconds. I will discuss an alternative, parallel-in-time, strategy - the Parallel Replica Dynamics (ParRep) method - that aims at addressing the timescale limitation of MD for systems that evolve through rare state-to-state transitions. While ParRep was introduced by Voter more than 16 years ago, it has significantly evolved over the last few years. I will review the theoretical foundations of the method and discuss recent results that show that ParRep can provide arbitrarily accurate results for any definition of the states, while providing a parallel speedup that can reach up to the number of replicas used. I will then demonstrate the practical usefulness of ParRep by presenting different examples of materials simulations where access to long timescales was essential to study the physical regime of interest.

Temperature-accelerated dynamics (TAD) [1] is a powerful ‘unbiased&’ method to accelerate the dynamical simulations of systems with infrequent events. In this method, a basin-constrained high-temperature molecular dynamics (MD) simulation is used to determine the system evolution at low temperature. In particular, in some cases TAD can accelerate molecular dynamics simulations by factors as large as 10⁶. As a result, for small systems TAD can be used to carry out realistic simulations of non-equilibrium processes such as thin-film growth over experimental time-scales. However, since the computational work associated with TAD scales approximately as N³ (where N is the number of atoms in the system) this places a severe limit on the size of systems which can be simulated. One solution to this problem is the use of parallel TAD simulations [2] based on spatial decomposition combined with our synchronous sublattice (SL) algorithm [3]. Using this method, we have demonstrated [2] that TAD simulations can be carried out over extended time-scales and for very large system sizes. However, while this significantly improves the scaling, the size of activated events is limited by the decomposition size. In addition, one needs to carefully control the cycle (communication) time between processors to avoid errors in the dynamics. Accordingly, it is desirable to develop additional methods to improve the scaling of serial TAD. Here we discuss recent work which we have carried out in order to characterize the relevant bottlenecks and also improve the scaling behavior. These include the use of a local saddle-point method to determine the activation barrier of high-temperature events [4] as well as a more recently developed local method for identifying [5] high-temperature transitions. By combining these methods with a previously developed method [6] to block previously observed transitions during the basin-constrained high-temperature MD we have found a significant improvement in the scaling of serial TAD. Possible methods for additional improvement will also be discussed.
[1]. M. R. Soslash;rensen and A. F. Voter, J. Chem. Phys. 112, 9599 (2000).
[2]. Y. Shim, J. G. Amar, B. P. Uberuaga, and A. F. Voter, Phys. Rev. B 76, 205439 (2007).
[3]. Y. Shim and J.G. Amar, Phys. Rev. B 71, 125432 (2005).
[4]. Y. Shim, N.B. Callahan, and J. G. Amar, J. Chem. Phys. 138, 094101 (2013).
[5]. Y. Shim and J.G. Amar, unpublished.
[6]. G. Subramanian and A.F. Voter, unpublished.
* Supported by NSF DMR-1410840

Some recent developments in saddle point finding and long time scale dynamics methods will be presented, including biased gradient squared optimization; the kappa-dimer method; solid state phase transitions; and basin-constrained molecular dynamics saddle-search based adaptive kinetic Monte Carlo.

Adaptive kinetic Monte Carlo (aKMC) is a powerful atomistic modeling tool, which has shown early success at capturing the time-evolution of materials under ion and neutron irradiation. Results and developments pertaining to its application to Ni and Ni-based alloys are presented. Modifications to the kinetic Activation Relaxation Technique (k-ART), an extension of aKMC, are introduced to facilitate the study of systems reaching 256000 atoms. The aging of these fcc materials between successive displacement cascades in a given volume is captured using this implementation of k-ART. Furthermore, k-ART simulations in smaller cells are performed to shine light on the mechanisms controlling diffusion of individual defect clusters. Finally, we discuss the challenges and limitations of adaptive KMC in dealing with concentrated alloys, which possess very rough potential energy surfaces. The work was supported as part of the “Energy Dissipation to Defect Evolution”, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science.

Lattice Kinetic Monte Carlo (KMC) simulations offer a powerful alternative to using ordinary differential equations for the simulation of complex chemical reaction networks. Lattice KMC provides the ability to simulate and analyze the local spatial behavior of chemicals in the reaction network, resulting in a more detailed description of the reaction pathway. Within the Lattice KMC framework, the transition probabilities of different reactions can differ by orders of magnitude, creating subsets of processes that occur more frequently or more rarely. Consequently, processes that have a high probability of occurring may be selected repeatedly without actually progressing the reaction network (i.e. the forward and reverse process for the same reaction). In order to avoid long simulation run times resulting from the repeated occurrence of such processes, it is necessary to throttle the transition probabilities in such a way that avoids altering the overall selectivity. Likewise, as the reaction progresses, new frequently occurring species and reactions may be introduced, making dynamic throttling of the transition probabilities a necessity. We present a dynamic steady-state detection scheme with the goal of accurately throttling rate constants in order to optimize the KMC run time without compromising the selectivity of the reaction network. The algorithm has been applied to a complex chemical reaction network, specifically that of methanol oxidative dehydrogenation, as well as additional pathways on CeO₂(111) ultimately leading to formaldehyde, CO, methanol, CO₂, H₂ and H₂O as gas products.

When complimented with accurate techniques for calculation of atomistic diffusion rates, the kinetic Monte Carlo (KMC) method can serve as an important tool for the simulation of temporal and spatial evolution of surface phenomena such as epitaxial growth, adatom-island diffusion, coarsening, and morphological transformations. To enhance its predictive capacity, we have developed a self-learning kinetic Monte Carlo (SLKMC) ¹ method, in which the standard KMC is combined with automatic generation of a table of microscopic events, facilitated by a pattern recognition scheme. Each time the system encounters a new configuration, the algorithm initiates a procedure for saddle point search. Nontrivial paths are thus selected and the fully characterized transition path is permanently recorded in a database for future usage. Once the data base of all possible single and multiple atom processes is built, the system evolves automatically and efficiently by picking diffusion mechanisms of its choice. I will present application of the method to the diffusion and coalescence of 2-dimensional Cu, Ni and Ag adatom and vacancy clusters on the (111) surface of these three transition metal, covering both homo- and hetro-systems. Of interest are multiple atom processes revealed in the simulation whose presence may have been ignored otherwise. For adatom clusters varying in size from 2 to 100, I will discuss the dependence of the diffusion coefficient, effective energy barriers, and dominant mechanism (periphery atom or concerted-cluster motion), on cluster size. I will highlight the role played by specific diffusion processes and show that a crossover from collective island motion to periphery diffusion takes place at critical sizes which are specific to the metallic system in question. Comparison will be made of process responsible for the diffusion of same-sized clusters in homo and hetero systems. Results will be compared with experiments, where available, and with those from KMC simulations based on a fixed catalogue of diffusion processes. I will also provide details of a recent extension of the techniques to three dimensions. Also, for the case of early stages of sub monolayer island coarsening, I will point to the kinetic stabilization of certain island sizes resulting from specifics of adatom detachment /attachment processes².
^*Work done in collaboration with O. Trushin, A. Kara, A. Karim, G. Nandipati and S. I. Shah and supported by US-NSF.
¹ S. I. Shah, G. Nandipati, A. Kara, and T. S. Rahman, Phys. Rev. B 88, 035414 (2013); A. Karim, A. et al., Phys. Rev. B 73, 165411 (2006).
² G. Nandipati, S.I. Shah, A. Kara, and T. S. Rahman, Phys. Rev. B 88, 115402 (2013).