Large Simulation Cells Research Articles

Integrating Newton's equations of motion in the reciprocal space.

We here present the normal dynamics technique, which recasts the Newton's equations of motion in terms of phonon normal modes by exploiting a proper sampling of the reciprocal space. After introducing the theoretical background, we discuss how the reciprocal space sampling enables us to (i) obtain a computational speedup by selecting which and how many wave vectors of the Brillouin zone will be considered and (ii) account for distortions realized across large atomic distances without the use of large simulation cells. We implemented the approach into an open-source code, which we used to present three case studies: in the first one, we elucidate the general strategy for the sampling of the reciprocal space; in the second one, we illustrate the potential of the approach by studying the stabilization effect of temperature in α-uranium; and in the last one, we investigate the characterization of Raman spectra at different temperatures in MoS2/MX2 transition metal dichalcogenide heterostructures. Finally, we discuss how the procedure is general and can be used to simulate periodic, semiperiodic, and finite systems such as crystals, slabs, nanoclusters, or molecules.

Spatially Resolved Uncertainties for Machine Learning Potentials.

Machine learning potentials have become an essential tool for atomistic simulations, yielding results close to ab initio simulations at a fraction of computational cost. With recent improvements on the achievable accuracies, the focus has now shifted on the data set composition itself. The reliable identification of erroneously predicted configurations to extend a given data set is therefore of high priority. Yet, uncertainty estimation techniques have achieved mixed results for machine learning potentials. Consequently, a general and versatile method to correlate energy or atomic force uncertainties with the model error has remained elusive to date. In the current work, we show that epistemic uncertainty cannot correlate with model error by definition but can be aggregated over groups of atoms to yield a strong correlation. We demonstrate that our method correctly estimates prediction errors both globally per structure and locally resolved per atom. The direct correlation of local uncertainty and local error is used to design an active learning framework based on identifying local subregions of a large simulation cell and performing ab initio calculations only for the subregion subsequently. We successfully utilized this method to perform active learning in the low-data regime for liquid water.

Atomistic modeling of idealized equal channel angular pressing process

Severe plastic deformation (SPD) processes are forming techniques that impose large plastic strains to achieve ultra-fine grained and nanocrystalline microstructures. The refined microstructure results in a significant increase in strength and in many cases without a considerable loss in ductility. Despite their long existence, a comprehensive understanding of the deformation and grain refinement mechanisms during SPD is still lacking. Atomistic simulations that are regularly used to obtain insights into material behavior have rarely been used to study SPD processes. In this work, we use large scale atomistic simulations of the molecular statics/dynamics kind to model the equal channel angular pressing (ECAP) process, a specific SPD process. Three different materials—Al, Ni and Cu—are used to model high, medium and low stacking fault energy materials. Large simulation cells are used to facilitate grain refinement at the nanoscale. The simulations are carefully analyzed in terms of stress–strain behavior, dislocation activity and grain refinement. All samples show the formation of deformation twins at the nanoscale, which generally detwin at increased strains via dislocation twin interactions. A significant number of stair-rod dislocations form via dislocation–dislocation interactions and are present at low angle grain boundaries in the early stages of deformation. Such stair-rod dislocations also result in the formation of stacking fault tetrahedra which are present in significant numbers in all samples. Finally, difficulties in studying grain refinement at very large strains in atomistic simulations are discussed.

Magnetic ground state of holmium nitride

Rare-earth mononitrides such as HoN exhibit a wide range of useful properties leading to potential applications as magnetic semiconductors, spintronic half-metals, or magnetocaloric refrigerants in hydrogen liquefaction systems. First-principle calculations of electronic structures and related properties of such materials should correctly reproduce their magnetic moment. First, we identify the unusually high number of unoccupied electronic states which guarantees that the energy minimum identified is the global one. Second, we develop a method which allows the experimentally relevant magnetisation to constitute an energy minimum, emphasising the favourable distribution of the spins in an exceptionally large simulation cell. Third, we examine the dependence of selected HoN characteristics on the cell size and on the magnetisation. The results provide a theoretical insight into the spin structure of rare-earth nitrides and allow one to use the correct methodology of similar calculations of properties of strongly correlated materials.

Enabling automated high-throughput Density Functional Theory studies of amorphous material surface reactions

Amorphous materials are increasingly considered for electronic and energy applications. However, modeling site-specific processes on amorphous materials is prohibitively expensive due to the large configuration space, and the large simulation cell needed to properly capture material properties. In this work, we develop a high-throughput workflow, powered by semiempirical methods in combination with Density Functional Theory, to simulate amorphous surfaces and utilize statistical sampling methods to assess surface reaction kinetics. We employ the Smooth Overlap of Atomic Positions (SOAP) to featurize surface species and apply a Bayesian Gaussian Mixture clustering model to perform a site reduction analysis. Finally, we automate the generation of input images for Improved-Tangent Nudged Elastic Band (IT-NEB) to simulate an etching reaction on the reduced sites. The resulting etching barriers are found to follow experimental etching trends for both amorphous silicon and amorphous carbon. Notably, our method reproduces etching barriers using a significantly reduced amount of sites (6 sites) out of all potential sites on the amorphous material surfaces. The successful site reduction method described in this study opens amorphous materials to high-throughput explorations of interfacial chemistry and surface properties.

Simulated migration behavior of metastable Σ3 (11 8 5) incoherent twin grain boundaries

Molecular dynamics simulations are used to examine the migration behavior of one incoherent twin, a Σ3 (11 8 5) / (8 11 5) grain boundary. The boundary is known to exhibit non-Arrhenius boundary migration that slows as temperature increases. This behavior is examined in 165 metastable structures of the same boundary and in a large simulation cell where smaller length scales can have less of an effect. The metastable boundaries show diverse migration behaviors from non-Arrhenius to Arrhenius, though the majority of them exhibit non-Arrhenius behavior. The large simulation cell sizes show no dependence on system size, eliminating concerns about the effect of facet nucleation on the migration in periodic simulation cells. However, facet structures play an important role in the migration of the metastable boundaries. Boundaries with larger facets typically migrate faster than those with smaller facets and with defects in the boundary structure.

Missing theoretical evidence for conventional room-temperature superconductivity in low-enthalpy structures of carbonaceous sulfur hydrides

To elucidate the geometric structure of the putative room-temperature superconductor, carbonaceous sulfur hydride (C-S-H), at high pressure, we present the results of an extensive computational structure search of bulk C-S-H at 250 GPa. Using the minima hopping structure prediction method coupled to the GPU-accelerated sirius library, more than 17 000 local minima with different stoichiometries in large simulation cells were investigated. Only 24 stoichiometries are favorable against elemental decomposition, and all of them are carbon-doped ${\mathrm{H}}_{3}\mathrm{S}$ crystals. The absence of van Hove singularities or similar peaks in the electronic density of states of more than 3000 candidate phases rules out conventional superconductivity in C-S-H at room temperature.

The phase behaviour of cetyltrimethylammonium chloride surfactant aqueous solutions at high concentrations: an all-atom molecular dynamics simulation study.

We explore the phase behaviour of aqueous solutions of the cetyltrimethyl ammonium chloride (CTAC) surfactant and in particular the transition from the micellar phase (L1) to the hexagonal columnar phase (H1) by employing all-atom (AA) molecular dynamics (MD) simulations for six CTAC concentrations in the range of 34.1 wt% to 70.5 wt%, at the temperature of 318 K and pressure of 1 atm. For the concentrations considered, we examine the spontaneous occurrence of the H1 phase by testing a number of plausible values for the linear density (molecules per unit length) along the cylindrical columns. Using large simulation cells and starting from random initial configurations, the MD simulations demonstrate that the micellar phase occurs for concentrations up to 50.0 wt%, with CTAC molecules self-assembling into a mixture of spherical and rod-like micelles. At even higher concentrations, the system self-organizes into the H1 phase in accordance with the available experimental data. For the analysis of the MD trajectories, we devise a clustering algorithm based on Voronoi tesselation which enables (a) the thorough characterization of the shape and structure of both molecules and assemblies, and (b) the investigation of the positional and orientational order in the system that are further scrutinised using radial pair correlation functions and X-ray diffraction patterns. Our work paves the way for the investigation of the phase behaviour at high concentrations of other surfactants.

Using Neural Network Force Fields to Ascertain the Quality of Ab Initio Simulations of Liquid Water.

Accurately simulating the properties of bulk water, despite the apparent simplicity of the molecule, is still a challenge. In order to fully understand and reproduce its complex phase diagram, it is necessary to perform simulations at the ab initio level, including quantum mechanical effects both for electrons and nuclei. This comes at a high computational cost, given that the structural and dynamical properties tend to require long timescales and large simulation cells. In this work, we evaluate the errors that density functional theory (DFT)-based simulations routinely incur into due time- and size-scale limitations. These errors are evaluated using neural-network-trained force fields that are accurate at the level of DFT methods. We compare different exchange and correlation potentials for properties of bulk water that require large timescales. We show that structural properties are less dependent on the system size and that dynamical properties such as the diffusion coefficient have a strong dependence on the simulation size and timescale. Our results facilitate comparisons of DFT-based simulation results with experiments and offer a path to discriminate between model and convergence errors in these simulations.

Screened configuration interaction method for open-shell excited states applied to NV centers

We present a computational approach, based on density functional theory and screened configuration interaction, able to accurately treat highly correlated electron spins localized around semiconductor defects, typically occurring in the context of qubit implementations. The method is computationally not more demanding than a usual density functional theory calculation, which makes it suitable for the calculation of isolated defects, or defect complexes, typically requiring large simulation cells. We illustrate the approach by applying it to the three different charge states of the nitrogen vacancy defect in diamond and obtain very good agreement with experiment.

Screw dislocation impingement and slip transfer at fcc-bcc semicoherent interfaces

Understanding of dislocation-interface interactions at the atomistic level could provide useful insights into the mechanisms behind slip transfer, a key information for higher length-scale formulations. A systematic study of screw dislocation impingement on both Fe/Al and Fe/Cu interfaces with Kurdjumov-Sachs orientation relationship has been carried out using atomistic simulations. For the “strong” Fe/Al interface, there is appreciable dislocation delocalization, in which the change in dislocation content provided by Nye tensor analysis. Interestingly, the delocalization is significant enough that large simulation cells must be employed in order to avoid boundary effects. When the impinging dislocations are far away from the boundaries, the large accumulation of shear strains nearby the impinging dislocations at the interface causes dislocation nucleation/emission at these places. While for the “weak” Fe/Cu interface, the dislocation delocalization at the interface is notably long-range, extending more than 20 nm away from where the initial dislocations entered the interface.

How Effective are Retro-Inverso Peptides? Insights from MD Supported by Paramagnetic NMR Data

Early MD simulations were largely limited to observation of fast structural fluctuations in protein molecules. With recent advent of GPU computing, it became possible to model far more relevant dynamic processes: protein folding, ligand binding, allostery, etc. In our study, we use long unbiased MD simulations to model binding of the so-called retro-inverso (ri) peptide to its protein target. Retro-inverso peptides are designed from native peptides by replacing L-amino acids with D-amino acids and inverting the amino-acid sequence. They are believed to retain the binding affinity of their parent peptides, while being resistant to proteolytic degradation. We investigated the binding of retro-inverso peptide ri-Sos1 (rrrppvpppv) to N-terminal SH3 domain from adaptor protein Grb2. The peptide was initially placed in a random position inside a large simulation cell. During the course of the simulation, we observed the binding of ri-Sos1 at the canonical SH3 binding site. The resulting complex features a number of key hydrogen bonds that are similar to the ones found in the native Sos1/Grb2 N-SH3 structure. However, the pattern of these interactions has undergone some noticeable rearrangements. We argue that the concept of retro-inversion lacks any solid theoretical justification and, therefore, ri-peptides are only poor imitators of their parent peptides. To validate our findings, we conducted a number of NMR experiments. In particular, we used the peptides crrrppvpppv and rrrppvpppvc tagged with paramagnetic label MTSL to measure paramagnetic relaxation enhancements (PRE) in Grb2 N-SH3. The results were compared with the PRE rates rigorously calculated from the MD trajectories involving the two respective MTSL-tagged peptides. The 1HN,15N-HSQC-based titration data showed that ri-Sos1 is a much weaker binder than the original Sos1, thus confirming the tenuous nature of ri-peptide ligands. We acknowledge grant support from SPbU 51142660.

Effect of non-Heisenberg magnetic interactions on defects in ferromagnetic iron

Fundamental flaws in the Heisenberg Hamiltonian are highlighted in the context of its application to bcc Fe, including the particular issues arising when modeling lattice defects using spin dynamics. Exchange integrals are evaluated using the magnetic force theorem. The bilinear exchange coupling constants are calculated for all the interacting pairs of atomic magnetic moments in large simulation cells containing defects, enabling a direct mapping of the magnetic energy onto the Heisenberg Hamiltonian and revealing its limitations. We provide a simple procedure for extracting the Landau parameters from density-functional theory calculations to construct a Heisenberg-Landau Hamiltonian. We quantitatively show how the Landau terms correct the exchange-energy hypersurface, which is essential for the accurate evaluation of energies and migration barriers of defects.

Aggregate-State Effects in the Atomistic Modeling of Organic Materials for Electrochemical Energy Conversion and Storage Devices: A Perspective.

Development of new functional materials for novel energy conversion and storage technologies is often assisted by ab initio modeling. Specifically, for organic materials, such as electron and hole transport materials for perovskite solar cells, LED (light emitting diodes) emitters for organic LEDs (OLEDs), and active electrode materials for organic batteries, such modeling is often done at the molecular level. Modeling of aggregate-state effects is onerous, as packing may not be known or large simulation cells may be required for amorphous materials. Yet aggregate-state effects are essential to estimate charge transport rates, and they may also have substantial effects on redox potentials (voltages) and optical properties. This paper summarizes recent studies by the author’s group of aggregation effects on the electronic properties of organic materials used in optoelectronic devices and in organic batteries. We show that in some cases it is possible to understand the mechanism and predict specific performance characteristics based on simple molecular models, while in other cases the inclusion of effects of aggregation is essential. For example, it is possible to understand the mechanism and predict the overall shape of the voltage-capacity curve for insertion-type organic battery materials, but not the absolute voltage. On the other hand, oligomeric models of p-type organic electrode materials can allow for relatively reliable estimates of voltages. Inclusion of aggregate state modeling is critically important for estimating charge transport rates in materials and interfaces used in optoelectronic devices or when intermolecular charge transfer bands are important. We highlight the use of the semi-empirical DFTB (density functional tight binding) method to simplify such calculations.

(Invited) Dynamics of Nitrogen Functionalization of Graphene

For graphene to be used efficiently in scientific and technological applications, its electronic properties need to be tailored to suit precise demands. The formation of a forbidden energy gap or the modification of the Fermi energy is achievable with functionalization techniques that simultaneously leave most of the atomic structure of graphene unchanged. Nitrogen functionalization is one of the most promising and studied types of treatment, in large part due to the proximity between carbon and nitrogen in the periodic table. It has already produced promising experimental results in the fields of high-sensibility biodetectors, field-effect transistors, solar cells, and supercapacitors, among others.To further understand this functionalization, we have performed an ab initio electronic structure study of different doping configurations and functionalization dynamics of nitrogen treatments of graphene. In particular, we have study the effect of simulation cell size used to carry our these calculations. We found that the inclusion of the special point K of the first Brillouin Zone in the momentum space sampling lowers formation energy results and should not be overlooked. Nudged Elastic Band calculations were also completed to study different incorporation mechanisms from an adsorbed state to in-plane doping. In the presence of native defects, low barriers of 0.55 eV and 0.46 eV were obtained. When no defects are present, much larger barriers between 3.70 eV to 4.38 eV were found, which suggests an external source of energy is required to complete the incorporation.In the second part of this presentation, we will report our result of using machine learning models to predict higher-level quantities. The goal of this work is to study the dependence of the Raman response of nitrogen functionalization of graphene. To simulate such a response, very large simulation cells need to be considered that would be out of reach for ab initio techniques, hence the use of machine learning methods to assist in the evaluation of the properties. Machine learning methods are now used more and more as a substitute for density functional theory calculations due to their low computational costs. However, in more advanced cases, relevant datasets are not always available, and the effort that would be necessary to generate the needed data suppress the advantages of using machine learning to speed up the calculations. Furthermore, the process of training a reliable model is not trivial and can also be expensive which should encourage the usage and diffusion of portable trained models. Figure 1

Role of hemibonding in the structure and ultraviolet spectroscopy of the aqueous hydroxyl radical.

The presence of a hemibond in the local solvation structure of the aqueous hydroxyl radical has long been debated, as its appearance in ab initio simulations based on density functional theory is sensitive to self-interaction error (favoring a two-center, three-electron hemibond) but also to finite-size effects. Simulations reported here use a mixed quantum mechanics/molecular mechanics (QM/MM) framework in a very large periodic simulation cell, in order to avoid finite-size artifacts and to facilitate testing of various density functionals, in order to probe the effects of delocalization error. The preponderance of hemibonded structures predicted by generalized gradient approximations persists in simulations using the hybrid functionals B3LYP and PBE0, but is reduced to a minor population if the fraction of exact exchange is increased to 50%. The hemibonded population is also small in simulations employing the long-range corrected functional LRC-ωPBE. Electronic spectra are computed using time-dependent density functional theory, and from these calculations emerges a consensus picture in which hemibonded configurations play an outsized role in the absorption spectrum, even when present as a minority species. An intense 1b2(H2O) → 2pπ(˙OH) charge-transfer transition in hemibonded configurations of the radical proves to be responsible for an absorption feature at 230 nm that is strongly shifted with respect to the gas-phase absorption at 307 nm, but this intense feature is substantially diminished in aqueous geometries where the hemibond is absent. Although not yet sufficient to quantitatively establish the population of hemibonded ˙OH(aq), these simulations do suggest that its presence is revealed by the strongly shifted ultraviolet absorption spectrum of the aqueous radical.

Harnessing machine learning potentials to understand the functional properties of phase-change materials

Abstract

A Ligand Field Molecular Mechanics Study of CO2‐Induced Breathing in the Metal–Organic Framework DUT‐8(Ni)

AbstractFlexible metal–organic Frameworks (MOFs) are an interesting class of materials due to their diverse properties. One representative of this class is the layered‐pillar MOF DUT‐8(Ni). This MOF consists of Ni2 paddle wheels interconnected by naphthalene dicarboxylate linkers and dabco pillars (Ni2(ndc)2(dabco), ndc = 2,6‐naphthalene–dicarboxylate, dabco = 1,4‐diazabicyclo‐[2.2.2]‐octane). DUT‐8(Ni) undergoes a volume change of over 140% upon adsorption of guest molecules. Herein, a ligand field molecular mechanics (LFMM) study of the CO2‐induced flexibility of DUT‐8(Ni) is presented. LFMM is able to reproduce experimental and DFT structural features as well as properties that require large simulation cells. It is shown that the transformation energy from a closed to open state of the MOF is overcompensated fivefold by the host–guest interactions. Structural characteristics of the MOF explain the shape of the energy profile at different loading states and provide useful insights to the interpretation of previous experimental results.

Computation of NMR Shifts for Paramagnetic Solids Including Zero-Field-Splitting and Beyond-DFT Approaches. Application to LiMPO4 (M = Mn, Fe, Co, Ni) and MPO4 (M = Fe, Co)

A protocol for the computation of NMR shifts of paramagnetic spin-coupled solids is described and applied to the olivine-type lithium-ion battery cathode materials LiMPO4 (M = Mn, Fe, Co, Ni) and to the related binary phosphates MPO4 (M = Fe, Co). The computational approach includes Fermi-contact (FC), pseudocontact (PC), as well as orbital contributions to the shifts and combines periodic density functional computations with multireference post-Hartree–Fock methods. A shift formalism based on EPR spin-Hamiltonian parameters corrected for residual spin couplings within the Curie–Weiss regime is used. Orbital shieldings as well as hyperfine couplings (HFCs) are obtained at DFT levels for large simulation cells using the CP2K code, taking advantage of hybrid functionals for the HFCs. g-Tensors and zero-field-splitting tensors are computed within an incremental cluster model, using CASSCF and NEVPT2 approaches, providing thus the first post-Hartree–Fock treatments of these properties for extended solids. Com...

Reducing uncertainty in simulation estimates of the surface tension through a two-scale finite-size analysis: thicker is better.

Recent simulation studies of the surface tension γ, and other properties of thin free-standing films, have revealed unexpected finite size effects in which the variance of the properties vary monotonically with the in-plane width of the films, complicating the extrapolation of estimates of film properties to the thermodynamic limit. We carried out molecular dynamics simulations to determine the origin of this phenomenon, and to address the practical problem of developing a more reliable methodology for estimating γ in the thermodynamic limit. We find that there are two distinct finite size effects that must be addressed in a finite size analysis of γ in thin films. The first finite size scale is the in-plane width of the films and the second scale is the simulation cell size in the transverse direction. Increasing the first scale enhances fluctuations in γ, measured by the standard deviation of their distribution, while increasing the second reduces γ fluctuations due to a corresponding increased ‘freedom’ of the film to fluctuate out of plane. We find that using progressively large simulation cells in the transverse direction, while keeping the film width fixed to an extent in which the full bulk liquid zone is developed, allows us to obtain a smooth extrapolation to the thermodynamic limit, enabling a reduction of the γ uncertainty to a magnitude on the order of 1% for systems having a reasonably large size, i.e., O (1 μm).