Method For Large-scale Simulation Research Articles

Efficient all-electron hybrid density functionals for atomistic simulations beyond 10 000 atoms.

Hybrid density functional approximations (DFAs) offer compelling accuracy for abinitio electronic-structure simulations of molecules, nanosystems, and bulk materials, addressing some deficiencies of computationally cheaper, frequently used semilocal DFAs. However, the computational bottleneck of hybrid DFAs is the evaluation of the non-local exact exchange contribution, which is the limiting factor for the application of the method for large-scale simulations. In this work, we present a drastically optimized resolution-of-identity-based real-space implementation of the exact exchange evaluation for both non-periodic and periodic boundary conditions in the all-electron code FHI-aims, targeting high-performance central processing unit (CPU) compute clusters. The introduction of several new refined message passing interface (MPI) parallelization layers and shared memory arrays according to the MPI-3 standard were the key components of the optimization. We demonstrate significant improvements of memory and performance efficiency, scalability, and workload distribution, extending the reach of hybrid DFAs to simulation sizes beyond ten thousand atoms. In addition, we also compare the runtime performance of the PBE, HSE06, and PBE0 functionals. As a necessary byproduct of this work, other code parts in FHI-aims have been optimized as well, e.g., the computation of the Hartree potential and the evaluation of the force and stress components. We benchmark the performance and scaling of the hybrid DFA-based simulations for a broad range of chemical systems, including hybrid organic-inorganic perovskites, organic crystals, and ice crystals with up to 30 576 atoms (101 920 electrons described by 244 608 basis functions).

Capture efficiencies of point defects by 1/2 〈111〉 edge dislocations in tungsten using atomistic simulations

We systematically study the interaction behaviors between a point defect (PD) and two types of 1/2 〈111〉 {110} and 1/2 〈111〉 {112} edge dislocations using molecular dynamics (MD), elastic dipole tensor, and kinetic Monte Carlo (KMC) methods in tungsten (W). We first determine the capture region of a vacancy and four types of self-interstitial atoms (SIAs) by calculating their binding energies to both edge dislocations using MD. We then find anisotropic distributions of the diffusion energy barriers of PDs in dislocation systems employing the MD and elastic dipole tensor methods. Subsequently, we utilize the KMC method to obtain the lifetime of PDs and calculate the capture efficiencies of both edge dislocations to PDs at various temperatures and dislocation densities. Results show the capture efficiency of a vacancy decreases with increasing temperature, while it first increases and then decreases with increasing temperature for SIAs. As dislocation density increases, the capture efficiency of PDs increases. Ultimately, the swelling rate is estimated based on the capture efficiency and falls within the range of experimental data. The swelling rate increases first and then decreases with increasing temperature, and the maximum swelling rate occurs around 1200 K, while it increases with increasing the dislocation density. The present results are critical for understanding the interaction behaviors between PDs and 1/2 〈111〉 edge dislocations in W, which provide an essential database for large-scale simulation methods such as rate theory and phase-field simulations.

Adaptive Subsystem Density Functional Theory.

Subsystem density functional theory (DFT) is emerging as a powerful electronic structure method for large-scale simulations of molecular condensed phases and interfaces. Key to its computational efficiency is the use of approximate nonadditive noninteracting kinetic energy functionals. Unfortunately, currently available nonadditive functionals lead to inaccurate results when the subsystems interact strongly such as when they engage in chemical reactions. This work disrupts the status quo by devising a workflow that extends subsystem DFT's applicability also to strongly interacting subsystems. This is achieved by implementing a fully automated adaptive definition of subsystems which is realized during geometry optimizations or ab initio molecular dynamics simulations. The new method prescribes subsystem merging and splitting events redistributing the resources (both for work and data) in an efficient way making use of modern parallelization strategies and object-oriented programming. We showcase the method with examples probing from moderate-to-strong inter-subsystem interactions, opening the door to using subsystem DFT for modeling chemical reactions in molecular condensed phases with a black box computational tool.

Impact of hydrogen coverage on silane adsorption during Si epitaxy from ab initio simulations

Epitaxy by Chemical Vapor Deposition (CVD) is a commonly used technique for the growth of Si alloys in microelectronic devices. In this work, we perform Density Functional Theory (DFT) simulations to study the impact of H coverage on the dissociative adsorption of silane on a Si(001) surface. Silane adsorption is systematically found thermodynamically favorable but its kinetics is limited by H2 desorption which exhibits an activation energy of 2.4 eV. In addition, our calculations suggest that hydrogenated surfaces tend to reduce the adsorption activation energies compared to uncovered surfaces. This work provides an atomistic description of the SiH4 adsorption mechanisms and associated energies for the modeling of the epitaxial deposition process using large scale simulation methods.

Accurate force field of two-dimensional ferroelectrics from deep learning

The discovery of two-dimensional (2D) ferroelectrics with switchable out-of-plane polarization such as monolayer $\alpha$-In$_2$Se$_3$ offers a new avenue for ultrathin high-density ferroelectric-based nanoelectronics such as ferroelectric field effect transistors and memristors. The functionality of ferroelectrics depends critically on the dynamics of polarization switching in response to an external electric/stress field. Unlike the switching dynamics in bulk ferroelectrics that have been extensively studied, the mechanisms and dynamics of polarization switching in 2D remain largely unexplored. Molecular dynamics (MD) using classical force fields is a reliable and efficient method for large-scale simulations of dynamical processes with atomic resolution. Here we developed a deep neural network-based force field of monolayer In$_2$Se$_3$ using a concurrent learning procedure that efficiently updates the first-principles-based training database. The model potential has accuracy comparable with density functional theory (DFT), capable of predicting a range of thermodynamic properties of In$_2$Se$_3$ polymorphs and lattice dynamics of ferroelectric In$_2$Se$_3$. Pertinent to the switching dynamics, the model potential also reproduces the DFT kinetic pathways of polarization reversal and 180$^\circ$ domain wall motions. Moreover, isobaric-isothermal ensemble MD simulations predict a temperature-driven $\alpha \rightarrow \beta$ phase transition at the single-layer limit, as revealed by both local atomic displacement and Steinhardt's bond orientational order parameter $Q_4$. Our work paves the way for further research on the dynamics of ferroelectric $\alpha$-In$_2$Se$_3$ and related systems.

An GPU-accelerated particle tracking method for Eulerian–Lagrangian simulations using hardware ray tracing cores

To address the high computational cost of particle tracking for realistic Eulerian–Lagrangian simulations, a novel efficient and robust particle tracking method (RT method) for unstructured meshes is presented. The method, for the first time, leverages both hardware ray tracing (RT) cores and GPU parallel computing technology to accelerate Eulerian–Lagrangian simulations. The method includes a hardware-accelerated hosting cell locator using bounding volume hierarchy tree (BVH) and a robust treatment of particle-wall interaction (multiple specular reflection) using an improved neighbor searching approach. The method is implemented in a GPU-accelerated open-source code, which is verified against a reference neighbor-searching particle-tracking method (NS method) and experimental observations. To evaluate the performance of our method, several numerical simulations of fluid-driven scalar transport problem are solved. Using a verification case, we show that the particle distribution simulated by our code is in a good agreement with an experimental observation. Tracking failures and stuck particles are not observed in any simulations. Benchmark results indicate that our RT method leads to a roughly 1.8−2.0× performance improvement compared to the reference NS method for large-scale simulations (millions of mesh cells and particles).

Fluid Flow Characterization Framework for Naturally Fractured Reservoirs Using Small-Scale Fully Explicit Models

Flow modelling challenges in fractured reservoirs have led to the development of many simulation methods. It is often unclear which method should be employed. High-resolution discrete fracture and matrix (DFM) studies on small-scale representative models allow us to identify dominant physical processes influencing flow. We propose a workflow that utilizes DFM studies to characterize subsurface flow dynamics. The improved understanding facilitates the selection of an appropriate method for large-scale simulations. Validation of the workflow was performed via application on a gas reservoir represented using an embedded discrete fracture model, followed by the comparison of results obtained from hybrid and dual-porosity representations against fully explicit simulations. The comparisons ascertain that the high-resolution small-scale DFM studies lead to a more accurate upscaled model for full field simulations. Additionally, we find that hybrid implicit–explicit representations of fractures generally outperform pure continuum-based models.

A level-set method for large-scale simulations of three-dimensional flows with moving contact lines

A new method is presented to perform three-dimensional simulations of two-phase flows with moving contact lines using level-set. To account for the full range of length scales involved in the physical problem under realistic conditions, without having to resolve the flow down to the smallest continuum scale, a dynamic contact angle model based on asymptotic theory is used in conjunction with the computational method. Contact-angle hysteresis is also represented in this methodology. The method is validated against simulations wherein the flow is fully resolved over all length scales, and experiments of spreading droplets and droplets sliding down an inclined substrate.

Implementing discrete element method for large-scale simulation of particles on multiple GPUs

Recently, the use of graphics processing units (GPUs) has become popular in scientific computations due to their low cost, impressive floating-point capabilities, high memory bandwidth, and low electrical power requirements. Discrete element method (DEM), closely related to molecular dynamics (MD), is a widely-accepted method for simulating particle systems, such as powders, bubbles, grains, etc. In this work, we present an efficient implementation of DEM on multiple GPUs. Linked-cell list method is used to accelerate recognition of contacting particles and is beneficial for communications among multiple GPUs. Moreover, asynchronous communication is used to reduce the time costs of communications. Physical validation is performed and parallel efficiency is evaluated, showing that our implementation has high computational efficiency and good scalability. The simulation of a large-scale dense particle system (up to 128,000,000) is achieved with a good parallel efficiency. The speed-up ratio in 16 GPUs vs 128 CPU cores is 10.39.

A physics-motivated Centroidal Voronoi Particle domain decomposition method

In this paper, we propose a novel domain decomposition method for large-scale simulations in continuum mechanics by merging the concepts of Centroidal Voronoi Tessellation (CVT) and Voronoi Particle dynamics (VP). The CVT is introduced to achieve a high-level compactness of the partitioning subdomains by the Lloyd algorithm which monotonically decreases the CVT energy. The number of computational elements between neighboring partitioning subdomains, which scales the communication effort for parallel simulations, is optimized implicitly as the generated partitioning subdomains are convex and simply connected with small aspect-ratios. Moreover, Voronoi Particle dynamics employing physical analogy with a tailored equation of state is developed, which relaxes the particle system towards the target partition with good load balance. Since the equilibrium is computed by an iterative approach, the partitioning subdomains exhibit locality and the incremental property. Numerical experiments reveal that the proposed Centroidal Voronoi Particle (CVP) based algorithm produces high-quality partitioning with high efficiency, independently of computational-element types. Thus it can be used for a wide range of applications in computational science and engineering.

Quantifying the exchange coupling in linear copper porphyrin oligomers.

Linear π-conjugated porphyrin oligomers are of significant current interest due to their potential applications as molecular wires. In this study we investigate electronic communication in linear butadiyne-linked copper porphyrin oligomers by electron paramagnetic resonance (EPR) spectroscopy via measurement of the exchange interaction, J, between the copper(ii) centers. The contributions of dipolar and exchange interactions to the frozen solution continuous wave (cw) EPR spectra of the compounds with two or more copper porphyrin units were explicitly accounted for in numerical simulations using a spin Hamiltonian approach. It is demonstrated that a complete numerical simulation of the powder spectrum of a large spin system with a Hamiltonian dimension of 26 244 and beyond can be made feasible by simulating the spectra in the time domain. The exchange coupling in the Cu2 dimer (CuCu distance 1.35 nm) is of the order of tens of MHz (Ĥ = -2JS1·S2) and is strongly modulated by low-energy molecular motions such as twisting of the molecule.

Charge-Transfer Modified Embedded-Atom Method for LiNi1-x-YCoxMnyO2 Cathode Materials

Benefiting from the high voltage and high capacity, transition metal (TM) layered oxides LiNi1-x-yCoxMnyO2 (NCM) constitute a very attractive family of compounds for the positive electrode of Li-ion batteries. Full quantum simulations such as density functional theory (DFT) method have provided many accurate atomic and electronic structure information, but these understanding is limited to a very small system size of 100-1000 atoms. For micro-scale structures simulations, e.g. , multiple phases and interfaces in-between, a large-scale simulation method that can simulate 1000-10000 atoms is required. However, high challenges are faced to develop such method for cathode materials, due to the complex TM-O bonding nature with both ionic and covalent components, the diverse interactions in a quinary system, and the continuous charge state variations during delithiation. In this work, we have developed a multi-scale simulation method for NCM cathode materials, based on Charge-Transfer Modified Embedded-Atom Method (CT-MEAM), in which both the long range Coulomb interaction due to charge transfer and the short range non-ionic interaction can be well described. Properties derived from this method match quite well with the corresponding DFT reference work and essential atomic structure behaviors have been modeled, including NCM interfaces, defect structures, size effects, and corresponding TM distribution interactions, etc. The presented simulation method can help facilitate large-scale atomic calculations required for the optimal design of NCM oxides applications.

Charge-transfer modified embedded-atom method for manganese oxides: Nanostructuring effects on MnO2 nanorods

In spite of the superior performance of nanostructured manganese oxides as materials for electro-catalysis, there has been low progress in their theoretical understanding because of the large scale of the required simulation cells. In this work, we have developed a large-scale simulation method for manganese oxides based on the Charge-Transfer Modified Embedded-Atom Method (CT-MEAM), which is not only capable of describing complex atomic bonding in transition metal oxides but also dynamic oxidation states, as those induced by the widely exposed surface area of nanostructures. The reliability of the method is confirmed by comparison with the corresponding ab initio reference results, and afterwards used to investigate the nanostructuring effects on manganese oxide nanorods. We found that nanostructuring introduces surface and edge structural distortions, driven by emerged unsaturated surface ions. These ions are highly mobile, and easily migrate to form new bonds with neighboring ions. Nanostructuring also weakens charge transfer between Mn and O ions located at the outer layers, thus bringing in more reductive Mn ions and enabling oxygen vacancy formation, which ultimately enhances the catalytic activity of manganese oxides. Our results also show that large scale surface defects, such as voids, could produce a similar effect than size reduction in improving the catalytic performance of the nanorod, thus avoiding undesirable amorphization instabilities, which are inevitable at reduced sizes. These findings will help to the optimization of the design and synthesis of nanostructured manganese oxides with enhanced catalytic activity.

First-principle optimal local pseudopotentials construction via optimized effective potential method.

The local pseudopotential (LPP) is an important component of orbital-free density functional theory, a promising large-scale simulation method that can maintain information on a material's electron state. The LPP is usually extracted from solid-state density functional theory calculations, thereby it is difficult to assess its transferability to cases involving very different chemical environments. Here, we reveal a fundamental relation between the first-principles norm-conserving pseudopotential (NCPP) and the LPP. On the basis of this relationship, we demonstrate that the LPP can be constructed optimally from the NCPP for a large number of elements using the optimized effective potential method. Specially, our method provides a unified scheme for constructing and assessing the LPP within the framework of first-principles pseudopotentials. Our practice reveals that the existence of a valid LPP with high transferability may strongly depend on the element.

A large-scale simulation method on complex ternary Li–Mn–O compounds for Li-ion battery cathode materials

To meet the requirement of large-scale simulation technics for Li-ion battery electrode materials, we introduce the charge-transfer modified embedded-atom method (CT-MEAM) in which the complex nature of the chemical bonding in transition metal (TM) oxides is described as a balance between metallic/covalent and ionic contributions by MEAM and a variable-charge model, respectively. The method is applied to Li2MnO3, and the parameterization is performed through fitting the energy–strain curves of Li2MnO3 under uniaxial, biaxial and hydrostatic strains to a training set from ab initio density-functional theory calculations. The CT-MEAM prediction of the critical physical properties such as charge states and redox potentials match quite well with the ab initio results in various Li–Mn–O compounds beyond Li2MnO3. The constructed Li–Mn–O phase diagram is also qualitatively consistent with the ab initio reference work. The excellent transferability ensures use of the present method for a wide range of oxidation states in complex ternary TM oxides. Therefore, it will facilitate large-scale atomistic calculations required for the optimal design of many TM oxide applications including lithium-ion battery cathode materials.

Sound Wave Propagation Anisotropy in Silver Nanoprisms: Characterization of Photoinduced Multiple Modes Using the Symmetric Molecular Dynamics Method

We present both theoretical and experimental studies of the sound velocity anisotropy and laser-induced multiple acoustic phonon oscillations in silver nanoprisms. A large-scale simulation method based on molecular dynamics and group theory for nanoprisms with C3v symmetry was developed. We observed the most prominent planar modes and the less studied vertical modes and found that the anisotropic vibrational properties of silver nanoprisms strongly depend on the particle’s aspect ratio (bisector length over thickness). For the case of a smaller aspect ratio, we identified these modes as the A1 modes with a larger in-phase radial atomic displacement. These modes can be classified into three main types: breathing modes and totally symmetric modes which are planar modes and thickness-related vertical modes. Due to the strong coupling with laterally confined modes, a broad band of thickness-related modes is found (in contrast to a single sharp peak for thin films), which leads to a fast decaying component in the oscillation spectrum. The calculated results and the sound velocity anisotropy along different crystal axes are in good agreement with our transient optical absorption experiments.

A distinct element method for large scale simulations of carbon nanotube assemblies

Mesoscale simulation techniques are becoming increasingly important due to the interest in complex mechanical problems involving nanomaterials. We propose applying the established macroscopic modeling concept of distinct spherical elements down to the mesoscale to simulate mechanical behavior of carbon nanotube systems. Starting from a microscopic description, the important interactions are encapsulated into two types of contact models that act simultaneously. Each individual nanotube is coarse-grained into a chain of spherical elements interacting by parallel-bonded contacts, representing the short-ranged covalent bonding. An anisotropic van der Waals model with aligning moments acts at the contact between elements located in different tubes to represent the long-ranged interactions. The promising potential of the proposed methodology to model large scale carbon nanotube assemblies is illustrated with several examples, including self-folding of individual nanotubes, mechanical testing of nanotube ropes, self-assembly of a high-porosity nanotube paper, and mechanical testing of a low-porosity nanotube paper.

Preconditioned iteration for saddle-point systems with bound constraints arising in contact problems

We introduce a new algorithm, which we shall call Mixed Iteration, or by the acronym Mix-It, for the solution of saddle point problems arising from minimization problems under equality and/or inequality constraints. Such constrained problems are found in many applications and we shall be specially interested in contact problems in solid mechanics. The algorithm is however general and can be applied to other situations. We present a real saddle point algorithm in the sense that the primal variables and Lagrange multipliers unknowns are solved simultaneously. The effectiveness of the proposed method for large-scale simulations is improved by a block preconditioner and an adequate treatment of inequality constraints. The block preconditioner involves the solution of two subsystems associated respectively with the primal and the dual variables. We present some numerical results to illustrate the performances of the proposed method.

Parallel Simulation of Cellular Flow in Microvessels Using a Particle Method

We developed a numerical method for large-scale simulations of cellular flow in microvessels. We employed a particle method, where all blood components were modeled using a finite number of particles. Red blood cell deformation was modeled by a spring network of membrane particles. A domain decomposition method was used for parallel implementation on distributed memory systems. In a strong scaling test up to 64 CPU cores, we obtained a linear speedup with the number of CPU cores, and demonstrated that our model can simulate O(10 3 ) red blood cells in vessels a few tens of micrometers in diameter. For quantitative validation, we analyzed the Fahraeus effect and the formation of a cell-depleted peripheral layer. Simulations were performed for tube hematocrit ranging from 20 to 45%, and microvessel diameters from 9 to 50 μm. Our numerical results were in good agreement with previous experimental results both for the discharge hematocrit and cell-depleted peripheral layer thickness.

On the modeling and simulation of non-hydrostatic dam break flows

The numerical simulation of three-dimensional dam break flows is discussed. A non-hydrostatic numerical model for free-surface flows is considered, which is based on the incompressible Navier---Stokes equations coupled with a volume-of-fluid approach. The numerical results obtained for a variety of benchmark problems show the validity of the numerical approach, in comparison with other numerical models, and allow to investigate numerically the non-hydrostatic three-dimensional effects, in particular for the usual test cases where hydrostatic approximations are known analytically. The numerical experiments on actual topographies, in particular the Malpasset dam break and the (hypothetical) break of the Grande-Dixence dam in Switzerland, also illustrate the capabilities of the method for large-scale simulations and real-life visualization.