Symposium DS02—Advanced Atomistic Algorithms in Materials Science

Usually when one thinks of calculating the properties of a material or molecule, the first calculation that comes to mind is electronic structure, which can be approximated on traditional computers using techniques such as density functional theory. Appropriately, most materials-related quantum algorithms developed so far have targeted electronic degrees of freedom. However, there are other types of materials simulation that would require the use of a quantum computer for accurate results. These include problems for which one is simulating vibrational, excitonic, or nuclear degrees of freedom. Though these types of simulations are relatively unexplored in the quantum computation community, they are relevant to many industrial applications. Here, we present recent work on designing quantum algorithms for important industrial applications such as microscopic energy transfer, spectroscopy, thermal properties, and analytical chemistry. We begin with an overview of general quantum computing methods for materials and chemistry simulation. Next, we discuss the design of efficient low-level quantum data structures for these applications. Metrics related to computational complexity are compared between electronic structure and these new application areas, leading us to postulate that simulating these alternative degrees of freedom will be possible before simulating electronic structure. Finally, we discuss the mathematical ways in which the problem goals are different for each application, with an eye toward guiding quantum algorithm design in these areas of materials simulation.

First-principles electronic structure calculations are the state-of-the-art for solid-state physics, chemistry, and materials science; however, their predictive capability is constrained by the limitations on simulated length- and time-scale due to computational cost and the unfavorable scaling of the algorithms. 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 the PAOFLOW framework for the generation of local basis representations for effective ab-initio tight-binding schemes. PAOFLOW is a software tool to efficiently post-process standard first principles electronic structure plane-wave pseudopotential calculations in order to promptly compute, from interpolated band structures and density of states, several quantities that provide insight on transport, optical, magnetic and topological properties such as anomalous and spin Hall conductivity, magnetic circular dichroism, spin circular dichroism, and topological invariants. After reviewing recent advances in the establishment of theoretical framework, including the generation of efficient atomic basis sets, I will illustrate the method with a few recent examples from our research on 2D ferroelectric materials, including the demonstration of a transistor prototype based on electrically switchable persistent spin helix and a memory device that exploits the domains of switchable spin textures and spin-to-charge conversion.

Fcc structured Fe-Ni alloys containing 36 at.% Ni are known as Invar alloys. They are characterized by an anomalous low, practically zero, thermal expansion at room temperature and ambient pressure, a finding that awarded C.E Guillaume the Nobel Prize in Physics already 1920. It is clear that the origin of the Invar effect is linked to the magnetic properties of the material but a complete explanation for the effect is not yet agreed on. Phenomena such as Fe-moment spin flips, high-to-low-spin transitions, non-collinear magnetism, local chemical environments, and strong electron correlations have been considered as possible origins, which together possess a challenge for the state-of-the-art first-principles methods for disordered magnetic alloys. [1,2]

The Invar effect is known to disappear when the Ni composition is even slightly increased.[2] However, it was observed that the Invar effect could be induced by high pressure also for Fe-Ni alloys with high Ni concentrations.[3]

In this study we investigate the origin of the pressure induced Invar effect in the disordered fcc Fe_1-xNi_x alloy for x = 0.36, 0.50, and 0.75 from first principles calculations. We perform a large pool of different non-collinear supercell electronic structure calculations, and let the magnetic state and atomic structure gradually relax towards their most preferred states. Using this approach we observe a transition from ferromagnetic states to noncollinear or almost collinear spin-flipped states as the volume decreases for all three compositions. The volumes at which the transitions occur correspond to the volumes observed to display volume-pressure anomalies in previous experiments. Our study confirms that the pressure induced Invar effect is a result of a magnetic transition from a ferromagnetic state to a complex magnetic state in Fe-Ni alloys.

[1] van Schilfgaarde, M., Abrikosov, I. & Johansson, B. Origin of the Invar effect in iron–nickel alloys. Nature 400, 46–49 (1999)

[2] E.F. Wassermann. The Invar problem. Journal of Magnetism and Magnetic Materials, 100(1):346-362, 1991

[3] Leonid Dubrovinsky, Natalia Dubrovinskaia, Igor A. Abrikosov, Marie Vennström, Frank Westman, Stefan Carlson, Mark van Schilfgaarde, and Börje Johansson. Pressure-induced Invar effect in Fe-Ni alloys. Phys. Rev. Lett., 86:4851-4854, 2001

Machine Learning (ML) is a promising approach to improve the accuracy of exchange-correlation (XC) functionals for Density Functional Theory (DFT). Two key developments are necessary to successfully leverage ML for XC functionals. First, because of the limited information provided by semi-local ingredients like the density and its gradient, computationally efficient, physically relevant, and nonlocal descriptors of the density or density matrix must be designed as input to the model. Building on our recently developed CIDER functional for the exchange energy, we explore possible improvements to the model as well as its computationally efficient implementation. We also investigate the application of ML to efficient, orbital-dependent ingredients like Rung 3.5 functionals [1]. Second, robust ML architectures and workflows must be developed to train functionals that are both accurate and transferable. We discuss how Gaussian Processes yield smooth and transferable models and how different regularization and training schemes impact the performance of parametric functionals.
[1] B. G. Janesko, J. Chem. Phys. 133, 104103 (2010); https://doi.org/10.1063/1.3475563

Electron dynamics in external electric fields governs the behavior of solid-state electronic devices. First-principles calculations enable precise predictions of charge transport in low electric fields. However, studies of high-field electron dynamics remain elusive due to a lack of accurate and broadly applicable methods. Here we develop an efficient approach to solve the real-time Boltzmann transport equation with both the electric field term and ab initio electron-phonon collisions. These nanosecond-long simulations provide field-dependent electronic distributions with a femtosecond resolution, allowing us to investigate both transient and steady-state transport in electric fields ranging from low to high (>10 kV/cm). The broad capabilities of our approach are shown by computing nonequilibrium electron occupations and velocity-field curves in Si, GaAs, and graphene, obtaining results in quantitative agreement with experiment. Our approach sheds light on microscopic details of transport in high electric fields, including dominant scattering mechanisms and valley occupation dynamics. Our results demonstrate quantitatively accurate calculations of electron dynamics in low-to-high electric fields. This method, implemented in our open-source code Perturbo [1], provides a timely and robust tool to advance the discovery and design of novel electronic materials.

[1] J.-J. Zhou, J. Park, I-Te Lu, I. Maliyov, X. Tong, M. Bernardi. Comput. Phys. Commun. 2021, 264, 107970.

We introduce electrochemical dynamics with implicit degrees of freedom (EChemDID), a model to describe electrochemical driving force in reactive molecular dynamics simulations. The method describes the equilibration of external electrochemical potentials (voltage) within metallic structures and their effect on the self-consistent partial atomic charges used in reactive molecular dynamics. An additional variable assigned to each atom denotes the local potential in its vicinity and we use fictitious, but computationally convenient, dynamics to describe its equilibration within connected metallic structures on-the-fly during the molecular dynamics simulation. This local electrostatic potential is used to dynamically modify the atomic electronegativities used to compute partial atomic changes via charge equilibration. Validation tests show that the method provides an accurate description of the electric fields generated by the applied voltage and the driving force for electrochemical reactions. We demonstrate EChemDID via simulations of the operation of electrochemical metallization and valence change cells. The simulations predict the switching of the device between a high-resistance to a low-resistance state as a conductive filament is formed and resistive currents that can be compared with experimental measurements. In addition to applications in nanoelectronics, EChemDID could be useful to model electrochemical energy conversion devices.

Magnetic Nanoparticles (MNPs) can organize into novel structures in solutions with excellent order and geometry that are important for many industrial and biomedical applications. Studies of self-assembly of smaller MNPs is challenging due to the complicated interplay between VDW, electrostatic, dipolar, steric, hydrodynamic, and external magnetic field interactions. An all-atom molecular dynamics (AMD) simulation technique can decipher such complex interactions and can account for the effect of the MNP core, surfactant layers, solvent molecules. Here we present a novel AMD simulation method to study the dynamics, self assembly and properties of MNPs structures with various particles size and shapes, surfactant and solvents. We demonstrate the use and effectiveness of the model by simulating the self-assembly of several commonly obtained equilibrium structures like ring, line, chain, clusters by oleic acid surfactants coated spherical and cubic Magnetite (Fe₃O₄) nanoparticles in vacuum and in explicit solvents and under the presence of a uniform external magnetic field. We found that the long-range electrostatic interactions can favor chain formation over rings, the VDW and steric interactions by the surfactants promote MNPs cluster growth and the interplay of solvent and increased surface energy can reduce the rotational diffusion of the MNPs. The algorithm has been parallelized to take advantage of multiple processors of a modern computer and can be used as a plugin for the popular LAMMPS simulation software to study the behavior of small magnetic nanoparticles and gain insights into the physics and chemistry of different magnetic assembly processes with atomistic details.

REACTER is a heuristic method for modeling chemical reactions in classical molecular dynamics simulations, implemented in LAMMPS as fix bond/react. The authors recently extended LAMMPS to support alphanumeric labels for atom types, bond types etc., which enables the pre- and post-reaction templates required by the REACTER protocol to be portable between different simulations and greatly simplifies the task of creating simulation-ready reaction templates. To further increase the generality of reaction templates, support for wildcard characters within atom types has been added, along with the automatic assignment of interaction types for new bonds, angles, etc. based on the involved atom types. In some cases, this feature can express a class of reactions with one pair of reaction templates, where previously dozens may have been required. Advanced reaction constraints have also been added, including an Arrhenius constraint to enforce an effective activation energy, a root-mean-square-deviation option for complex geometrical constraints, as well as a custom constraint that leverages LAMMPS’ powerful built-in variable framework. Other new features include variable support for various inputs (e.g. to allow reaction rates or cutoffs to be dependent on overall conversion), on-the-fly update of molecule IDs, and the ability to create new atoms positioned with respect to the reaction site. The new features are applied to modeling polymeric, thermosetting and composite materials, and advanced applications of the new reaction constraints are demonstrated. For example, REACTER is shown to accurately reproduce mechanically-induced bond breaking, as characterized by third-order DFT-based tight-binding (DFTB3) simulations, via a constraint on the total potential energy of the involved atoms.

Grain structures impact Cadmium Telluride/ Cadmium Sulfide (CdTe/CdS) solar cells since grain boundaries can act as recombination centers for carriers. Computer simulations such as molecular dynamics (MD) can be a convenient and cost-effective method for investigating the evolution of grain structures during the growth of materials. Recently, MD simulations have been successfully applied to study the growth of polycrystalline CdTe/CdS structures, revealing complex zincblende and wurtzite grains along with highly disordered atoms that are often observed in experiments. However, little efforts have been made to quantify the grain structures observed in MD simulations in literature. In this work, the grain tracking algorithm originally developed by Panzarino et al. has been extended for polycrystalline CdTe/CdS structures. Specifically, we focus on the face-centered-cubic (fcc) sublattice of the zincblende crystal structure. The common neighbor analysis and centrosymmetric parameter values are used to calculate the orientation of each atom in the grain tracking algorithm. We then apply this modified approach to analyze the CdTe/CdS films obtained from our MD polycrystalline growth simulations. We demonstrate that our method can provide a variety of useful information such as grain domains, grain orientations, plane indices, and sample texture. Moreover, kinetic analysis of microstructure evolution of the materials can be performed. These results will add to the understanding of granular microstructure that can describe physical, chemical, mechanical, and electronic properties.

Atomistic modeling of metallic systems needs the inclusion of electron-ion interactions in order to have realistic dynamics. For example, it is generally known that around 20% of the projectile's kinetic energy could be absorbed by the electrons and redistributed within the metallic system outside the collision cascade. Capturing this in classical molecular dynamics is rather challenging as these simulations usually do not include explicit electrons. We present a novel model that captures both the electron-ion interactions as well as energy transport within the electronic system. Our Langevin based stochastic model is able to catpure the non-adiabatic collective dynamics of phonons. Moreover, we have demonstrated that we can also apply the same model to the high energy events, such as radiation damage. In this work, we present an important aspect of electronic heat transport - we have developed an atom based diffusion model where energy is transferred directly between atoms, taking into account convective transport of electronic energy. Our model reproduces the macroscopic heat diffusion behaviour, but is also able to incorporate environment dependent and hydrodynamics effects. Furthermore, the model is coupled to the Langevin based electron-ion model, thus creating a unified approach for enabling classical molecular dynamics with electronic effects. We will present the theoretical aspects of the model as it is implemented in an extension for LAMMPS code. Finally, we present an application of compression shockwave dynamics through a system consisting of a metal and an insulator. This example demonstrates the strength of our approach, as electronic heat is transported both by diffusion and convection.

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

The local structure of conjugated polymers is vital for many optoelectronic properties, such as charge conduction and charge photogeneration at donor-acceptor interfaces. We perform molecular dynamics simulation in Gromacs to study the local arrangement of these polymers. Conjugated polymers have aromatic rings as repeat units. Due to the stiffness and large size of these polymers, atomistic simulations are computationally expensive. Hence, to increase simulation efficiency, we exploit the stiffness of aromatic moieties using virtual sites to represent multiple atoms within each ring. These virtual sites do not have mass but contain electronic properties; hence, we significantly reduce degrees of freedom while retaining the shapes and dipole moments of the aromatic rings. As a result, these simulations are ~10 times faster compared to atomistic simulations. We simulate poly(3-hexylthiophene-2,5-diyl) (P3HT) and rhodanine-benzothiadiazole-coupled indacenodithiophene (O-IDTBR), as well as mixtures of these two materials. Our simulations predict a persistence length of 3.8 nm, a density of 0.91 kg/m³ for P3HT, and an amorphous scattering profile that compare well to experimental results. In addition, the crystal structure of IDTBR is stable when arranged in experimentally-determined unit cell parameters. Our simulations allows us to examine equilibrated molecular configurations at the interface of P3HT and IDTBR, and reveal preferences in terms of contacts between specific moieties.

Modeling glasses requires extended timescales—since their simulation involves the slow cooling of a liquid into an out-of-equilibrium disordered solid. As such, reliable, yet computationally efficient empirical interatomic forcefields are key to facilitate the modeling of glasses. Although analytical forcefields offer optimal computational efficiency, their parametrization is challenging as the high number of parameters renders traditional optimization methods inefficient or subject to bias. Here, we present a new parameterization method based on machine learning, which combines ab initio molecular dynamics simulations, Gaussian Process Regression, and Bayesian optimization. By taking the examples of silicate and chalcogenide glasses, we show that our method yields new interatomic forcefields that offer an unprecedented agreement with ab initio simulations, at a fraction of computational cost. This method offers a new route to efficiently parametrize new interatomic forcefields for disordered solids in a non-biased fashion.

The last few years have seen considerable advances in the development of machine-learned interatomic potentials. A very important aspect of the parameterization of transferable potentials is the generation of training sets that are sufficiently diverse, yet compact enough to be affordably characterized with high-fidelity reference methods. We formulate the generation of a training set as an optimization problem where the figure of merit is the entropy of the distribution of atom-wise descriptors. This can be used to create a fictitious potential to explicitly drive the generation of new configurations that maximally improve the diversity of the training set. I will show how this strategy can provide an automated and scalable solution to generate large training sets without human intervention.

In recent years, machine learning has been used to improve the training and simulation efficiency of “bottom-up” coarse grained potentials, which guarantee the preservation of thermodynamic quantities in the coarse grained space [1]. While these potentials can be used to capture the physics of systems on large length and time scales that are traditionally inaccessible with molecular dynamics, they are known to generally have a high degree of complexity.

It has been shown that Gaussian process regression can be used as a framework for active learning of force fields from ab initio data [2], and we extend this idea to coarse grained models. In particular, we study the utility of active learning over passive learning schemes in the efficient modeling of coarse grained potentials. Moreover, we examine how Gaussian processes can be used as a tool to better understand the importance of model complexity in coarse grained problems.

[1] W. G. Noid, Jhih-Wei Chu, Gary S. Ayton, Vinod Krishna, Sergei Izvekov, Gregory A. Voth, Avisek Das and Hans C. Andersen. “The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models,” J. Chem. Phys. 128, 244114 (2008).
[2] 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,” NPJ Computational Materials 6 (2020).

Many Materials Science problems of interest involve a number of atomic species and this often makes deriving potentials for these systems very challenging. In this work we will present simulations on dopants and interfaces in systems based on ZnO that have been modelled using a ReaxFF potential approach, using machine learning methods to fit the potentials. This is an important transparent conducting oxide system and the doping of it and the interfaces of it and other materials are often critical in the mechanical performance of devices that utilise this. We will also present work on using a machine learning potential to study high entropy ceramic materials and compare this with ab-initio calculations on these systems.

Machine-learned interatomic potentials are far more expressive than traditional physically motivated interatomic potentials like Lennard-Jones, Stillinger-Weber, Embedded Atom Potentials, etc. While they can be far more accurate, they are also more likely to be completely wrong outside of the training domain, are more difficult to train reliably, and are computationally expensive. We have developed MLIPs for the Hf-Ni-Ti shape memory alloy. We share cautionary tales, best practices for generating training sets, and demonstrate how community tools make for "easy entry" to realistic thermodynamic modeling with these potentials.

We introduce an unsupervised machine learning workflow to autonomously parameterize Tersoff Bond Order Potential using two-dimensional (2D) materials as a representative example. This is a 14-dimensional space and even with the state-of-the-art techniques, the decision of fixing the search boundaries still remains a daunting task. This has been tackled by continually sampling the potential energy surface using the Latin Hypercube Sampling (LHS) technique starting with an arbitrary search bounds. The high dimensional parameter space is projected into a low-dimensional space using Principal Component Analysis (PCA) and are tagged with corresponding energy values estimated from molecular dynamic simulations using LAMMPS. We next employ clustering techniques to identify the range of parameter sets that fall in the desired energy space. A simplex optimization is performed next to identify the parameter set with energies close to the actual energies and guide the range selection for the next iteration. We demonstrate that the search can autonomously achieve target minima on the potential energy surface - a bond order model developed for Carbon 2D polymorphs outperforms existing parameterizations. Our workflow represents an elegant method to overcome issues with selection of parameter ranges and weights for multiobjective high-dimensional optimization of force-field parameters.

Graph neural networks are powerful tools for representing a wide variety of subjects. In recent years these tools have become state of the art for the social network domains and elsewhere. However, they are relatively new in the field of materials science. Prior approaches utilizing graphs such as CGCNN or MolCLR do not utilize the breadth of feature tools available in social networks and these could offer advantages when representing compounds. In this work, we apply graph attention networks (a variant of graph neural networks) on a graph representation of inorganic compounds including unique global, edge, and node features to create representation vectors for the whole graph. Emphasis is placed on learning relations of bonds rather than nearest neighbors by relying on chemical domain knowledge such as dipoles and electronegativity. Node features are extracted through clustering algorithms and neighbor-based encoding. Global features will include bonding motifs, structural topology characteristics, composition-based feature vector, clustering coefficients, and thermodynamic descriptors. Models are constructed to predict material properties and the attention map attribute of the graph attention network is used to learn and visualize the importance of the information obtained from neighboring nodes.

In the development of new machine learning potentials, it is essential to guarantee the correct description of the properties of interest. Our goal is to correctly describe long time scale properties of fusion materials to simulate the surface degradation caused by the contact with the plasma. In order for that to be accomplished, it is necessary to estimate the errors made for those properties by candidate potentials. In this work we developed a method to estimate confidence intervals for the calculation of energy barriers using the SNAP potential.

Collective variables (CVs), a low dimensional representation of atomic structures, are the cornerstone of many enhanced sampling techniques that enable efficient free energy landscape explorations. Conventionally, CVs are simple geometric descriptors, such as atomic bonds and angles. However, it often needs intuition and trial-and-error testing to determine which CV can lead to accurate free energy calculations and efficient sampling.

In this work, we demonstrate a data-driven method to learn CV with a multitask learning framework. Short molecular dynamics simulations around the transition states and basins are used to train a multitask network consisting of an encoder that encodes atomic structures to a low dimensional latent space and two separate downstream parts that map the latent space to predictions of basin class labels and potential energies. The trained latent space can be used as an effective CV for umbrella sampling bias and free energy calculations. This learning framework is demonstrated on model systems, small molecules, and metal surfaces.

[1] L Sun et al., Multitask machine learning of collective variables for enhanced sampling of rare events, arXiv preprint arXiv:2012.03909

Fundamental understanding of the nanotube-protein interactions at the interface is a crucial issue for the development and rational design of the various biocompatible and hybrid bionanomaterials for applications in bio-sensing, biocatalysis, and drug delivery. We investigated the interactions of bovine serum albumin (BSA) with multi-walled carbon nanotubes (MWCNT) by all-atom molecular dynamics (MD) simulation. MD simulations indicated that hydrophobic, van der Waals (VDW), electrostatic, and hydrogen bonding interactions contribute to the total binding, and these components have been separated out in our work. As a cutoff distance of 5A^o used to see the bio-nano surface interactions observed between CNT and BSA surface for 80 ns MD simulation. Analysis of binding components of the total interaction energy (-63 Kcal/mol) shows that the driving force of the binding is mainly hydrophobic (-30 Kcal/mol) interactions by the non-polar side chain residues (valine314, proline302 & 303, ala311). CNT interacts with positively charged residues (Arg and lysine312) negatively charged residues (Glu 299 & 310, Asp 311) in proteins through the cation-π and anionic-π interactions, respectively. Notably, the BSA molecular structure and function preserved largely after its adsorption on the CNT surface. By investigation of the protein adsorption on the CNT surface would be helpful to better understand the stability and bio-nanosurface interactions at a molecular level.

Keywords: protein corona, carbon nanotubes molecular dynamics, hybrid bionanomaterials

Directly computing mass transport coefficients in stochastic models requires integrating over time the equilibrium correlations between atomic displacements. Here, we show how to accelerate the computations by decoupling short- and long-time correlations, which statistically amounts to conditioning the mass transport estimator through a law of total diffusion. We illustrate the approach with kinetic path sampling simulations of atomic diffusion in alloys where percolating solute clusters trap the mediating vacancy. There, Green functions serve to generate first-passage paths escaping the traps and to propagate the long-time dynamics. When they also serve to estimate mean-squared displacements via conditioning, colossal reductions of statistical variance are achieved.

Conjugated polymers (CPs) possessing the characterization of relatively rigid conjugation backbone and peripheral flexible side chain have attracted considerable attention toward the application of organic electronic and optoelectronic devices. The computational prediction of the thermomechanical behavior of CPs to serve the needs of materials design and prediction of their functional performance is a grand challenge due to the prohibitive computational times of all-atomistic (AA) molecular dynamic (MD) simulations. Atomistically informed coarse-grained (CG) modeling is an essential strategy for making progress on this problem. In this work, focusing on the widespread explored regioregular poly(3-hexylthiophene) (P3HT), we develop a temperature-transferable CG model based upon the energy-renormalization (ER) approach that allows for modeling of CPs dynamics over a wide temperature range and accessing greater spatiotemporal scales. The results show that the CG model faithfully and accurately capturing the short- and long-timescale dynamics, i.e., Debye-Waller factor and segmental relaxation time. Systematically exploration of the mechanical properties, such as tensile and shear, reveals that the cohesive interaction exhibits a great impact on the mechanical response. Our work highlights the necessity and feasibility of ER when modeling a conjugated polymer that possesses branched architecture, and sheds new light on the rational design of the temperature-transferable CG model of complex polymers for their thermomechanical and conformational prediction.

We tailored and calibrated a sequential molecular dynamics (MD)/time-stamped force-bias Monte Carlo (tfMC) simulation methodology to investigate the growth mode, the resulting crystalline structure, and orientation of Cu thin films obtained by sputter deposition on epitaxial TiN substrates. Our studies reveal for the first time that at the very early stage of growth, BCC-Cu grows pseudomorphically on the TiN(001) substrate as a very thin continuous film with the BCC-Cu[001]//TiN[001] growth direction. As the Cu thin film thickness increases the film transforms, through the Nishiyama-Wasserman mechanism, from BCC into predominantly FCC-Cu with abundant nanotwins, which is the same type of structure obtained in the experiment conducted here via a dc magnetron sputter deposition technique to grow Cu on TiN001 at 105 °C [1]. The MD/tfMC simulations also reveal that on the N-terminated TiN(111), Cu shows a very poor wettability and the the FCC-Cu(111) film grows vertically in the form of tall 3D islands. On Ti-terminated TiN(111) surface, however, FCC-Cu(111) initially grows in the form of 2D islands with high wettability. With additional Cu deposition, a triangular misfit dislocation network is generated at the Cu(111)//Ti terminated TiN(111) interface with subsequent formation of a two-layer nanotwin with its twinning plane parallel to the surface substrate. The present simulation results provide deep insights into some complexities of the growth mechanisms operating during Cu growth on TiN and shed light on related thin film growth in weakly-interacting metal/substrate systems.

[1] X.M. Zhang, S. Shao, A.S.M. Miraz, C.D. Wick, B.R. Ramachandran, W.J. Meng, Low temperature growth of Cu thin films on TiN(001) templates: Structure and energetics, Materialia, 12, 100748 (2020).

* Work funded in part by the NSF EPSCoR program, awards OIA-1541079 and OIA-1946231.

The discovery of correlated insulating states at half filling at a twist angle of 1.08⁰ known as t the magic angle has stirred tremendous interest in twisted bilayer graphene (TBG). Such exotic physics and electronic properties are intrinsically dependent on the atomic structure of the material. In TBG, a periodic variation of different stacking regions emerges from the relative rotation between the graphene layers. At large twist angles (> 1.47⁰), the positions of the atoms are well-described by considering a rigid rotation between the individual graphene layers. However, the symmetry changes significantly for low twist angles ( < 1.47⁰) due to atomistic relaxations that expand the stable AB stacking region while shrinking the size of the AA and SP stacking regions. While this atomic relaxation is very well-described in the literature, the effect of the symmetry change on the charge density of TBG remains unexplored since electronic structure calculations at small twist angles require large, and often computationally infeasible, supercells. Therefore, we develop a computationally efficient framework that enables the exploration of the charge density symmetry of low twist angle TBG. Our approach is based on the evolution of the high intensity Bragg peaks in the Fourier representation of the charge density, obtained from density functional theory, as the twist angle decreases. The structural relaxation is found to correspond to low intensity satellite peaks in the Fourier representation. By accounting for these satellite peaks, it is possible to observe the transformation of symmetry in the charge density distribution from high to low twist angle TBG. One striking outcome of our framework is the observation of electron localization in the AA region of low twist angle TBG. The framework shows clearly the effects of atomistic relaxations on the symmetry of the charge density of TBG at small twist angles, which may be useful for better understanding of the exotic physics of TBG.

Coarse-graining (CG) approaches reduce the total number of particles in molecular dynamics (MD) simulations and often enable the use of longer integration time steps. Both of these effects serve to reduce computational expense and therefore provide access to larger length and time scales for a given system. The Martini model is a widely used CG force field, flexible and transferrable enough to find applications ranging from biological to polymer-based systems, including reactive CG MD simulations. However, CG models like Martini typically combine several atoms into a representative 'bead,' meaning that individual reaction sites effectively disappear as atomic identities are absorbed into the Martini beads. In this work we use this loss of specificity as a tool that enables the simulation of different chemical reaction mechanisms that may arise from the same CG model. This approach therefore isolates reaction mechanism as the variable in these simulations.

As proof of concept, we consider monomers that share a similar molecular mechanical description but polymerize by different mechanisms: radical chain-growth polymerization and radical step-growth polymerization. We represent these physically similar monomers with the same CG model, a 'Universal Monomer,' and use reactive CG MD simulations to prepare highly crosslinked thermosets of each monomer. Since the monomers are represented by precisely the same model, differences observed in the polymerization reactions or in the mechanical properties of the resulting thermosets are a function of the reaction mechanism. We compare on-the-fly simulated polymerization, network topology, and simulated mechanical testing of the two reactive CG MD systems and discuss related experimental results. Further, we suggest that this Universal Monomer approach is a reasonable survey of 'mechanism-space' with application to the rational design of chemical feedstocks for photopolymerization-based additive manufacturing techniques.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344, IM release number LLNL-ABS-823466.

Identifying the aluminum oxide/aluminum (Al₂O₃/Al) interfacial structure is important to improve the performance of the wide range of applications that employ this interface, such as to mitigate sources of noise in superconducting electronic devices and to improve the chemical and mechanical properties of corrosion coatings and structural composites. It will also help gain deeper understanding of oxide propagation into the metal in this fundamental system. However, the structure of the interface has not been resolved at the atomic scale, and the stages of growth of the oxide into the metal are unknown. Experimental studies claim that the interface is sharp, coherent and indicate the presence of undercoordinated aluminum atoms, but provide no information on the atomic arrangement at the interface. Computational studies have primarily only looked into the stability of different terminations of the oxide on the metal, computed the work of adhesion, electronic structure, etc. Given that the interface presents an environment vastly different from either the bulk oxide or metal, it is likely that it has different structural features and properties. In this study, we utilize ab initio Grand Canonical Monte Carlo, a physically motivated, bias-free method, to accurately search the interface configuration space and understand the structure and evolution of the interface. We show that 12.5% aluminum vacancies are present at the interface at the oxide/metal equilibrium. As the oxygen chemical potential is gradually increased, the interface propagation appears to occur in a layer-by-layer fashion. Oxygen is first incorporated at the FCC hollow sites at the interface, and then at the tetrahedral interstitial sites. This resembles oxygen incorporation at the Al(111) surface, after which the stacking sequence of Al₂O₃ is regained. We also discuss the trends in bond lengths, coordination environment and electronic structures of these interfacial configurations. The results of this study have implications for interfacial defect transport, charge trapping and adhesion.

FCC metal thin films produced by atom-by-atom deposition methods can be formed with dense arrays of parallel coherent twin boundaries (CTB’s) having average spacings in the nanometer regime. There has been great interest in these nanotwinned structures due to their expected unusual properties. To access those properties, it must be possible to control the average CTB spacing. Yet demonstrations of experimental control remain elusive and existing models (molecular dynamics simulations, phase field and analytical models) provide conflicting predictions. We have simulated nanotwin formation as a function of deposition rate and temperature on (111) surfaces in FCC metal thin films using a rejection-free kinetic Monte Carlo (KMC) approach. In our simulations, atoms are added to a 2-D hexagonal lattice and moved according to an energy landscape that accounts for the differences in transitions between twin and FCC sites using an adaptation of the variable step size method. Rates and temperatures that mimic experimental conditions are used. Layer growth occurs by nucleation, growth, and island interaction. However, interactions with simulation cell boundaries strongly influence results and a choice between cell sizes small enough to run many times for good statistics and large enough to obtain realistic behavior limited predictive capabilities. We find that for experimentally relevant rates, temperatures, and grain sizes, nanotwin formation takes place in the nucleation-dominated regime and developed an analytical model based on the KMC methods to describe this. Predictions for nanotwin formation are described.

Accelerated molecular dynamics methods have the capability to increase the time scales probed in rare-event simulations by orders of magnitude, with no loss in accuracy. However, these methods have not been widely applied because of several drawbacks. In hyperdynamics, one of these drawbacks is the appropriate choice of the boost potential. If the boost is too high, the simulation can be inaccurate. If the boost is too low, the method is ineffective. Obtaining an optimal boost requires detailed knowledge of the underlying free-energy surface, which is computationally demanding. However, with advances in parallel computing, this problem is becoming tractable. In this work, we introduce an adaptive hyperdynamics method for global hyperdynamics that can predict optimal parameters for the bond-boost potential. By quantifying the probability distribution for a rare event to occur within a time t as a function of the boost, we observe a stair-like property. We can correlate stair steps with features of the boosted potential surface in a way that allows us to avoid inaccuracies associated with too large of a boost. We have implemented this method using a Python wrapper that allows us to test boosts for many replicas of a system in parallel, so a near-optimal boost can be identified relatively rapidly. The method is currently applicable to rare events governed by a potential-energy barrier, but extensions are possible in the future. We demonstrate the method in simulations of metal nanostructures in solvent with adsorbed capping molecules, a problem relevant to the shape-selective growth of metal nanocrystals.

Building a chemical kinetics model traditionally requires conducting many experiments over the course of several years. Recently developed atomistic simulations had made it possible to develop kinetics models within months, using the large amount of data generated by these simulations. These kinetics models allow scientists to predict the chemical evolution of complex systems that can involve several thousands of reactions and species, while reducing computation time by up to four orders of magnitude compared against atomistic simulations. When building these kinetics models, there are usually two frameworks depending on the size of the species that are being studied. First, a kinetics model usually used for gas phase systems: this type of kinetics model has been thoroughly studied and can involve several thousands of molecules and reactions. The reactions are described using molecules, but this description prevents such a kinetics model from studying larger molecules, because these kinetics models fail to predict the creation and the reactivity of unobserved molecules. The second type is used for long chains or solids, like in polymer growth or soot formation. These models are usually simplified to include only a few reactions occurring at some reacting sites. This description is not adapted to smaller molecules, and the reduced number of reactions cannot always reflect the complexity that these processes have. These two frameworks have shown successes throughout the years, but there is a lack of kinetics models that would span across all these different systems, starting from small molecules and going all the way to solid particles.

In this work, we propose a new description for reactions in kinetics models that allows us to study small particles and large clusters: the atomic-level features. In our proposed approach, the reactions are described by the atoms surrounding the reactive site. This description resonates with the approach used for long chains and solids, but it is also applied to small molecules, resulting in a model that can describe both sizes. To demonstrate the efficiency of our approach, we study carbon particle growth in the atmosphere of icy giant planets. In this application, the initial composition is liquid methane which then grows into long hydrocarbons chains. The mechanism here is not fully understood, and a kinetics model that could describe accurately the three different states that this system goes through is needed. We demonstrate that the atomic description of reactions outperforms the molecular description on the prediction of the larger chain growth, while maintaining similar accuracy on the small molecules. This is achieved because the atomic features allow the creation of molecules unobserved in the training set, in addition to accurately predicting their reactivity. Furthermore, the atomic-level model requires fewer unique reactions than a conventional molecular-level model. This induces a more accurate estimation of the reaction rates and the ability to parametrize the kinetics model with less data. A corollary of this ability is that it is possible to extrapolate in time and in size using the atomic-level features. Finally, the atomic kinetics model can be transferred to different starting molecules while maintaining a high level of accuracy. Therefore, we find that our proposed atomic featurization of reactions enables simultaneous description of small molecules and large clusters.

Thermal protection systems (TPS) in hypersonic vehicles are primarily made from carbon-based materials, due to the high heat capacity and high energy of vaporization of carbon, which protects the vehicles from damage in reentry conditions. However, at very high temperatures, reactions with oxygen cause the erosion of the TPS through formation of CO and CO₂. In this work, we model defect generation in multi-layer graphene due to adsorption and diffusion of oxygen using an atomic scale multi-lattice KMC approach. Simulations based on this KMC model reveal details of the pitting process driven by atomic oxygen in multi-layer graphene for specific temperatures and partial pressures of atomic and molecular oxygen. The model is developed using separate, interacting lattices for the carbon and oxygen sites. The carbon lattice includes all possible carbon sites in a multilayer graphene structure representative of an HOPG material. In every time step of the simulation, a change in the state of a site in the oxygen lattice affects the possible states of the neighboring sites in the carbon lattice. The use of separate lattices makes it possible to effectively distinguish the oxygen and carbon sites and allow them to communicate throughout the simulation as various chemical interactions are tracked. The model implemented here consists of 20 surface reactions and 23 reactions on graphene edge sites, for a total of 43 reactions. Surface groups such as lactone-ethers and ether-lactone-ethers are formed on the surface as a result of reactions between neighboring epoxies and are shown to have a low energy barriers to CO and CO₂ formation. These groups are thus crucial to the initial formation of defects as they lead to formation of CO and CO₂ from initially defect free surfaces. For computational efficiency, only active oxygen sites around the previously executed site in the last timestep are checked for available events and are updated to the existing list of available events. We use the approach to show the effects of high symmetry grain boundaries and terraces on the formation and growth of defects in the oxidation process, and we observe pit formation enhanced by the presence of the exposed edges. Using the multi-lattice KMC model it is possible to simulate much longer times than in direct MD simulation. Scales of milliseconds are easily reached on a laptop computer without significant effort to optimize efficiency. Thus, the approach enables analyses of relatively large domains over experimental laboratory time scales.

Self-interstitial atoms (SIAs) are almost exclusively observed in irradiated materials. In part because of their presence, unique materials science phenomena occurs during irradiation of materials, such as radiation-induced segregation and radiation-induced (dis)ordering. While SIAs in pure metals typically diffuse very rapidly, within the timescale accessible by molecular dynamics, SIA kinetics in disordered alloys can be inhibited by configurational traps.
In this talk, I will present two examples. First, I will present the case of SIA diffusion in Ni(x)Fe(1-x) concentrated solid solution alloys. Second, I will present the case of SIA-based reordering of irradiated Ni3Al alloys.
In the case of the Ni(x)Fe(1-x) alloys, we use a mix of micro-second scale molecular dynamics, the kinetic Activation Relaxation Technique and accelerated lattice kinetic Monte Carlo to explain why the SIA diffusion coefficient is a non-monotonic function of alloy chemical composition.
In the case of Ni3Al, we use molecular dynamics and temperature-accelerated Markov models with
Bayesian estimation of rates (TAMMBER ), which show that even though the SIA diffusion coefficient increases as it re-orders the system, the re-ordering kinetics, counterintuitively, slow down. We demonstrate that a relatively simple mean-field rate theory is able to quantitatively predict this behavior.

Atomistic/Continuum coupling methods have now been developed for over 25 years. However, methods in which non-elastic behavior is occurring in the continuum domain remain computationally challenging due to the need for full 3d solutions in the continuum domain using finite elements. Greens function methods have the potential to greatly reduce the computational costs and enable large multiscale simulations. Here, progress toward a GFM method is described. First, the Flexible Boundary Condition Method of Sinclair et al., which is a GFM for an atomistic region in an infinite elastic space, does not give unique solutions even for linear problems. Second, efficient/usable GFM methods require a transition from a full lattice GF to a continuum GF, and this is shown to create numerical errors. However, these errors can be controlled and reduced. Third, starting from a fully-atomistic description, a coarse-graining of the outer (continuum) boundary is presented that greatly reduces the computational costs with essentially no loss of accuracy. These results show that an emerging GFM method will solve the bottleneck of current multiscale methods.

In this work, we propose a new general method for analyzing and calculating diffusivity and ionic conductivity in media with strong ionic correlations from molecular trajectories [1]. Previously, the only two options were the dilute uncorrelated approximation, which converges rapidly but it is inaccurate, and the exact Green-Kubo method that is prohibitively expensive for complex and large systems. The herein presented approach automatically extracts and utilizes the collective diffusion eigenmodes of the displacement correlation matrix to denoise the calculation of the transport properties. The proposed approach is universally applicable, simple, and provably superior to previously available methods, exhibiting speed ups of several orders of magnitude. Thus, it opens wide opportunities to study correlated diffusion in previously inaccessible complex electrolytes.

[1] N Molinari, Y Xie, I Leifer, A Marcolongo, M Kornbluth, and B Kozinsky, "Spectral denoising for accelerated analysis of correlated ionic transport", Phys. Rev. Lett.

A method is presented for time independent density functional calculations of excited electronic states using direct orbital optimization [1,2]. The calculations are variational and converge on stationary states, so they provide estimates of atomic forces in the excited state. Calculations of materials can be calculated using either a real space grid or a plane wave basis set. The functional can be orbital density dependent (ODD) as in self-interaction corrected functionals or be of the Kohn-Sham (KS) form as in the generalized gradient approximation. The implementation for KS functionals involves two nested loops: (1) An inner loop for finding a stationary point in a subspace spanned by the occupied and a few virtual orbitals corresponding to the excited state; (2) an outer loop for minimizing the energy in a tangential direction in the space of the orbitals. For ODD functionals, a third loop is used to find the unitary transformation that minimizes the energy functional among occupied orbitals only. Combined with the maximum overlap method, the algorithm converges in challenging cases where conventional self-consistent field algorithms tend to fail. The benchmark tests presented include charge-transfer excitations, defect states in diamond and excitons in silica.

[1] 'Variational Density Functional Calculations of Excited States via Direct Optimization', G. Levi, A. Ivanov and H. Jónsson, J. Chem. Theory Comput. 16, 6968 (2020).

[2] 'Direct Optimization Method for Variational Excited-State Density Functional Calculations Using Real Space Grid or Plane Waves', A. V. Ivanov, G. Levi, E. Ö. Jónsson and H. Jónsson, J. Chem. Theory Comput. (in press, 2021). Manuscript available on arXiv.

The large number of grid points per atom required for accurate real-space Kohn-Sham Density Functional Theory (DFT) calculations restricts their efficiency. In this work, we present an approach to accelerate such calculations severalfold, without loss of accuracy, by systematically reducing the cost of the key computational step: the determination of the Kohn-Sham orbitals spanning the occupied subspace. This is achieved by systematically reducing the dimension of the discrete eigenproblem that must be solved, through projection into a highly efficient discrete discontinuous basis. In calculations of quasi-1D, quasi-2D, and bulk metallic systems, we find that accurate energies and forces are obtained with 8–25 basis functions per atom, reducing the dimension of full-matrix eigenproblems by 1-3 orders of magnitude.

In this talk, the speaker will present SPARC: Simulation Package for Ab-initio Real-space Calculations. SPARC can perform Kohn-Sham density functional theory calculations for isolated systems such as molecules as well as extended systems such as crystals and surfaces, in both static and dynamic settings. It is straightforward to install/use and highly competitive with state-of-the-art planewave codes, demonstrating comparable performance on a small number of processors and increasing advantages as the number of processors grows. Notably, SPARC brings solution times down to a few seconds for systems with O(100-500) atoms on large-scale parallel computers, outperforming planewave counterparts by an order of magnitude and more.

MnGa alloys possess a set of properties that give them a wide range of highly interesting applications: from magnetic tunnel junctions [1] to semiconductor-magnetic hybrid devices [2]. In the present study, we report a first-principles investigation on a surface reconstruction based on the L1₀-ordered MnGa (001) with Cu atoms incorporated at the surface. We present a comprehensive analysis of the viability of all the possible reconstructions, as well as the structural, electronic, and magnetic properties of the most stable one: a 1 × 2 reconstruction with a Cu-by-Ga substitution at the surface, similar to the one reported in a previous study [3]. We propose that this structure may pose as a magnetic catalyst, adding yet another interesting application to MnGa alloys.
Acknowledgements
The authors would like to thank the support provided by the following DGAPA-UNAM projects: IN101019 and IA100920, as well as the CONACYT concession A1-S-9070. The calculations were performed in the supercomputing center DGCTIC-UNAM, through the projects LANCAD-UNAM-DGTIC-051 and LANCAD-UNAM-DGTIC-368. We thank Laboratorio Nacional de Supercómputo del Sureste de México, CONACYT, member of the national web of laboratories for providing computational resources and technical support. We would also like to thank the computing center THUBAT KAAL IPICYT for providing its computational resources. Finally, we thank Eduardo Murillo and Aldo Rodríguez Guerrero for their technical support.
References
[1] Q. L. Ma, T. Kubota, S. Mizukami, X. M. Zhang, H. Na-ganuma, M. Oogane, Y. Ando, and T. Miyazaki, Phys. Rev. B 87, 184426 (2013).
[2] M. Tanaka, J. P. Harbison, T. Sands, B. Philips, T. L.Cheeks, J. De Boeck, L. T. Florez, and V. G. Keramidas, Applied Physics Letters 63, 696 (1993).
[3] J. Corbett, J. Guerrero-Sanchez, A. Richard, D. Ingram,N. Takeuchi, and A. Smith, Applied Surface Science 422, 985 (2017).

Defects in crystalline solids play a crucial role in determining properties of materials at the nano, meso- and macroscale, such as the coalescence of vacancies at the nanoscale to form voids and prismatic dislocation loops or diffusion and segregation of solutes to nucleate precipitates. First principles Density Functional Theory (DFT) simulations can provide a detailed understanding of these phenomena, however the number of atoms needed to correctly simulate these systems is often beyond the reach of current DFT codes. We present an accurate and efficient finite-difference formulation and parallel implementation of Linear Scaling Kohn-Sham Density (Operator) Functional Theory (DFT) for non-periodic systems embedded in a bulk environment. We discuss the parallel scalability of the framework, acceleration on the state-of-the-art Graphics Processing Units, and will present the results of applying the framework to study isolated defects in lightweight Magnesium alloys. Finally, we will talk about integrating Machine Learning methods into the framework to further extend the scope of simulations.

Acknowledgements:
Part of this research was carried out when S.G. held a position at the California Institute of Technology. Some of the computations were carried out on the Caltech High Performance Cluster partially supported by a grant from the Gordon and Betty Moore Foundation. This research was also funded in part by the U. S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division and it used resources of the Oak Ridge Leadership Computing Facility, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.

Accurate description of the temporal evolution of correlated electron materials, when perturbed by an external field, is a challenging and sought-after problem in condensed matter physics and material science. This is particularly the case for the ultrafast (femtosecond) response of the system since it may be dominated by electron correlations, unhindered by interactions, for example with the lattice. Suitable ab initio computational methods would thus need to extend beyond the standard density functional theory (DFT). Time-dependent density-functional theory (TDDFT) when equipped with appropriate exchange-correlation (XC) functional offers a way forward. This necessarily implies development of XC functionals beyond the standard local density or the generalized gradient approximations because of the need to treat on-site orbital-resolved correlations and/or time dependent interactions. In this presentation, we first provide some details of our new methodology in which we combine dynamical mean-field theory (DMFT) and TDDFT to examine the ultrafast response of materials in which electron-electron correlations play an important role. We obtain a non-adiabatic XC kernel (potential for nonlinear response) from the electron charge susceptibility (self-energy) by solving an effective Hubbard model using DMFT. This XC kernel or potential when implemented in TDDFT provides a computationally feasible approach for examination of nonequilibrium properties of systems in question. We next demonstrate the viability of the approach through computation of the ultrafast demagnetization of Ni [1], for which we calculate the pre-thermalization of the electron and hole subsystems and related spin-resolved properties. We show that the system magnetic moment reaches its lower, transient-state value (demagnetized) in tens of femtoseconds, and that the charge (electron and hole) pre-thermalization occurs after ~10 femtoseconds, both results in agreement with experimental data. Since at such short time scales electron-electron interaction is the only source of scattering in the system, these results capture the power of the TDDFT+DMFT methodology for isolating the femtosecond physics of strongly correlated materials.
This work was supported in part by DOE grant DE-FG02-07ER46354.
[1] S. R. Acharya, V. Turkowski, G-p Zhang, and T. S. Rahman, Phys. Rev. Lett. 125, 017202 (2020).
https://doi.org/10.1103/PhysRevLett.125.017202

The small-polaron hopping model has been used for several decades for modeling electronic charge transport in oxides. Despite its significance, the model was developed for binary oxides, and its accuracy has not been rigorously tested for higher-order oxides. To investigate this issue, we chose the Mn_xFe_3-xO₄ spinel system, which has exciting electrochemical and catalytic properties, and mixed cation oxidation states that enable us to examine the mechanisms of small-polaron transport. Using a combination of experimental results and DFT+U calculations, we find that the charge transport occurs only between like-cations (Fe/Fe or Mn/Mn). And due to asymmetric hopping barriers and formation energies, we find that the polaron is energetically preferred to the polaron, resulting in an asymmetric contribution of the Mn/Mn pathways. The effect of vacancies near and far from the polaron will be emphasized in this talk.

Reference:
A. Bhargava, R. Eppstein, J. Sun, M. A. Smeaton, H. Paik, L. F. Kourkoutis, D. G. Scholm, M. Caspary Toroker*, R. D. Robinson*, “Breakdown of the small-polaron hopping model in higher-order spinels”, Adv. Mat., 2004490 (2020).

Electronic structure calculations, especially those using density functional theory (DFT), have been very useful in understanding and predicting a wide range of materials properties. Despite the success of DFT, and the tremendous progress in theory and numerical methods over the decades, the following challenges remain. Firstly, the state-of-the-art implementations of DFT suffer from cell-size and geometry limitations, with the widely used codes in solid state physics being limited to periodic geometries and typical simulation domains containing a few thousand electrons. This limits the complexity of materials systems that are accessible to DFT calculations. Further, there are many materials systems (such as strongly-correlated systems) where the widely used model exchange-correlation functionals are inaccurate. Addressing these challenges will enable large-scale quantum-accuracy DFT calculations, and can significantly advance our predictive modeling capabilities of complex materials systems.

This talk will discuss our recent advances towards addressing the aforementioned challenges. In particular, the development of computational methods and numerical algorithms for large-scale real-space DFT calculations using adaptive finite-element discretization will be presented, which form the basis for the recently released DFT-FE open-source code. The computational efficiency, scalability and capability of DFT-FE will be presented, and contrasted with other widely used codes. Further, ongoing efforts, and related thoughts, on developing a framework for a data-driven approach to improve the exchange-correlation description in DFT will be discussed.

Hubbard U corrections are a widely used approach for overcoming the characteristic band gap underestimation in approximate Density Functional Theory (DFT) calculations, without incurring significant computational costs. However, Hubbard U corrections for transition metal oxides (TMOs) are often chosen on a somewhat arbitrary basis and tend to over-correct. This has led to a search for improvements to the U correction such as DFT+U+J, which includes a correction based on the Hund's coupling parameter J that accounts for the dependence of interactions on the spin orientation.

In our research, a standardised first-principles algorithm is implemented for calculating U and J for both metal and oxygen subspaces, based on the linear response to an applied potential. This procedure is benchmarked with several TMO materials to assess the effect of DFT+U+J on the band gap, effective mass, bond lengths and unit cell volume.

We find that in certain materials, first-principles DFT+U+J can give band gap values that are close to experimental values, out-performing PBE and standard DFT+U and, in some cases. even outperforming hybrid functional results. It is also shown that DFT+U+J does not appear to introduce significant distortions in material ionic geometry or band-structures, suggesting that this methodology could be useful for the computationally efficient modelling of TMOs.

One-dimensional nanostructures such as nanotubes, nanowires, nanohelices, nanospirals, and nanocoils have received increased attention over the past three decades due to their fascinating elastic, electronic, vibrational, optical, and thermal properties. In this work, we present a cyclic and helical symmetry-adapted formulation and large-scale parallel implementation of real-space Kohn-Sham DFT for 1D nanostructures, with application to the mechanical and electronic response of carbon nanotubes subject to torsional deformations. Specifically, we derive symmetry-adapted variants for the energy functional, electronic ground state’s variational problem, Kohn-Sham equations, atomic forces, and axial stress, all posed on the fundamental domain, while employing a semilocal exchange-correlation functional and a local electrostatic formulation. Within this framework, we develop a representation for nanotubes of arbitrary chirality subject to external twists. We develop a high-order finite-difference parallel implementation capable of performing cyclic and helical symmetry-adapted Kohn-Sham calculations in both the static and dynamic settings. Using this implementation, we study the mechanical and electronic response of carbon nanotubes to twist-controlled deformations, at both small and large deformations.

Thermodynamics dictates that a crystal under pressure will prefer the structure of lowest enthalpy, provided the temperature is sufficiently low. Structure prediction techniques like random structure searching rely on this insight, optimizing candidate structures to balance the applied pressure, then ranking them by enthalpy. This procedure becomes more complicated when the applied stress is nonhydrostatic. Numerous formulations of enthalpy for nonhydrostatic environments appear in the literature, tracing back to early work by Gibbs, but many are unsuitable for large deformations or have other undesirable features. We illustrate the various difficulties and present our solution in a form suitable for random structure searching, illustrating it with a few simple examples.

Rational design of novel materials for organic light emitting diodes (OLED) requires deep understanding of the processes accompanying device operation at the molecular scale, such as charge and energy transfer, generation and relaxation of excited states, photochemical reactions leading to the degradation of the material, etc. Modern computational chemistry provides versatile tools to study morphology and properties of the OLED materials that often cannot be accurately addressed with experimental techniques.
We present an overview of theoretical methods applied for computational OLED materials development and intend to cover the following subjects:
● prediction of photoluminescence efficiency and CIE coordinates for cyclometalated Platinum and Iridium complexes, which requires evaluation of radiative and non-radiative decay constants, Frank-Condon progressions [1, 2]
● simulations of functional materials using quantum mechanically derived force fields, including modeling of sublimation process of large organometallic compounds, which is an essential step in OLED device fabrication [3-4]
● computational prediction of photochemical stability of OLED materials and more general problem of OLED devices operational stability, including the isotope effects and the problem of charge recombination rate prediction in OLED emitting layers [5-7]
● computational design of OLED materials using generative artificial neural networks: application to triplet-triplet fusion materials

References
[1] Kim. I. et al. “Predicting phosphorescence quantum yield for Pt(II)-based OLED emitters from correlation function approach” J. Phys. Chem. C 123(17), 11140-11150 (2019).

[2] Osipov A. et al. “Tetradentate Pt(II) Phosphors – a Computational Analysis of Nonradiative Decay Rates and Luminescent Efficiency” J. Phys. Chem. C, 124(22), 12039–12048 (2020).

[3] Yakubovich A. et al. “Accurate vapor pressure prediction for large organic molecules utilized in organic light-emitting diodes (OLED)”, J. Chem. Theory Comput. 2020, 16(9), 5845–5851

[4] Odinokov. et al. “Exploiting the Quantum Mechanically Derived Force Field for Functional Materials Simulations”, submitted (2021).

[5] Odinokov A. et al. "Charge recombination in polaron pairs: a key factor for operational stability of blue-phosphorescent light-emitting devices”, Adv. Theory and Sim. 3(8), 2000028 (2020).

[6] Freidzon A. Ya. et al. “Predicting the operational stability of phosphorescent OLED host molecules from first principles: a case study” J. Phys. Chem. C 121, 22422-22433 (2017).

[7] Hyejin Bae et. al. “Protecting benzylic C–H bonds by deuteration doubles the operational lifetime of deep blue Ir-phenylimidazole dopants in phosphorescent OLEDs“, Adv. Opt. Mat. Doi: 10.1002/adom.202100630 (2021).

Accurate molecular simulation requires computationally expensive quantum chemistry models that makes simulating complex material phenomena or large molecules intractable. The relatively recent trend of "data-driven" models intended to overcome this barrier. I will focus on a specific class, based on the Atomic Cluster Expansion (Drautz, 2019) and explain how such models are constructed from symmetry-adapted "features" which can be used in regression schemes including classical linear regression, kernel ridge regression or ANNs, and show some concrete examples highlighting successes and potentials.
I will then discuss some of the more theoretical aspects and try to explain why the ACE approach (and ML in general) has been so succesful: a necessary condition is that there must be a relatively low-dimensional representation of the potential energy landscape. I will focus on three aspects, locality, low body-order and self-consistency, and try to explain rigorously how these occur in the potential energy landscape of ab initio models, and connect this back to the ACE model.

For expansive structures, the individual energy calculations per optimization step becomes non-trivial and therefore rapidly expands the computational resources required for full optimization. In comparison, ML predictions can be agile even on small machines as many computationally expensive matrix diagonalization algorithms are substituted with neural network layers. The model was trained on a wide range of materials, which were systematically stretched/compressed to simulate unoptimized structures. In order to speed up structural optimization, a gradient descent is applied and used to continually adjust atomic positions in order minimize the total energy of the material. This presentation will discuss the universality of the optimizer generated and the overall accuracy achieved.
There are several potential applications of this optimizer, such as speeding up costly energy relaxation calculations. Furthermore, it can be used to as a pre-cleaning tool for novel materials before being further analyzed by other ML models. In addition, a similar strategy of employing stretched/compressed materials was used to improve the performance of bandgap prediction. This presentation will demonstrate the capabilities of the optimizer in improving bandgap predictions for inorganic materials, and how this can be applied in the search for novel blue-UV emitting materials.

One of the main open challenges in computational materials science is the construction of interatomic potentials which achieve quantitative agreement with fully—but infeasible—ab initio models for large-scale problems, possibly involving tens of thousands of atoms. Empirical potentials, arguably the most popular type of interatomic potentials, generally fail in making quantitative predictions. Therefore, they largely remain inadequate for a predictive modeling of multicomponent systems and, thus, for a computational discovery of new materials.

The advent of machine learning interatomic potentials (MLIPs) for small systems of order 10–100 yet holds promise that overcoming these limitations appear within sight (see [1] for a review). On the other hand, a well-known drawback of state-of-the-art MLIPs is their poor ability to extrapolate beyond the training domain. This demands for carefully chosen, problem-dependent training data which can hardly be defined by the user prior to a simulation due to the vast amount of possible atomic neighborhoods to be considered. Therefore, solution algorithms for atomistic problems using MLIPs have to incorporate error estimators which measure the degree of extrapolation in order to “find” the right training data during a simulation.

Here, we present such an algorithm based on active learning and apply it to screw dislocation motion in bcc tungsten [2,3]. Using structural identification methods, we locate the dislocation line, check the extrapolation grade of the MLIP around its core, and construct new training configurations, sufficiently small to be computable with plane-wave density functional theory, if the extrapolation grade exceeds some threshold. We show that the algorithm reproduces classical results from the literature (i.e., core structure, Peierls barrier/stress), in addition to its anticipated application to more complex problems involving interactions of dislocations with other types of defects, e.g., interstitials, which are out of scope with ab initio methods solely.

[1] Y. Zuo et al., 2019. Performance and Cost Assessment of Machine Learning Interatomic Potentials. J. Phys. Chem. A 124, 731–745
[2] E. Podryabinkin, A. Shapeev, 2017. Active learning of linearly parametrized interatomic potentials. Comput. Mater. Sci. 140, 171–180
[3] M. Hodapp, A. Shapeev, 2020. In operando active learning of interatomic interaction during large-scale simulations. Mach. Learn.: Sci. Technol. 1 045005

Although ab initio methods are indispensable tools for predicting properties of materials and simulating chemical processes, the tradeoff between computational efficiency and predictive accuracy limits their application to large systems and long simulation times. The Ultra-fast Force Fields (UF3) framework is an open-source software package for fitting machine-learning models of the energy landscape that are fast to train and evaluate, interpretable in form, and accurate even for relatively small training sets. Specifically, these interatomic potentials are based on a general many-body expansion combined with a cubic B-spline basis and regularized linear regression. We benchmark speed and accuracy using a bcc tungsten dataset, including melting point calculations with thousands of atoms. We find that our potential is two to four orders of magnitude faster than state-of-the-art machine-learning potentials while achieving comparable accuracy. Finally, we discuss the application of the introduced potentials in accelerating structure prediction.

Using machine learning model to yield a machine learning force field (ML-FF) from ab initio calculation data has becoming a powerful approach to simulate material properties with both size and temporal scales much larger than what can be handled by direct ab initio methods. The common ML-FF approach describes the energy of one atom as a function of its neighboring atom positions. This however cannot describe the long range electrostatic interactions. Furthermore, the point charge based electrostatic model in classical force field might not be accurate enough due to the lack of dipole moment, as well as the polarization. Besides, the nonbonding interaction can be rather complicated in terms of its atomic configuration, which makes direct machine learning training difficult. Here, we present an electron charge density based machine learning model to deal with the nonbond interactions. In this model, the electron charge density is fittied, and used for nonbond interactions directly, including electrostatic and exchange-correlation interactions. We also developed an accurate polarization model, describing the polarization energy of a molecule under external potential. We have used battery organic electrolyte as our testing example for this development.

The Local Yield Stress (LYS) method involves directly probing local regions of atoms in atomistic simulations to quantify the incremental stress necessary to induce plastic rearrangements [1,2]. This quantity is key to informing meso-scale and continuum models of plastic deformation and flow in glasses and other amorphous solids. In doing so the LYS seeks to quantify scalar functions of a tensor quantity, the normalized applied shear strain, for a given assemblage of a few tens to hundreds of atoms. While the LYS algorithm is scalable and highly generalizable, it is also computationally intensive. This is particularly true in 3D, which is most relevant for practical application. For this reason, we seek to develop a machine learning methodology to utilize the direct simulation data obtained on a given material to enable efficient and highly generalizable quantification of the LYS. We begin by developing a novel distance metric to quantify the relatedness of atomistic assemblages. We then apply a dimension reduction scheme to the atomistic data in order to determine the most relevant degrees of freedom for our learning algorithm. We make a direct comparison between the atomistic quantification and the learned results in order to validate the method.

[1] Patinet, S., Vandembroucq, D. and Falk, M.L., 2016. Connecting local yield stresses with plastic activity in amorphous solids. Physical review letters, 117(4), p.045501.
[2] Barbot, A., Lerbinger, M., Hernandez-Garcia, A., García-García, R., Falk, M.L., Vandembroucq, D. and Patinet, S., 2018. Local yield stress statistics in model amorphous solids. Physical Review E, 97(3), p.033001.

Recently, quantum dots (QDs) have acquired much attention as promising materials for optical applications. To fully realize the potential of QD, it requires to exhibit a sharp luminescence spectrum. Therefore, many efforts have been taken into minimizing the luminescence line widths of QDs. While it is well known that the luminescence linewidth of a single QD arises from exciton-phonon coupling (EPC), the mechanism and origin of EPC are not fully resolved yet. For example, it is still controversial which mode between acoustic and optical phonon dominates EPC in QDs. A theoretical investigation by density functional theory (DFT) calculations can provide fundamental insights into EPC as well as its impact on the line broadening at an atomic scale. In our previous study, we reported an ab initio method to evaluate luminescence line shapes of QDs due to EPC, combining the DFT and Frank-Condon approximation [1]. However, DFT calculations require large computational time so that the size of QD is restricted to ~500 atoms while real-size QDs are much larger (usually consist of 1,000-100,000 atoms).
In this study, we combine DFT and machine-learning potential to simulate real-size InP/ZnSe QD which consists of ~4000 atoms. The calculated lineshape and its temperature dependence are in excellent agreement with the experiments. In addition, we also find that acoustic phonon modes significantly couple with an exciton in core/shell QDs, which is in stark contrast to the widely accepted perception that exciton-optical phonon coupling dominates the luminescence line broadening.
[1] S. Kang, et al. ACS Appl. Mater. Interfaces, 12, 22012 (2020).

Ab initio modelling of dislocations is essential to gain quantitative insight into a vast range of important material phenomena, from elementary glide or climb mechanisms to predictive models of solute solution strengthening or post irradiation annealing. The simulation data can be used directly in phenomenological micro-structural models or in the construction of empirical interatomic potentials in a multi-scale framework. With some notable exceptions, ab initio size limitations mean that most dislocations cannot be contained in periodic supercells, requiring the use of flexible boundary techniques which either couple to an elastic Green's function or an empirical interatomic potential in QM/MM simulations. In the latter case, the potential should ideally have identical elastic properties as the ab initio medium. A key limitation of these flexible boundary methods is that whilst the ionic positions and forces can be extracted, the total energy becomes ill-defined.
In recent work [1,2], we showed how this information could be used to extract energetic barriers for glide or segregation through the principle of virtual work. A very good agreement was found with special cases treatable with periodic supercells, whilst in general our ab initio results deviate significantly from predictions of empirical potentials, particularly for prismatic (edge) dislocations. In this work, we show how linear-in-descriptor machine learning (ML) potentials [3] can be used to exactly match ab initio elastic properties, typically leading to efficiencies over "traditional" empirical potentials during QM/MM relaxations. We then go further, using the ionic positions and forces to train the ML potential, showing that the resultant "classical" migration barriers are in excellent agreement with that obtained from QM/MM. Our work shows that QM/MM methods can qualitatively expand the range of atomic structures in the training database of machine learning interatomic potentials, with significant impact on their transferability, or predictive power, for materials modelling.

References:
[1] Swinburne, T. D., & Kermode, J. R. (2017). Computing energy barriers for rare events from hybrid quantum/classical simulations through the virtual work principle. Phys. Rev. B, 96(14), 144102.
[2] Grigorev, P., Swinburne, T. D., & Kermode, J. R. (2020). Hybrid quantum/classical study of hydrogen-decorated screw dislocations in tungsten: Ultrafast pipe diffusion, core reconstruction, and effects on glide mechanism. Phys. Rev. Materials, 4(2), 23601.
[3] Dérès J, Goryaeva A. M., Lapointe C., Grigorev P., Swinburne, T. D., Kermode, J. R., Ventelon L., Baima J., and Marinica M.-C. Efficient and transferable machine learning potentials for the simulation of crystaldefects in bcc Fe and W. Submitted to Phys. Rev. Materials

Many present day technological challenges related to battery materials, electro-catalysis, fuel cells, corrosion and others, centre around understanding and controlling processes at solid/liquid interfaces. A targeted design of desired functionalities critically relies on understanding and quantifying the underlying fundamental mechanisms. However, understanding them is equally challenging to both theoretical modelling and experimental characterisation.

The presentation will discuss our recent advances achieved by developing a novel potentiostat design [Surendralal et al., Phys. Rev. Lett. 120, 246801 (2018)] and a canonical thermopotentiostat [Deißenbeck et al., Phys. Rev. Lett. 126, 136803 (2021)] approach. These enable us to control the electrode potential of the system by tuning the excess charge of the working electrode and use ab-initio molecular dynamics simulations to study solid/liquid interfaces under realistic conditions of applied bias. The opportunities opened by the availability of these new approaches will be discussed for a few prototype examples related to electro-catalysis and corrosion: The mechanisms leading to the onset of the hydrogen evolution reaction (HER) in the H/Pt/H₂O system [Surendralal et al., Phys. Rev. Lett. 126, 166802 (2021)] and corrosion at the Mg/water interfaces.

We illustrate how the Hill relation and the notion of quasi-stationary distribution can be used to analyse the error introduced by many algorithms that have been proposed in the literature, in particular in molecular dynamics, to compute mean reaction times between metastable states for Markov processes. The theoretical findings are illustrated on various examples demonstrating the sharpness of the error analysis as well as the applicability of our study to elliptic diffusions. This is a joint work with Manon Baudel and Arnaud Guyader. Reference: M. Baudel, A. Guyader, and T. Lelièvre, On the Hill relation and the mean reaction time for metastable processes, https://arxiv.org/abs/2008.09790 .

Understanding hydrogen transport is vital to industries focused on discovering new materials for energy storage and corrosion mitigation. However, knowledge of the physical nature of hydrogen’s diffusion pathways is often limited, especially for materials that exhibit a multitude of phases/defect classes. These materials typically have three rate-limiting structural domains: (1) bulk (2) grain boundaries, and (3) surfaces. In this work we have chosen to study hydrogen diffusion through titania due to its importance in the aforementioned application spaces. Here, we aim to understand diffusion through titania grain boundaries, which are approximated via the amorphous phase. Density functional theory (DFT) was used to calculate thousands of activation energies of hydrogen diffusion in the amorphous phase via nudged elastic band (NEB) calculations using an automated hydrogen pathway generation scheme. Amorphous structures, on the order of tens of nanometers, were generated using a machine learning force field via classical molecular dynamics (MD). Markov chains (MC) were generated using the MD-derived atomic structures as their reference. The MC edge weights, which employ the activation energy between hydrogen hopping sites, were calculated via a Gaussian process regression model, which was trained on the DFT-NEB data. Kinetic Monte Carlo simulations were then performed over a variety of temperatures and amorphous structures. Both stoichiometric and non-stoichiometric amorphous structures were chosen to understand the relationship between hydrogen diffusion and oxygen concentration. The machine learning-assisted adaptive Markov Chain Kinetic Monte Carlo scheme employed in this work allows us to quantitatively describe hydrogen diffusion rates for systems without pre-defined sites and/or activation energies. This work sets the stage for one to perform long time-scale and/or length-scale simulations of transport phenomena for practically any material system, removing the need to simplify the system either structurally and/or kinetically.

Global warming caused by the excess use of fossil fuels has become a critical issue in modern society. The proton-exchange membrane fuel cell (PEMFC) has attracted wide attention as the environmentally friendly renewable energy device that is used in hydrogen-fueled vehicles. The performance of PEMFC relies on the nanoscale electrocatalyst on the cathode electrode. Over the last several decades, experimental and computational researches have focused on the sluggish oxygen reduction reaction (ORR) and the cost of a precious metal catalyst. However, as the activity of the catalyst increases, their low durability problem also rises. Therefore, an atomistic picture of the degradation of nanoscale catalysts is needed to improve the stability of nanoscale catalysts.
While several works describe the long-term degradation evolution with kinetic Monte Carlo (kMC) simulation, they are limited due to the inaccurate classical potential or computationally expensive quantum calculation. Recently, neural network potential (NNP) has been applied to various systems, delivering the accuracy of quantum calculation at a much lower cost.
In this talk, we combine NNP and kMC simulation to explore the degradation of various shapes of Pt alloy nanoparticles with 3–5 nm under electric bias. We found that our simulation results agree with the experimental observation. This work would pave the way to design more durable nanoparticle electrocatalysts.

Titanium nitride (TiN) belongs to the family of refractory transition metal nitrides and represents useful properties such as chemical inertness, thermodynamic stability, extreme hardness, good thermal and electrical conductivity. Atomic layer deposition (ALD) of TiN has been used for protecting layers in tungsten as the minimum feature size of the devices decreased and a high aspect ratio was required. Accordingly, experimental and theoretical research on the TiN-ALD thin film process using various reducing agents (NH₃, H₂) based on the TiCl₄ precursor was published. However, observations on the sub-nanometer scale surface reaction are limited due to the characteristics of the ALD process that grows at the atomic scale. In the case of previous theoretical studies, there is a gap in understanding the actual process because it is limited to indirect analysis through thermodynamic energy based on density functional theory (DFT) calculation. So in this study, we introduce the multi-scale modeling for TiN-ALD. We investigate the overall growth of the TiN thin films grown via TiCl₄/NH₃. Multi-scale modeling with DFT and on-lattice kinetic Monte Carlo (kMC) is adopted to reproduce the whole ALD process. To determine the most probable pathway, all individual reactions are considered and calculated with DFT, and a set of corresponding activation energies are implemented in the kMC code as a transition probability for the surface reactions on the TiN (111) surface. To verify the compatibility with the experiment ALD process, the growth rate and mass gain according to the process conditions (temperature, pulse time, pressure) were compared with the simulation results. The temperature-dependent growth profile obtained from our model is in good qualitative agreement with the experimental data. These kinetic results should help to explain TiN ALD growth rates observed at different reactant exposures and substrate temperatures.

Atomistic structure search for organic/inorganic heterostructures is made complex by the many degrees of freedom and the need for accurate but costly density-functional theory (DFT) simulations. To accelerate structure determination in heterogenous functional materials across lengthscales, we employ active learning with the Bayesian Optimization Structure Search (BOSS) approach [1].

Here, Bayesian optimization is employed to build N-dimensional surrogate models for the energy or property landscapes and infer global minima. The models are iteratively refined by sequentially sampling DFT data points that are promising and/or have high information content. Representing heterogenous materials with compact chemical ‘building blocks’ allowed us to build in prior knowledge and reduce search dimensionality. The uncertainty-led exploration/exploitation sampling strategy delivers global minima with modest sampling, but also ensures visits to less favorable regions of phase space to gather information on rare events and energy barriers.

We applied this active learning scheme to study adsorption at the organic/inorganic interfaces: C₆₀ on TiO₂(101) anatase [1] and camphor on Cu(111) [2,3]. At larger lengthscales, we investigated the substrate-decoupled F4TCNQ-TTF charge transfer complexes. Global minima configurations in 6 dimensions are identified with reasonable computational efficiency. BOSS produces chemically-intuitive adsorption energy landscapes that are parsed for local minima and the minimum energy paths between them, allowing us to extract complex barrier-related atomistic pathways. With a recent batch implementation for active learning, BOSS can make use of exascale computing resources to solve largescale structural problems without sacrificing quantum-mechanical accuracy.

[1] M. Todorović, M.U. Gutmann, J. Corander and P. Rinke, ‘Efficient Bayesian Inference of Atomistic Structure in Complex Functional Materials’, npj Comput. Mater. 5, 35 (2019).
[2] J. Järvi, P. Rinke and M. Todorović, ‘Detecting stable adsorbates of (1S)-camphor on Cu(111) with Bayesian optimization’ Beilstein J. Nanotechnol. 11, 1577 (2020).
[3] J. Järvi, B. Aldritt, O. Krejčí, M. Todorović, P. Liljeroth and P. Rinke, ‘Integrating Bayesian Inference with Scanning Probe Experiments for Robust Identification of Surface Adsorbate Configurations’, Adv. Func. Mater. 2010853 (2021).

Uncertainty quantification (UQ) methods such as Bayesian methods have been recently used to quantify the uncertainty associated with the force-field parameters and their predictions. However, the reliability and accuracy of these UQ algorithms depend on the sampling methods used to perform the analysis. Here, as a test case, we have performed UQ on a new all-atom (AA) embedded atom method (EAM) potential of gold developed in our group. The AA model itself was aimed at reproducing physical, thermodynamic, and mechanical properties of interest obtained through experimental methods and/or density functional theory. The robustness and parametric sensitivity of the gold model was analyzed using Sobol sequencing and calculating sensitivity indices. Bayesian Parameter Estimation with UQ was performed using different sampling methods including Metropolis-Hastings (MH), Ensemble-slice (ESS), Uniform, and Grid as implemented in CheKiPEUQ, an open-source Bayesian parameter estimation package. This allowed us to investigate the performance and dependence of Bayesian predictions on these sampling methods. The sampling process was accelerated through surrogate modeling by developing Gaussian Process Regression (GPR) models trained with data obtained from molecular dynamics (MD) simulations. The Bayesian framework was used to obtain the highest posterior density interval and the highest probability parameter set, which better reproduces the properties with respect to the confidence levels of experimental and/or DFT observations.

Modeling diffusion from the atomic scale is necessary because experimental studies are only feasible at high temperatures and in most systems, it is not possible to measure the full Onsager matrix. The problem is mostly solved in pure and dilute alloys, for instance using the KineCluE code. Computing the full Onsager matrix in concentrated alloys remains a challenge due to the vastness and complexity of the configuration space and the connected ones associated with jump sequences. Current approaches rely either on mean-field approximations or atomic kinetic Monte Carlo simulations, however each of them has its own limitations; with the former, kinetic correlations are not fully taken into account while with the latter statistical errors arise, especially when calculating off-diagonal transport coefficients. Therefore, a comprehensive deterministic diffusion model is needed. In this work, we aim to extend the kinetic cluster expansion formalism to concentrated alloys in order to get an efficient evaluation of transport coefficients, time life and mobility of clusters in concentrated alloys based on atomic-scale data. Our formalism relies on the Self-Consistent Mean Field Theory and it is implemented in the KineCluE code. We treat the local configuration of the diffusing species exactly, while further away atoms are replaced by a homogeneous mean-field. This allows us to treat short-range kinetic correlations exactly and to construct the long-range kinetic paths under controlled approximations while reducing the computational load and keeping the calculations feasible. In the case of a random alloy, we compare our model’s performance to previous analytical theories and Monte Carlo simulations. The model will also be tested in the case of Fe-Cr concentrated alloys where thermodynamic short-range order cannot be neglected, and the results will be compared to Monte Carlo simulations as well as experimental studies.

In the context of the high-dimensionality and multimodality of the potentials encountered in chemical physics, massive efforts have been devoted into the development of efficient simulation methodologies. Since the first works of Molecular dynamics (MD) and Markov-chain Monte Carlo (MCMC) methods, bidimensional particle systems, in particular with hard-core interactions, have played the role of a challenging test bed. In spite of their apparent simplicity, they exhibit a rich behavior where topological defects bind or unbind themselves in pair, ruling by doing so the transition phases. I will first review successful MCMC simulation strategies in sphere systems. These methods are based on the introduction of non-reversibility in the underlying stochastic sampling processes, upgrading them into piecewise deterministic Markov processes (PDMP) and competitors to their MD counterparts. Then I will present how to adapt these strategies in the context of anisotropic particles by exploiting recent progress in PDMP-based sampling.

Monte Carlo modeling of secondary electron emission from metal surfaces usually requires information on material properties obtained from experimental measurements, in particular the energy loss function, a quantity derived from the dielectric function is needed, as well as the Fermi energy and the work function. In addition, the momentum space representation of the energy loss function is often approximated with the use of extrapolation schemes. Here we present a new approach and a first open-source Monte Carlo code for secondary electron yield modeling [1], where all the necessary quantities are calculated from first principles, density functional theory (DFT) methods. In particular, the energy and momentum dependent dielectric function is obtained from time-dependent DFT calculations explicitly, removing the necessity for experimentally measured quantities as well as momentum space extrapolation. In addition DFT calculated density of states is used in secondary electron generation, and DFT obtained work functions are used. This approach allows for predictive modeling of systems where experimental measurements are not available or are very challenging to obtain with high accuracy, such as metal alloys [2]. Successful verification of the approach on well known systems such as Cu and Al, encouraged calculations of secondary electron emission from elemental metals across the whole periodic table. The calculations revealed clear periodic trends in the maximum secondary electron yield, which has not been shown previously, and provided the largest computational database of secondary electron yield curves for elemental metals up to date.

[1] M. P. Polak and D. Morgan, “MAST-SEY: MAterial Simulation Toolkit for Secondary Electron Yield. A monte carlo approach to secondary electron emission based on complex dielectric functions”, Computational Materials Science 193, 110281 (2021), https://doi.org/10.1016/j.commatsci.2021.110281.
[2] R. E. Gutierrez, I. Matanovic, M. P. Polak, R. S. Johnson, D. Morgan, and E. Schamiloglu, "First principles inelastic mean free paths coupled with Monte Carlo simulation of secondary electron yield of Cu-Ni, Cu-Zn, and Mo-Li", Journal of Applied Physics 129, 175105 (2021), https://doi.org/10.1063/5.0049522

Quasi-stationary distributions can be seen as the first eigenvector associated with the generator of the stochastic differential equation at hand, on a domain with Dirichlet boundary conditions (which corresponds to absorbing boundary conditions at the level of the underlying stochastic processes). Many results on the quasi-stationary distribution hold for non degenerate stochastic dynamics, whose associated generator is elliptic. The case of degenerate dynamics is less clear. In this work, together with T. Lelièvre and J. Reygner (Ecole des Ponts, France) we generalize well-known results on the probabilistic representation of solutions to parabolic equations on bounded domains to the so-called kinetic Fokker Planck equation on bounded domains in positions, with absorbing boundary conditions. Furthermore, a Harnack inequality, as well as a maximum principle, is provided for solutions to this kinetic Fokker-Planck equation, as well as the existence of a smooth transition density for the associ- ated absorbed Langevin dynamics. The continuity of this transition density at the boundary is studied as well as the compactness, in various functional spaces, of the associated semigroup. This work is a cornerstone to prove the consistency of some algorithms used to simulate metastable trajectories of the Langevin dynamics, for example the Parallel Replica algorithm.

In recent years we have been working on adapting a novel Bayesian statistical approach, called nested sampling, for studying atomistic systems. Nested sampling automatically generates all the relevant atomic configurations, unhindered by high barriers, and one of its most appealing advantages is that the global partition function can be calculated very easily, thus thermodynamic properties, such as the heat capacity or compressibility becomes accessible.

Nested sampling is a top-down sampling technique, starting from the high-energy region of the potential energy landscape (gas phase) and progressing towards the ground state structure through a series of nested energy levels, estimating the corresponding phase space volume of each. This way the method samples the different basins proportional to their volume, and instead of providing an exhaustive list of the local minima, it identifies the thermodynamically most relevant states without any prior knowledge of the structures or phase transitions. This means that unlike other methods, nested sampling may be fully automated, allowing high-throughput calculations of novel materials.

The use and advantages of the nested sampling method has been demonstrated in sampling the potential energy landscape of simple model system, calculating the pressure-temperature phase diagram of several metals and alloys[1]. Studying phase transitions and structural properties of small Lennard-Jones[2] and water clusters[3] highlighted the importance of identifying those local minimum basins which contribute the most to the total partition function. Our recent work on the core softened double ramp potential, the Jagla-model, has demonstrated that nested sampling plays an important role in discovering high temperature thermodynamically stable phases, redrawing what we knew about the system’s liquid-liquid transitions[4].

[1] RJN Baldock, LB Partay, AP Bartok, MC Payne and G Csanyi, Phys. Rev. B 93, 174108 (2016); RJN Baldock, N Bernstein, KM Salerno, LB Partay and G Csanyi, Phys. Rev. E 96, 043311 (2017); LB Partay, Comp. Mat. Sci. 149 153 (2018)
[2] LB Partay, AP Bartok and G Csanyi, J. Phys. Chem. B 114 10502 (2010)
[3] J Dorrell and LB Partay, Phys. Chem. Chem. Phys. 21 7305 (2019)
[4] AP Bartok, G Hantal and LB Partay, arXiv:2103.03406

Widely considered a workhorse in materials simulation, MD has struggled to maintain strong-scalability as it only weak-scales with traditional parallelization methods. One method aimed towards improving the strong-scalability of MD is Parallel Trajectory Splicing (ParSplice); a specialized MD method designed to accelerate the dynamics of rare-event systems without compromising the accuracy of MD. It works by generating a large number of independently-generated trajectory "segments'' in such a way that they can later be assembled into a single statistically-correct state-to-state trajectory. Due to the impossibility of deterministically forecasting the future evolution of the trajectory (which would be required to perfectly assign the execution of segments), ParSplice relies on stochastic speculation as it assigns segments via a KMC-like model that is parameterized on-the-fly. In this talk we present a method for extracting an estimate of the likelihood that any given segment will be used as part of the simulation. Then, using these probabilities, we derive an optimization framework for dynamically allocating resources among the set of possible tasks (i.e, potential segments). We detail the improvements in performance resulting from our optimized framework as compared to more naive allocation schemes including the maximum-task-throughput and minimum-time allocations. Performance results are illustrated through the use of several models of varying complexity, lending to fully interpretable results. These models enable an understanding as to when and why other allocation schemes fail, thus further highlighting the versatility/importance of our optimized allocation method. We observe that this new strategy can improve throughput up to 20x, thus greatly improving the scalability of a traditionally scale-limited application.

Irving and Kirkwood (IK) derived the transport equations from the principles of classical statistical mechanics using the Dirac delta to define local densities. Thereby, formulas for fluxes (local stress and heat flux) were obtained as point functions in terms of molecular variables. The IK formalism has inspired numerous formulations. Many of the later developments, however, believed it more rigorous to replace the Dirac delta with a continuous volume-weighted averaging function and subsequently define fluxes as a volume density. Although these volume-averaged (VA) flux formulas have dominated the literature for decades and are widely implemented in popular molecular dynamics (MD) software, they are a fundamental departure from the well-established physical concept of fluxes.
In this talk, we review the historical developments that led to the unified physical concept of transport fluxes. We then use MD simulations to show that replacing the Dirac delta with a volume weighting function changes the fundamental nature of fluxes as a surface density. As a result, the dynamic balance between the rate of change of the total momentum or energy in a volume element and the fluxes across the surface boundary cannot be established. This then leads to the failure of VA flux formulas in satisfying the momentum and energy conservation laws or typical transport boundary conditions. We also use two different approaches to derive fluxes for general many-body potentials. The results show that atomistic formulas for fluxes can be fully consistent with the physical definitions of fluxes and conservation laws.

VLSI (Very-Large Scale Integration) industry is shifting the paradigm from microelectronics to nanoelectronics so that Moore’s law can be extended. Molecular electronics directly utilizes individual molecules to construct circuits and systems, leading to ultra-high efficiency and device density. For digital molecular electronics, a suitable molecule with controllable and stable electrical switching between "ON" and "OFF" states is desired. Among various potential candidates, bistable [2]rotaxane has attracted extensive attention from researchers due to its repeatable switching characteristics and the ability to be driven electrically. Many studies have been performed on characterizing the [2]rotaxane molecule, its I-V characteristics, and molecular dynamics mechanism. In this paper, we reported our findings of molecular simulation of [2]rotaxane utilizing Nanoscale Molecular Dynamics (NAMD) and Visual Molecular Dynamics (VMD) platforms. The simulation illustrates the switching mechanism of [2]rotaxane between its two stable states, and what factors might affect the switching process. The simulation shows the switching process of [2]rotaxane is repeatable and controllable, resulting it to be an excellent candidate for molecular memory application. The signal “crosstalk” – interference between neighboring [2]rotaxane molecules, is also investigated with simulation. A mathematical model is proposed to describe the minimum pitch of the memory array to ensure the reliability of the data stored in the memory. Based on this, a set of design rules, the analogy to DRC (Design Rule Check) in CMOS VLSI, are proposed to tolerate the process variation and signal interference of the [2]rotaxane-based molecular memory. The simulation results verify [2]rotaxane as a promising candidate for molecular memory application.

The need to develop environmentally friendly and energy efficient technology for gas separations can be fulfilled by polymer mixed matrix membranes (MMMs) of metal-organic frameworks (MOFs). In these MMMs, polymers, MOFs and the polymer-MOF interface are known to play an important role in determining the gas separation mechanism. However, fundamental understanding of the parameters determining the compatibility of the polymer-MOF interface is still unknown. In this study, we have performed all-atom (AA) molecular dynamics (MD) simulation of a model polymer-MOF MMMs of hydrophobic/hydrophilic polymers with IRMOF-1 and IRMOF-10 as fillers. These systems include the polymer matrix of hydrophobic polymers viz. polyvinylidene fluoride (PVDF), polyethylene (PE), and hydrophilic polymers viz. polymethylmethacrylate (PMMA), polyethylene glycol (PEG). To understand the effect of MOF functionalization on the polymer-MOF interface, we decorated IRMOF-1 with the four functional groups viz. methyl (-CH₃), amine (-NH₂), hydroxy (-OH) and fluorine (-F). A stronger structural correlation for the polymer-MOF interface was observed in the systems with functionalized MOFs compared to unfunctionalized MOFs embedded in hydrophilic (or hydrophobic) polymer matrix. Thus, different functional groups can be added to control the compatibility of polymer-MOF interfaces.

Molecular dynamics simulations are currently undergoing a revolution thanks to modern machine learning tools that can provide very accurate short-range interatomic potentials. Unfortunately, without a dramatic increase in the computational cost, these new techniques fail to account for the critically important long-range electrostatic interactions and their associated polarization and charge relaxations – rendering them practically useless for many real-world problems. In this talk I will show how the computational overhead associated with the long-range charge equilibration problem in molecular dynamics simulations based on non-linear charge relaxation models can be avoided. The technique relies on the construction of approximate shadow energy functionals whose relaxed ground state charge distributions closely follow the fully equilibrated solutions of the corresponding exact non-linear charge relaxation models. This makes it possible to perform accurate molecular dynamics simulations with non-linear charge relaxation models using only a single Coulomb summation per time step.

In the last couple of decades, high entropy alloys (HEAs) or multi-element alloys are emerging as potential candidates for many diversified structural applications. Molecular dynamics based simulations were performed in conjunction with 2NN MEAM potential to study the fracture behaviour in ternary, quaternary, and quinary alloy configurations. Simulations were performed after aligning the crack plane with three different principle planes of FCC crystal. Deformation in each of the orientations was primarily governed by dislocations, twinning, and stacking faults emanating front the crack tips or from the surface of the crack plane. Medium and high entropy alloys have shown higher resistance to fracture and were capable of retaining strength up to a larger extent as compared to low entropy alloy and pure Ni crystal. The higher dislocation density and uniform distribution of dislocation in medium and high entropy alloys lead to plasticity in the crystal after the onset of yielding. Blunting of the crack tip in medium and high entropy alloys was predicted in orientation 1 and 2, whereas crack propagates in orientation 3. This study will help in elucidating the capabilities of medium and high entropy alloys in resisting the opening mode of crack propagation, which can be further utilized to enhance the diversified structural applications of these emerging multi-elemental alloys.

Beam damage in High-Resolution Transmission Electron Microscopy (HRTEM) is poorly understood theoretically, and not yet well described phenomenologically. For example, it is not yet known whether it is the total dose (electrons per area) or the dose rate (electrons per area per time) that is most critical for inducing beam damage [1,2]. In order to study the dependence of the total dose, a molecular dynamics simulation can be an appropriate tool to understand the atomic level behavior, when coupled to experimental observations.
Using molecular dynamics, we study the beam damage in gold nanoparticles caused by high-energy electron beam irradiation. We assume that the main damage mechanism in metallic nanoparticles is knock-on damage, where an energetic electron transfers momentum and energy to an atom through Rutherford scattering.
In this work, we focus on the simplest mechanism for beam damage: direct transfer of momentum between the energetic electrons and the nuclei through Rutherford scattering. Several molecular dynamics studies have previously focused on electron head-on collision in carbon nanomaterials [3,4]. Here, we extend these methods to include collisions with non-zero impact parameter, and build a model including both complete removal of atoms as well as beam-induced surface diffusion on the nanoparticle, in the context of HRTEM. We use a stochastic approach for the momentum transfer to the atoms in the nanoparticle, focusing on momentum transfer to the easily removed surface atoms. In this way, for each nanoparticle size, we get a rate of atom removal and a rate of induced surface diffusion proportional to the beam intensity. We compare these rates with observation of beam induced surface diffusion and beam induced damage in gold nanoparticles on a cerium dioxide (CeO₂) substrate, as observed in HRTEM. The analysis of the experimental data leverages convolutional neural networks to identify when the structure changes [5].
If the momentum transferred to an atom is insufficient (such as non-zero impact parameters) to remove the atom from the nanoparticle, the induced motion will be primarily in the beam direction. In this case the nanoparticle may heal through thermal diffusion processes. We model this qualitatively through a short molecular dynamics annealing step between each energy transfer event. This may lead to a dependency on not only the total dose, but also on the dose rate.
In addition, effects of the dose rate are expected to be shown predominantly through the heating of the nanoparticle due to the increase of the events per unit of time as we increase the dose rate. In addition to heating from the many Rutherford scattering events with too low momentum transfer to inflict direct damage, there may also be a heating effect due to interactions between the beam electrons and the conduction electrons of the metals. Simple models can be made for these effects, and rough estimates of the heating of the nanoparticle can then be used to extend the current molecular dynamics simulations to higher temperatures, where beam damage will occur more readily.




REFERENCES:

[1] Kisielowski, Christian, et al. Physical Review B 88.2 (2013): 024305.
[2] Kisielowski, C., et al. Advanced structural and chemical imaging 2.1 (2016): 1-14.
[3] Yasuda, Masaaki, et al. MRS Online Proceedings Library 1700.1 (2014): 29-35.
[4] Pregler, Sharon K., and Susan B. Sinnott. Physical Review B 73.22 (2006): 224106.
[5] Madsen, Jacob, et al. Advanced Theory and Simulations 1.8 (2018): 1800037.

We simulate structure in the vicinity of different size nanovoids using a new variant of the Molecular Statics, wherein atomic structure in the vicinity of nanovoids and the parameters that define the displacements of atoms placed in elastic continuum around main computation cell are determined in a self-consistent manner. Atomic structure are calculated in vicinity of nanovoids for bcc and fcc metals. Should pay attention that in bcc structures displacements of atoms located in the direction <100> or close to it are positive, whereas the atom displacements located in other directions is negative [1,2,3]. For fcc structures all displacements are negative but may vary significantly in value depending on the direction. Strain fields calculated based on the displacement fields lead to anisotropy of vacancy fluxes. Then we simulate structure of the surface of different crystallographic planes. Kinetic equations for shifting rate of void surface elements are obtained [2] that take into account the dependence of the vacancy flux on strain fields and the surface energies of the crystallographic planes that located normal to the void surface. The shifting rates of void surface elements are calculated for different crystallographic directions in a wide range of temperatures and different supersaturations of vacancies. Rates of displacements in different crystallographic directions for bcc and fcc metals differ significantly, and for the <100> direction they are greatly lower than for other directions. This effect can lead to a change in the shape of the initially spherical nanovoids and cause their transformation into pores of a cuboid shape. The influence of the calculated values of the surface energies for various crystallographic planes on the shifting rates in the corresponding directions is also evaluated.
[1] A.V. Nazarov, I. V. Ershova, and Y. S. Volodin, KnE Materials Science, p. 451, (2018) DOI 10.18502/50.
[2] Andrei V. Nazarov, Alexander A.Mikheev, Aleksey P. Melnikov, Journal of Nuclear Materials, 532, 152067 (2020).
[3] A V Nazarov et al 2020 IOP Conf. Ser.: Mater. Sci. Eng. 1005 012026-6
DOI 10.1088/1757-899X/1005/1/012026

At present, Molecular Dynamics (MD) is the most suitable approach to gain information at the atomic scale. However, a simulation of diffusion by MD faces some difficulties in implementation. In particular, we noticed that when modeling diffusion by MD and molecular statics (MS) using the same potentials, the found values of the migration energy differ markedly [1]. In our work we discuss features of atomic jump simulation. Then we present a new approach (approach of natural thermostat) that fully takes into account the fluctuating character of atomic jumps at the constant temperature MD simulation. In addition the combination of MD and modified method of MS [2] allowed within this model to take into account thermal expansion of the lattice and long-range elastic displacement of the atoms in the vicinity of defects. With the developed model, we study features of atom diffusion in iron by vacancy mechanism at different temperatures using a many-body potential [3]. Migration energy and pre-exponential factor are obtained and a comparison of the calculated diffusion parameters with the results of the other authors’ work whose used various modeling methods is carried out. We then developed our approach to take into account the effect of pressure on diffusion jumps of atoms. This approach made it possible to calculate the temperature dependence of the migration volume. In addition, the model developed by us makes it possible to select and construct the trajectories of motion of atoms making a diffusion jump, and therefor to analyze the features of various types of diffusion jumps, in particular, to detect such events as double jumps of atoms and jumps leading to the return of an atom to a vacancy. Analysis of a large number of diffusion jump trajectories allows us to estimate the contribution of double jumps in the diffusion mobility at different temperatures.

[1] Mikhin A.G., Osetsky Yu.N., J. Phys.: Condens. Matter, –1993 –Vol.5 – p.9121.
[2] Valikova I.V., Nazarov A.V., Defect Diffus. Forum, – 2008 – Vol. 277 – p.125.
[3] Ackland G.J., Bacon D.J., Calder A.F., Harry T. , Phil. Mag. A. – 1997. – Vol. 75. – N 3. –p. 713.

Interface diffusion and interfacial slip are two key processes controlling the high-temperature creep deformation of metal-matrix composites. Fundamental understanding of these processes remains poor. Atomistic simulations are a promising route toward understanding, but reliable methodology for measuring kinetic coefficients and disentangling competing mechanisms has not been developed in this context. In this work, we compare interface diffusion and sliding coefficients in pure Al and Si grain boundaries to those in AlSi phase boundaries using direct molecular dynamic simulations at high and intermediate homologous temperatures. Diffusion along AlSi phase boundaries is found to be quite fast, depending on interface orientation and structure, with activation energies on par with GB diffusion in pure Al. Clusters of highly mobile Si atoms at these phase boundaries often have activation energies similar to Al. Interface diffusion and interface slip are found to be distinct, but correlated processes with characteristic fractal dimensions, cluster size distributions, and nucleation and growth dynamics. Universal and system dependent features of atomic motion during diffusion and sliding are identified over a large parameter space including multiple interface geometries, temperatures and strain rates. Key results are verified with an accurate PINN potential.

While it is well-established that atoms diffuse along interfaces much faster than in the perfect lattice, atomistic understanding of this “short-circuit” diffusion phenomenon remains incomplete. This talk will discuss several topics related to interface diffusion: atomic mechanisms of grain boundary (GB) diffusion, the effect of GB segregation on the GB diffusion rate, and the role of GB diffusion in the solute drag phenomenon. These effects are studied by atomistic computer simulations combining molecular dynamics, Monte Carlo methods, and other computational approaches. It is shown that solute segregation can cause either acceleration or retardation of GB diffusion, depending on the thermodynamics of the system, GB structure, and the segregation mechanism. In a solute-drag study, the drag force has been investigated systematically as a function of temperature, alloy composition, and the GB velocity. It is shown that solute diffusion along the GB can dramatically impact the GB mobility by redistributing the solute atoms non-uniformly in the form of nano-clusters or islands and causing additional pinning. GB diffusion coefficients and diffusion mechanisms are compared with those for metal/semiconductor phase boundaries with different orientation relationships between the phases.

Dynamic charge distribution and polarization effects are often difficult to describe using traditional force field models, especially for crystalline and amorphous systems. Reactive molecular dynamics has become a panacea for describing the potential chemistry between atoms and molecules, and the environment-dependent charge distribution and polarization effects. Charge equilibration as well reactive force field models such as ReaxFF can include these effects in large systems. However, the numerical implementation of these models include Coulomb interaction up to a “short-distance” cutoff for computational speed, limiting their application to devices with unrealistic small sizes. Although designing an effective atomic charge model to reproduce the electrostatic fields in electrochemical systems has always been a challenge, ignoring long-range Coulomb interaction affects the ability to describe charge distribution and polarization effects in electrochemical systems. We include long-range effects to partial charge calculations (QEq) using shielded long-range Coulomb potential and then evaluate long-range effect between charged particles. By including long-range effect, we anticipate better account of Coulomb interaction for dynamic charge distribution and polarization in atomic and molecular systems. The accuracy of our model to compute charges on a number of metal organic frameworks and simple systems compare to conventional QEq and quantum mechanics (QM) shows that charges computed without long-range are slightly overestimated. A test of our QEq technique against standard QEq on a simple capacitor configuration under external potential shows some significant differences between charges computed with/without long-range Coulomb. The difference results from long-range influence on ions of the capacitor. The difference indicates clearly that the excluded ions via short-distance cut-off treatment of conventional QEq coupled molecular dynamics affect the interaction potential. Further improved investigation of our approach shows that inclusion of orbital overlap function and/or bond-dependent electronegativity as in split charge (SQE) method correctly exhibits dissociation property for largely separated atoms. The later improvement would enable the atomic description of electrochemical systems involving bond breaking while preserving long-range Coulomb interaction effect on charged electrodes throughout the dielectric layer such as in battery and emerging solid-state memory.

The present study aims to investigate the effect of two-dimensional (2D) inorganic nanofiller boron nitride nanosheet (BNNS) on the shock compressive response of polyethylene (PE) based nanocomposites. A non-equilibrium molecular dynamics (NEMD) based simulations were performed in conjunction with reactive force field (Reaxff) parameters to capture the interatomic interactions within 2D nanofiller BNNS and polyethylene chains, whereas interfacial properties were triggered with the help of non-bonded interactions. Nanocomposite reinforced with monolayer BNNS is predicted to be a superior nanofiller for attenuating the shock strength across the interface BNNS/PE nanocomposites as compared to pristine PE. In continuation with the monolayer reinforced nanocomposites, compressive shock hugoniot response of dispersed BNNS reinforced nanocomposites were also predicted and compared with the monolayer BNNS reinforced nanocomposites. It was reported that dispersed BNNS significantly increases the energy dissipation capability via deformation and temperature rise of the dispersed BNNS nanosheet. These present study outcomes will help developing inorganic nanofiller reinforced nanocomposites with superior shock mitigation fin lightweight armor and space structures.

We report a theoretical development of a fully anharmonic method for the calculation of local vibrational properties of transition metal-containing zeolites modelled by the hybrid QM/MM approach. The scheme is closely related to the vibrational self-consistent field (VSCF) method [1] applied to localised modes in gas-phase molecules [2] but includes the effects of the environments and accurate numerical solutions of the VSCF equations in real space using Schrödinger solvers. The method is implemented as a new module/library of the Py-Chemshell (PCS) [3]. These zeolites calculations are integrated with the QM/MM models, which have been reported previously [4] and implemented in PCS. This work continues the recent extension in PCS to support calculations of infrared and Raman intensities based on the harmonic approximations. Our VSCF implementation tackles anharmonicities, coupling and delocalisation of normal modes - these are recognised as key challenges for first-principle computations of vibrational properties. We carry on from the VSCF calculations to talk about how to assign band positions for intermediates of a selective catalytic redox cycle over Cu-SSZ-13 [5]. We conclude with a discussion of further intermodal correlation and evolution effects.



References:

[1] J. Bowman, et al. Molecular Physics. 2010, 106, 16-18, 2145-2182.
[2] X. Cheng, et al. J. Chem. Phys. 2014, 141.10, 104105.
[3] L. You, et al. J. Chem. Theory Comput. 2019, 1317-1328.
[4] P. Sherwood, et al. Faraday Discuss. 1997, 106, 79-92.
[5] A. Greenaway, et al. Chem. Sci. 2020, 447-455.

I will review recent progress on the development and application of advanced atomistic algorithms to simulate chemomechanical systems where local chemistry and long-range stress are tightly coupled, e.g. at the tip of a propagating crack or the core of a dislocation. I will discuss two general approaches (i): hybrid quantum/classical approaches where bond-breaking is treated at the DFT level embedded within a large-scale classical atomistic model to capture elastic relaxation, including recent applications to tungsten [1] and polymer nanocomposites [2]; (ii) atomistic-to-continuum modelling for crack propagation, either through extracing effective continuum models from large-scale atomistic simulations of fracture [3], or by flexibly embedding an atomistic domain within a continuum model, using a new algorithm to compute bifurcation diagrams for fracture systems [4].

[1] P. Grigorev, T. Swinburne and J. R. Kermode, QM/MM study of hydrogen-decorated screw dislocations in tungsten: ultrafast pipe diffusion, core reconstruction and effects on glide mechanism,
Phys. Rev. Materials 4 023601 (2020)
[2] J. R. Golebiowski, J. R. Kermode, P. D. Haynes and A. A. Mostofi, Atomistic QM/MM simulations of the strength of covalent interfaces in carbon nanotube-polymer composites, Phys. Chem. Chem. Phys 22 12007 (2020)
[3] S. M. Khosrownejad, J. R. Kermode, and L. Pastewka, Quantitative Prediction of the Fracture Toughness of Amorphous Carbon from Atomic-Scale Simulations, Phys. Rev. Materials 5, 023602 (2021)
[4] M. Buze and J. R. Kermode, Numerical-Continuation-Enhanced Flexible Boundary Condition Scheme Applied to Mode-I and Mode-III Fracture, Phys. Rev. E 103, 033002 (2021)

This talk discusses recent rapid solidification experiments on thin film aluminum samples that reveal the presence of lattice orientation gradients within crystallizing grains. To study this phenomenon, a new phase-field crystal (PFC) model that captures the properties of solid, liquid, and vapor phases is proposed to model pure aluminum quantitatively. A coarse-grained representation of this model is used to study rapid solidification in samples approaching micrometer scales. We present recent results of this model that reproduce the experimentally observed orientation gradients within crystallizing grains grown at experimentally relevant rapid quenches, discussing a causal connection between defect formation and orientation gradients.

Using 2D model glasses as a theoretical base for studying amorphous materials, the Local Yield Stress (LYS) method has ranked highly among available predictors for determining the location of plastic events. By probing local regions of atoms, the LYS method quantifies the amount of incremental stress necessary to induce plastic rearrangements. In 2D, the computed local yield stress (Δτ_c) is dependent on the direction in which the stress was applied on the local region. Thus, a more accurate prediction is made by sampling multiple loading directions in each local region. In 3D, in order to fully quantify the local yield surface, we sample 432 distinct loading orientations. The resulting yield surface can be used to calculate the projected local yield stress (Δτ_y) with respect to the loading on the boundary. We then quantify the difference in predictive capabilities of the LYS method when using the full local yield surface, compared to predictions made using only a single local probing direction aligned with the loading on the boundary.

The quasicontinuum (QC) method is a multiscale method that extends the length and time scales accessible to fully-atomistic simulations by orders of magnitude. This talk will describe the development over recent years of a high-performance implementation of QC for three-dimensional (3D) polycrystalline materials of arbitrary multilattice crystal structures. Features include finite temperature support (hot-QC), temporal acceleration via hyperdynamics (hyper-QC), and support for the KIM API (allowing usage of interatomic potentials archived on https://openkim.org). The code is optimized for a heterogeneous computing infrastructure including a mix of shared memory (OpenMP) and distributed memory (MPI) parallelism. A graph partitioning load balancing algorithm has been developed for efficient resource usage. State-of-the art software engineering techniques are applied including Git revision control, test-driven development (TDD), and continuous integration (CI) to minimize errors in this complex software platform (120K lines of code). The code features, structure and performance will be discussed as well as example applications to fracture, nanoindentation and friction in semiconductors, metals and thermal shielding materials.

Simulating athermal fiber-networks involves the computation of cohesive energies between fiber segments. A general fiber in three-dimensional space can be represented as a connected sequence of straight-line segments or as an assembly of Bezier curves. The former can result in poor representation of fiber curvature and has been addressed in previous work. The latter is treated in this work. The computation of energy (and force) between two Bezier curves can be implemented using a distance distribution function. This distribution function lacks a closed form expression and is computed using a mix of contouring techniques, interval trees and quadtrees. This method reduces some of the artificial friction associated with the discretization of continuous fibers.

Non-equilibrium material behavior is still riddled with many open questions: is it possible to predict the far-from-equilibrium material behavior under arbitrary loading from equilibrium data? Is it possible to predict the effect of element transmutation in the mechanical behavior of materials? In this talk, a statistical mechanics framework for non-equilibrium material response will be presented to address these questions.

Polymer-peptide conjugates have been of great interest over the years because of their proposed applications in biotechnology, medicine and nanotechnology. Conjugation of peptides and polymers can result in hybrid materials with increased biological functions like enzymatic action or receptor recognition as well as suppressed surface activity to prevent its degradation by proteolytic enzymes. Thermosensitive polymers like poly(N-isopropylacrylamide) (PNIPAM) that have a lower critical solution temperature (LCST = ~305 K) can undergo a reversible coil-to-globule transition in response to changes in temperature. When conjugated with peptides, PNIPAM has been shown to create temperature-induced phase switches to alter the substrate activity of the peptide, self-assembled mesoglobules etc. However, precise molecular details of these PNIPAM-peptide architectures as well as the effects of PNIPAM chain lengths and peptide type are largely unexplored, especially using computational methods. In this study, we employ coarse grained (CG) molecular dynamics (MD) simulations to simulate PNIPAM-peptide systems using PNIPAM and amino acid models previously developed in the group. We explore a number of peptides formed using combinations of hydrophobic (Leu, Phe, Trp etc.), hydrophilic (Arg, Lys, Asn etc.) and neutral (Gly, Ser, His etc.) amino acids like Gly-Arg-Lys-Phe-Gly. By performing detailed structural analysis, we investigate the molecular-level structure of these conjugates and solvent at the interface.

New coarse-grained (CG) embedded atom method (EAM) potentials have been developed to describe the interatomic interaction for face-centered cubic (FCC) metals viz. Gold (Au), Palladium (Pd), Nickel (Ni), Copper (Cu), Aluminum (Al) and Platinum (Pt). These CG metal models could reproduce several physical, mechanical, and thermodynamic properties obtained from experimental and/or quantum methods. The reliability and robustness of each EAM potential is assessed by performing sensitivity analysis using Sobol sampling. Subsequently, the model was refined with the Bayesian parameter estimation (BPE) and further uncertainty quantification (UQ) using CheKiPEUQ, a Bayesian parameter estimation tool. The refinement and UQ process was accelerated by integrating Gaussian Process (GP) based machine learning (ML) models. The Bayesian refined CG model predicted the metal properties within reasonable confidence levels of experimental data. This improved computational framework can be extended to develop the reliable and accurate CG interatomic potentials for other hard- and soft-materials.

The self-assembly of porous structures including metal-organic frameworks (MOFs) is a complex process. Even though MOF-5 is one of the most widely studied structures in the literature, very little is known about its self-assembly. In this study, the growth of 3-D structure of MOF-5 is successfully investigated using coarse-grained (CG) molecular dynamics (MD) simulations by employing different growth pathways. Interactions between metal nodes, linker, and solvent molecules were developed to reproduce their structure obtained from all-atom MD simulations. These building blocks were self-assembled by employing four different pathways from which only one resulted in a 3-D MOF-5 structure that resembled experimentally reported structure. We will discuss these approaches and our results in detail. This computational approach can be used to study self-assembly of other MOFs to provide insights into their growth and self-assembly process.

Metal nanoparticles (NPs) and their self-assembled functionalized hybrid structures have received huge interest due to their wide application in catalysis, sensors, and optics. In general interactions between these NP hybrids and solvents play an important role in their processing as well as in determining their self-assembled structures. To understand the role of these interactions at the molecular-level, coarse-grained (CG) molecular dynamics (MD) simulations can be employed. However, a key to generate reliable results is to use accurate and robust force-field (FF) parameters of all the components involved. In this study, we will present our methodology and FF parameters to define the interactions between CG models of metals (e.g. gold, silver, copper, etc.) with solvents (water, DMF, etc.) developed in our group. These interactions were developed to reproduce experimental (macroscopic contact angle, interfacial energies, surface tensions) and computational (binding energies) properties.

It is challenging to describe the fluid flow behavior in natural nanoporous materials due to inherent complexities of the porous network topology. In nanoporous medium, collisions and interactions between particles and walls are as important as the collisions and interactions between particles as the molecular mean free path is comparable to the characteristic length scale of the medium. In this work, we employ Many-body Dissipative Particle Dynamics (mDPD) simulations to investigate the fluid flow process through realistic bicontinuous nanoporous media, which are representative models for a wide class of nanoporous materials. Bicontinuous nanoporous models are constructed considering a defined morphology at 0.65 porosity level and varying pore sizes from 30 to 120 nm. The models provide a stochastic description of the morphology and pore size distribution and allow for a direct investigation of the dependence of permeability on average pore-size. Simulation results from fluid obtained by imposing different pressure differences on the nanoporous model by the action of two confining pistons, indicate a linear pressure drop within the nanoporous model, regardless of pore size. The steady state fluid flow through the nanoporous models is proportional to the pressure gradient applied in agreement with the Darcy’s law, in spite of the complexities and different scales of the porous media considered. The predicted pore-size dependence of the permeability is well described by the Hagen-Poiseuille law considering a single shape correction factor that accounts for the flow resistance due to the complex nanoporous morphology.

Shape Failure in amorphous solids arises from local rearrangements in the material's atomic configurations and how these respond to stress. To this end, the Local Yield Stress (LYS) method was recently developed to probe local regions and quantify their susceptibilities to plasticity by measuring the incremental stress required to trigger a local rearrangement. While this methodology shows great potential for enhancing predictive plastic theories, it is not a practical means for characterizing materials structure due to its computationally intensive nature. We propose a manifold learning-based framework to extract microstructural descriptors from atomistic configurations. More specifically, we deploy Diffusion Maps (a nonlinear manifold learning technique) to systematically extract the structural information features from the high-dimensional data (cartesian coordinates of local clusters of atoms and relate this structure to LYS). Due to the nature of the problem, a “point” corresponds to a geometric conformation of atoms, and thus meaningful similarity measure between configurations must be devised for capturing the actual distance between two points. Therefore, we utilize the Grassmann manifold to measure the similarity between the cluster of atoms. We start with a “proof of concept” example crystalline system where our aim is to identify known defects (vacancies and dislocation) using the machine learning approach trained to the geometric conformation of atoms before moving on to amorphous materials to demonstrate its practicability. Finally, our application is a two-dimensional binary glass-forming system, deformed with an incremental Athermal Quasistatic Shear (AQS) method.

Fibrinogen is a plasma-soluble glycoprotein responsible for blood clot formation. After an injury to the blood vessels, fibrinogen molecules are polymerized into fibrin through binding and interactions at their αC domains, forming a fibrin mesh that completes a blood clot. Previous studies indicate that excessive coagulation can lead to harmful diseases and health conditions like strokes and pulmonary embolism. P12, a peptide derived from fibronectin, can mitigate fibrinogen formation and reduce thrombogenesis. In this study, we investigated the properties of the distortion of the αC domain by P12, which inhibits fibrinogen polymerization and prevents the formation of a fiber network. To test this, an in silico model of the αC domain was and potential models of the P12-αC interaction were collected from multiple docking servers. The top models from each server were further evaluated using protein-ligand affinity prediction software. The top scores were all located within residues 518 to 584, which we hypothesized to be the binding site of P12 in the αC domain. To understand the binding interactions between P12 and the αC domain, hydrogen bonding sites were mapped through Cytoscape detection. Six residues displayed hydrogen bonding with the P12 ligand, and VMD analysis indicated that two salt bridges broke due to the change in the structure of the αC domain. Furthermore, the introduction of P12 increased the overall flexibility of the αC domain complex, as indicated by higher RMSF values. The hydrophobicity of the αC domain and P12 complex was mapped using UCSF Chimera to observe the forces that keep the complex intact and determine the degree of hydrophobicity at the P12 binding sites on the αC domain. On both the P12 molecule and the αC domain, the specific sites of interaction were shown to be very hydrophobic. This research suggests that P12 is an effective inhibitor of fibrinogen polymerization due to its interactions with the αC domain.

Atomistic models of materials are always approximate, giving rise to uncertainty. One source is from the cohesive model, due to e.g. the treatment of electronic correlation in density functional theory or the functional form / descriptor basis of an empirical interatomic potential. In many cases, some form of error bound can be derived, and a variety of propagation schemes have been developed.

An equally common, but often unquantified, uncertainty arises when targeting some observable whose mechanism is a priori unknown. As this requires computational exploration, via sampling methods such as acclerated molecular dynamics or saddle search routines, any estimate will be vulnerable to additional sampling data (namely, the discovery of new mechanisms). This complicates how results from atomistic simulations can be confidently used for direct predictions or the parametrization of higher-scale models.

I will discuss some recent efforts to quantify sampling uncertainties when system trajectories can be considered metastable, allowing a rigorous Markov chain / kinetic Monte Carlo representation[1]. Using examples of point defect diffuson in alloys and structural transformations of atomic clusters, I will show how the bounding and propagation of sampling uncertainty depends critically on both the quantity of interest and the sampling method. In particular, errors from dynamic sampling methods can be rigorously bounded[2], whilst static methods have approximate bounds that benefit from a committee approach[3]. Some of these ideas are exploited by the TAMMBER code[4], which optimally manages thousands of simulation tasks in parallel, autonomously distributing (dynamic) sampling effort such that the target uncertainty, propagated through the Markov chain, reduces maximally fast. If time allows I will also discuss how this same framework can be used to propagate and target model form uncertainty in empirical potentials.

[1] TD Swinburne, Comp. Mat. Sci. 2021, 10.1016/j.commatsci.2020.110256
[2] TD Swinburne and D Perez, NPJ Comp. Mat. 2020, 10.1038/s41524-020-00463-8
[3] TD Swinburne and DJ Wales, JCTC 2020, 10.1021/acs.jctc.9b01211
[4] github.com/tomswinburne/tammber

The Bethe-Salpeter equation (BSE) of many-body perturbation theory (MBPT) is the state of-the-art approach to determine optical, UV, and x-ray absorption spectra of crystalline materials. Likewise, time-dependent density functional theory (TDDFT) enables the study of excitations in molecules and solids either via time propagation or by computing their response to external perturbations in the linear regime. Such high-level approaches are, however, compute-intensive, often hampering their application to complex systems. This is in particular so, if highly accurate basis sets are employed to solve the underlying DFT problem and to expand the non-local operators appearing in excited-sate methodology. In this talk, we will discuss current limitations and present various algorithms to overcome some of the bottlenecks. All this is realized in the all-electron linearized augmented planewave + local orbital (LAPW+lo) code exciting [1]. The first example concerns (i) an iterative solver for diagonalizing huge sparse matrices that appear when treating low-dimensional systems in codes that apply periodic boundary conditions. Making use of this method, allows us to reach microHartree precision for the energies of atoms and molecules [2]. Building on this, we are (ii) extending this approach towards the time domain in TDDFT [3]. Moreover, we make use of a (iii) recently developed algorithm for substantially speeding up the solution of the BSE [4], and (iv) demonstrate the successful implementation in exciting by applications to complex systems. Finally, we will discuss how to bring high-level approaches to exascale as envisioned in the European Center of Excellence NOMAD [5].


[1] A. Gulans, S. Kontur, C. Meisenbichler, D. Nabok, P. Pavone, S. Rigamonti, S. Sagmeister, U. Werner, and C. Draxl, exciting: a full-potential all-electron package implementing density-functional theory and many-body perturbation theory, J. Phys: Condens. Matter (Topical Review) 26, 363202 (2014).
[2] A. Gulans, A. Kozhevnikov, and C. Draxl, Microhartree precision in density functional theory calculations, Phys. Rev. B 97, 161105(R) (2018).
[3] R. Rodrigues Pela and C. Draxl, All-electron full-potential implementation of real-time TDDFT in exciting, Electronic Structure (2021); in print.
[4] F. Henneke, L. Lin, C. Vorwerk, C. Draxl, R. Klein, and C. Yang, Fast optical absorption spectra calculations for periodic solid state systems, Commun. Appl. Math. Comp. Sci. 15, 89 (2020).
[5] NoMAD Center of Excellence, https://nomad-coe.eu. NOMAD CoE receives funding from the European Union’s Horizon 2020 research and innovation program under the grant agreement 951786.

We have developed a method, the combines atomistic spin dynamics and ab-initio molecular dynamics (ASD-AIMD) to investigate magnetic materials at high temperature.[1] We apply the method to study the mutual coupling of spin fluctuations and lattice vibrations in paramagnetic CrN and elemental Fe. The two degrees of freedom can be dynamically coupled, and our method allows for observations of nonadiabatic effects. In CrN, those effects is found to suppress the phonon lifetimes at low temperature compared to an adiabatic, fast-magnetism approximation. The dynamic coupling identified here provides an explanation for the experimentally observed unexpected temperature dependence of the thermal conductivity of CrN and other magnetic semiconductors at and above the magnetic ordering temperature.

Our method is based on an alternating scheme of ASD and AIMD steps. In every step of AIMD, constrained non-collinear magnetic moments are used in the density functional theory calculations. The obtained forces are used to update the atomic positions by a single-step propagation of the atoms. With the new atomic positions, new distance dependent magnetic exchange interactions are obtained and used to propagate the magnetic state for the same time-step in the ASD simulation together giving a coupled ASD-AIMD step that are repeated until convergence of studied quantaties.

The dynamics of the CrN spin system impacts atomic lattice vibrations and the quantities related to it, in particular phonon lifetimes. The magnetic state and its dynamics are influenced by the lattice vibrations. Slightly above the transition temperature (300K), the dynamical coupling, is found to significantly reduce the phonon lifetimes of the acoustic modes. In contrast, at high temperatures (1000 K) the dynamical coupling and its impact on phonon lifetimes are largely reduced.

Our ASD-AIMD method provides access to magnetic materials at high temperature with AIMD- level accuracy previously only existing for non-magnetic materials. We demonstrate its reliability by deriving the free energy differences between alpha- gamma- and delta-Fe as a function of temperature and investigate paramagnetic Fe under Earth-core conditions.

[1] Anomalous Phonon Lifetime Shortening in Paramagnetic CrN Caused by Spin-Lattice Coupling: A Combined Spin and Ab Initio Molecular Dynamics Study
Irina Stockem, Anders Bergman, Albert Glensk, Tilmann Hickel, Fritz Körmann, Blazej Grabowski, Jörg Neugebauer, and Björn Alling
Physical Review Letters 121, 125902 (2018)

When 2D materials are stacked to form a van der Waals (vdW) heterostructure with a small angular and/or lattice mismatch, a moiré superlattice may appear, which could lead to novel quantum phenomena, driven by an interplay of flat energy bands and strong electron-electron interactions. In particular, the moiré superlattices comprised of 2D semiconductors, such as transition metal dichalcogenides (TMDs), could host localized, long-lived and valley-polarized moiré excitons, which are envisioned as single-photon emitters in quantum information and optoelectronic applications. In this talk, we provide an overview of our recent effort to elucidate the properties of moiré excitons in twisted 2D vdW heterostructures from first principles. We predict that moiré excitons can be formed in twisted 2D organic–inorganic halide perovskites, group IV monochalcogenides, in addition to the well-known TMD bilayers. In conjunction with localized moiré excitons, nearly flat bands are also formed in these vdW heterostructures. 1D flat bands, anisotropic and 1D moiré excitons are also predicted to emerge in some of the vdW heterostructures, specially, when defects are present. We explore how vertical electric field can be tuned to control the position, polarity, emission energy, and hybridization strength of the moiré excitons. In the second part of the talk, we will shed light on exciton dynamics in TMD vdW heterostructures. In particular, we unravel the competition of charge and energy transfer in type-II WS₂/MoSe₂ heterostructures. Incorporating excitonic effect in non-adiabatic electronic dynamics, our first-principles study uncovers a two-step process in ultrafast energy and charge transfer, unravel their relative efficiencies and explore the means to control their competition. The excitonic effect is revealed to drive the efficient energy and charge transfer in the heterostructures.