Large scale simulation of pressure induced phase-field fracture propagation using Utopia.

Non-linear phase field models are increasingly used for the simulation of fracture propagation problems. The numerical simulation of fracture networks of realistic size requires the efficient parallel solution of large coupled non-linear systems. Although in principle efficient iterative multi-level methods for these types of problems are available, they are not widely used in practice due to the complexity of their parallel implementation. Here, we present Utopia, which is an open-source C++ library for parallel non-linear multilevel solution strategies. Utopia provides the advantages of high-level programming interfaces while at the same time a framework to access low-level data-structures without breaking code encapsulation. Complex numerical procedures can be expressed with few lines of code, and evaluated by different implementations, libraries, or computing hardware. In this paper, we investigate the parallel performance of our implementation of the recursive multilevel trust-region (RMTR) method based on the Utopia library. RMTR is a globally convergent multilevel solution strategy designed to solve non-convex constrained minimization problems. In particular, we solve pressure-induced phase-field fracture propagation in large and complex fracture networks. Solving such problems is deemed challenging even for a few fractures, however, here we are considering networks of realistic size with up to 1000 fractures.

Materials Cloud, a platform for open computational science.

Materials Cloud is a platform designed to enable open and seamless sharing of resources for computational science, driven by applications in materials modelling. It hosts (1) archival and dissemination services for raw and curated data, together with their provenance graph, (2) modelling services and virtual machines, (3) tools for data analytics, and pre-/post-processing, and (4) educational materials. Data is citable and archived persistently, providing a comprehensive embodiment of entire simulation pipelines (calculations performed, codes used, data generated) in the form of graphs that allow retracing and reproducing any computed result. When an AiiDA database is shared on Materials Cloud, peers can browse the interconnected record of simulations, download individual files or the full database, and start their research from the results of the original authors. The infrastructure is agnostic to the specific simulation codes used and can support diverse applications in computational science that transcend its initial materials domain.

CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations.

CP2K is an open source electronic structure and molecular dynamics software package to perform atomistic simulations of solid-state, liquid, molecular, and biological systems. It is especially aimed at massively parallel and linear-scaling electronic structure methods and state-of-the-art ab initio molecular dynamics simulations. Excellent performance for electronic structure calculations is achieved using novel algorithms implemented for modern high-performance computing systems. This review revisits the main capabilities of CP2K to perform efficient and accurate electronic structure simulations. The emphasis is put on density functional theory and multiple post-Hartree-Fock methods using the Gaussian and plane wave approach and its augmented all-electron extension.

Dynamics of the Bulk Hydrated Electron from Many-Body Wave-Function Theory.

The structure of the hydrated electron is a matter of debate as it evades direct experimental observation owing to the short life time and low concentrations of the species. Herein, the first molecular dynamics simulation of the bulk hydrated electron based on correlated wave-function theory provides conclusive evidence in favor of a persistent tetrahedral cavity made up by four water molecules, and against the existence of stable non-cavity structures. Such a cavity is formed within less than a picosecond after the addition of an excess electron to neat liquid water, with less regular cavities appearing as intermediates. The cavities are bound together by weak H-H bonds, the number of which correlates well with the number of coordinated water molecules, each type of cavity leaving a distinct spectroscopic signature. Simulations predict regions of negative spin density and a gyration radius that are both in agreement with experimental data.

Microcanonical RT-TDDFT simulations of realistically extended devices.

In this paper, real-time time-dependent density functional theory (RT-TDDFT) calculations of realistically sized nanodevices are presented. These microcanonical simulations rely on a closed boundary approach based on recent advances in the software package CP2K. The obtained results are compared to those derived from the open-boundary Non-equilibrium Green's Function (NEGF) formalism. A good agreement between the "current vs. voltage" characteristics produced by both methods is demonstrated for three representative device structures, a carbon nanotube field-effect transistor, a GeSe selector for crossbar arrays, and a conductive bridging random-access memory cell. Different approaches to extract the electrostatic contribution from the RT-TDDFT Hamiltonian and to incorporate the result into the NEGF calculations are presented.

Decisive Role of Perimeter Sites in Silica-Supported Ag Nanoparticles in Selective Hydrogenation of CO2 to Methyl Formate in the Presence of Methanol.

Methyl formate synthesis by hydrogenation of carbon dioxide in the presence of methanol offers a promising path to valorize carbon dioxide. In this work, silica-supported silver nanoparticles are shown to be a significantly more active catalyst for the continuous methyl formate synthesis than the known gold and copper counterparts, and the origin of the unique reactivity of Ag is clarified. Transient in situ and operando vibrational spectroscopy and DFT calculations shed light on the reactive intermediates and reaction mechanisms: a key feature is the rapid formation of surface chemical species in equilibrium with adsorbed carbon dioxide. Such species is assigned to carbonic acid interacting with water/hydroxyls on silica and promoting the esterification of formic acid with adsorbed methanol at the perimeter sites of Ag on SiO2 to yield methyl formate. This study highlights the importance of employing combined methodologies to verify the location and nature of active sites and to uncover fundamental catalytic reaction steps taking place at metal-support interfaces.

Machine Learning Adaptive Basis Sets for Efficient Large Scale Density Functional Theory Simulation.

It is chemically intuitive that an optimal atom centered basis set must adapt to its atomic environment, for example by polarizing toward nearby atoms. Adaptive basis sets of small size can be significantly more accurate than traditional atom centered basis sets of the same size. The small size and well conditioned nature of these basis sets leads to large saving in computational cost, in particular in a linear scaling framework. Here, it is shown that machine learning can be used to predict such adaptive basis sets using local geometrical information only. As a result, various properties of standard DFT calculations can be easily obtained at much lower costs, including nuclear gradients. In our approach, a rotationally invariant parametrization of the basis is obtained by employing a potential anchored on neighboring atoms to ultimately construct a rotation matrix that turns a traditional atom centered basis set into a suitable adaptive basis set. The method is demonstrated using MD simulations of liquid water, where it is shown that minimal basis sets yield structural properties in fair agreement with basis set converged results, while reducing the computational cost in the best case by a factor of 200 and the required flops by 4 orders of magnitude. Already a very small training set yields satisfactory results as the variational nature of the method provides robustness.

Efficient algorithms for large-scale quantum transport calculations.

Massively parallel algorithms are presented in this paper to reduce the computational burden associated with quantum transport simulations from first-principles. The power of modern hybrid computer architectures is harvested in order to determine the open boundary conditions that connect the simulation domain with its environment and to solve the resulting Schrödinger equation. While the former operation takes the form of an eigenvalue problem that is solved by a contour integration technique on the available central processing units (CPUs), the latter can be cast into a linear system of equations that is simultaneously processed by SplitSolve, a two-step algorithm, on general-purpose graphics processing units (GPUs). A significant decrease of the computational time by up to two orders of magnitude is obtained as compared to standard solution methods.

Nuclear Quantum Effects on Aqueous Electron Attachment and Redox Properties.

Nuclear quantum effects (NQEs) on the reduction and oxidation properties of small aqueous species (CO2, HO2, and O2) are quantified and rationalized by first-principles molecular dynamics and thermodynamic integration. Vertical electron attachment, or electron affinity, and detachment energies (VEA and VDE) are strongly affected by NQEs, decreasing in absolute value by 0.3 eV going from a classical to a quantum description of the nuclei. The effect is attributed to NQEs that lessen the solvent response upon oxidation/reduction. The reduction of solvent reorganization energy is expected to be general for small solutes in water. In the thermodynamic integral that yields the free energy of oxidation/reduction, these large changes enter with opposite sign, and only a small net effect (0.1 eV) remains. This is not obvious for CO2, where the integrand is strongly influenced by NQEs due to the onset of interaction of the reduced orbital with the conduction band of the liquid during thermodynamic integration. We conclude that NQEs might not have to be included in the computation of redox potentials, unless high accuracy is needed, but are important for VEA and VDE calculations.

Catalyst support effects on hydrogen spillover.

Hydrogen spillover is the surface migration of activated hydrogen atoms from a metal catalyst particle, on which they are generated, onto the catalyst support. The phenomenon has been much studied and its occurrence on reducible supports such as titanium oxide is established, yet questions remain about whether hydrogen spillover can take place on nonreducible supports such as aluminium oxide. Here we use the enhanced precision of top-down nanofabrication to prepare controlled and precisely tunable model systems that allow us to quantify the efficiency and spatial extent of hydrogen spillover on both reducible and nonreducible supports. We place multiple pairs of iron oxide and platinum nanoparticles on titanium oxide and aluminium oxide supports, varying the distance between the pairs from zero to 45 nanometres with a precision of one nanometre. We then observe the extent of the reduction of the iron oxide particles by hydrogen atoms generated on the platinum using single-particle in situ X-ray absorption spectromicroscopy applied simultaneously to all particle pairs. The data, in conjunction with density functional theory calculations, reveal fast hydrogen spillover on titanium oxide that reduces remote iron oxide nanoparticles via coupled proton-electron transfer. In contrast, spillover on aluminium oxide is mediated by three-coordinated aluminium centres that also interact with water and that give rise to hydrogen mobility competing with hydrogen desorption; this results in hydrogen spillover about ten orders of magnitude slower than on titanium oxide and restricted to very short distances from the platinum particle. We anticipate that these observations will improve our understanding of hydrogen storage and catalytic reactions involving hydrogen, and that our approach to creating and probing model catalyst systems will provide opportunities for studying the origin of synergistic effects in supported catalysts that combine multiple functionalities.

First-Principles Monte Carlo Simulations of Reaction Equilibria in Compressed Vapors.

Predictive modeling of reaction equilibria presents one of the grand challenges in the field of molecular simulation. Difficulties in the study of such systems arise from the need (i) to accurately model both strong, short-ranged interactions leading to the formation of chemical bonds and weak interactions arising from the environment, and (ii) to sample the range of time scales involving frequent molecular collisions, slow diffusion, and infrequent reactive events. Here we present a novel reactive first-principles Monte Carlo (RxFPMC) approach that allows for investigation of reaction equilibria without the need to prespecify a set of chemical reactions and their ideal-gas equilibrium constants. We apply RxFPMC to investigate a nitrogen/oxygen mixture at T = 3000 K and p = 30 GPa, i.e., conditions that are present in atmospheric lightning strikes and explosions. The RxFPMC simulations show that the solvation environment leads to a significantly enhanced NO concentration that reaches a maximum when oxygen is present in slight excess. In addition, the RxFPMC simulations indicate the formation of NO2 and N2O in mole fractions approaching 1%, whereas N3 and O3 are not observed. The equilibrium distributions obtained from the RxFPMC simulations agree well with those from a thermochemical computer code parametrized to experimental data.

Combining Linear-Scaling DFT with Subsystem DFT in Born-Oppenheimer and Ehrenfest Molecular Dynamics Simulations: From Molecules to a Virus in Solution.

In this work, methods for the efficient simulation of large systems embedded in a molecular environment are presented. These methods combine linear-scaling (LS) Kohn-Sham (KS) density functional theory (DFT) with subsystem (SS) DFT. LS DFT is efficient for large subsystems, while SS DFT is linear scaling with a smaller prefactor for large sets of small molecules. The combination of SS and LS, which is an embedding approach, can result in a 10-fold speedup over a pure LS simulation for large systems in aqueous solution. In addition to a ground-state Born-Oppenheimer SS+LS implementation, a time-dependent density functional theory-based Ehrenfest molecular dynamics (EMD) using density matrix propagation is presented that allows for performing nonadiabatic dynamics. Density matrix-based EMD in the SS framework is naturally linear scaling and appears suitable to study the electronic dynamics of molecules in solution. In the LS framework, linear scaling results as long as the density matrix remains sparse during time propagation. However, we generally find a less than exponential decay of the density matrix after a sufficiently long EMD run, preventing LS EMD simulations with arbitrary accuracy. The methods are tested on various systems, including spectroscopy on dyes, the electronic structure of TiO2 nanoparticles, electronic transport in carbon nanotubes, and the satellite tobacco mosaic virus in explicit solution.

Spin-Unrestricted Second-Order Møller-Plesset (MP2) Forces for the Condensed Phase: From Molecular Radicals to F-Centers in Solids.

To obtain consistent geometries for the computation of properties, nuclear gradients are essential. Here, we report a fully periodic Γ-point, massively parallel implementation of spin-unrestricted second-order Møller-Plesset (MP2) forces. It is based on the resolution-of-identity and Gaussian and plane waves approach to calculate electron repulsion integrals and is made available in the CP2K program. The algorithm is optimized for modern supercomputer architectures and is capable of employing both CPUs and GPUs. The asymptotic computational scaling is O(N(5)), which is observed for the systems containing more than 100 second-row atoms with triple-ζ quality basis sets. For smaller systems, the energy and forces can be evaluated within minutes; for bigger systems within hours, provided that thousands of processor units are available. Spin-unrestricted MP2 calculations are, computationally, ∼3 times more demanding than spin-restricted ones, but exhibit a better parallel performance. As a bonus, the gradient implementation allows for computing spin-density distributions not available in energy calculations. The method has been successfully used for computing the crystal structure of the TEMPO radical and to obtain the spin density distribution of an F-center in lithium fluoride at a consistent geometry. We expect this computationally efficient tool to be useful for the chemistry, materials science, and solid-state physics communities.

Calculation of Electrochemical Energy Levels in Water Using the Random Phase Approximation and a Double Hybrid Functional.

Understanding charge transfer at electrochemical interfaces requires consistent treatment of electronic energy levels in solids and in water at the same level of the electronic structure theory. Using density-functional-theory-based molecular dynamics and thermodynamic integration, the free energy levels of six redox couples in water are calculated at the level of the random phase approximation and a double hybrid density functional. The redox levels, together with the water band positions, are aligned against a computational standard hydrogen electrode, allowing for critical analysis of errors compared to the experiment. It is encouraging that both methods offer a good description of the electronic structures of the solutes and water, showing promise for a full treatment of electrochemical interfaces.

A generalized Poisson solver for first-principles device simulations.

Electronic structure calculations of atomistic systems based on density functional theory involve solving the Poisson equation. In this paper, we present a plane-wave based algorithm for solving the generalized Poisson equation subject to periodic or homogeneous Neumann conditions on the boundaries of the simulation cell and Dirichlet type conditions imposed at arbitrary subdomains. In this way, source, drain, and gate voltages can be imposed across atomistic models of electronic devices. Dirichlet conditions are enforced as constraints in a variational framework giving rise to a saddle point problem. The resulting system of equations is then solved using a stationary iterative method in which the generalized Poisson operator is preconditioned with the standard Laplace operator. The solver can make use of any sufficiently smooth function modelling the dielectric constant, including density dependent dielectric continuum models. For all the boundary conditions, consistent derivatives are available and molecular dynamics simulations can be performed. The convergence behaviour of the scheme is investigated and its capabilities are demonstrated.

Forces and stress in second order Møller-Plesset perturbation theory for condensed phase systems within the resolution-of-identity Gaussian and plane waves approach.

The forces acting on the atoms as well as the stress tensor are crucial ingredients for calculating the structural and dynamical properties of systems in the condensed phase. Here, these derivatives of the total energy are evaluated for the second-order Møller-Plesset perturbation energy (MP2) in the framework of the resolution of identity Gaussian and plane waves method, in a way that is fully consistent with how the total energy is computed. This consistency is non-trivial, given the different ways employed to compute Coulomb, exchange, and canonical four center integrals, and allows, for example, for energy conserving dynamics in various ensembles. Based on this formalism, a massively parallel algorithm has been developed for finite and extended system. The designed parallel algorithm displays, with respect to the system size, cubic, quartic, and quintic requirements, respectively, for the memory, communication, and computation. All these requirements are reduced with an increasing number of processes, and the measured performance shows excellent parallel scalability and efficiency up to thousands of nodes. Additionally, the computationally more demanding quintic scaling steps can be accelerated by employing graphics processing units (GPU's) showing, for large systems, a gain of almost a factor two compared to the standard central processing unit-only case. In this way, the evaluation of the derivatives of the RI-MP2 energy can be performed within a few minutes for systems containing hundreds of atoms and thousands of basis functions. With good time to solution, the implementation thus opens the possibility to perform molecular dynamics (MD) simulations in various ensembles (microcanonical ensemble and isobaric-isothermal ensemble) at the MP2 level of theory. Geometry optimization, full cell relaxation, and energy conserving MD simulations have been performed for a variety of molecular crystals including NH3, CO2, formic acid, and benzene.

Probing the structural and dynamical properties of liquid water with models including non-local electron correlation.

Water is a ubiquitous liquid that displays a wide range of anomalous properties and has a delicate structure that challenges experiment and simulation alike. The various intermolecular interactions that play an important role, such as repulsion, polarization, hydrogen bonding, and van der Waals interactions, are often difficult to reproduce faithfully in atomistic models. Here, electronic structure theories including all these interactions at equal footing, which requires the inclusion of non-local electron correlation, are used to describe structure and dynamics of bulk liquid water. Isobaric-isothermal (NpT) ensemble simulations based on the Random Phase Approximation (RPA) yield excellent density (0.994 g/ml) and fair radial distribution functions, while various other density functional approximations produce scattered results (0.8-1.2 g/ml). Molecular dynamics simulation in the microcanonical (NVE) ensemble based on Møller-Plesset perturbation theory (MP2) yields dynamical properties in the condensed phase, namely, the infrared spectrum and diffusion constant. At the MP2 and RPA levels of theory, ice is correctly predicted to float on water, resolving one of the anomalies as resulting from a delicate balance between van der Waals and hydrogen bonding interactions. For several properties, obtaining quantitative agreement with experiment requires correction for nuclear quantum effects (NQEs), highlighting their importance, for structure, dynamics, and electronic properties. A computed NQE shift of 0.6 eV for the band gap and absorption spectrum illustrates the latter. Giving access to both structure and dynamics of condensed phase systems, non-local electron correlation will increasingly be used to study systems where weak interactions are of paramount importance.

Efficient preconditioning of the electronic structure problem in large scale ab initio molecular dynamics simulations.

We present an improved preconditioning scheme for electronic structure calculations based on the orbital transformation method. First, a preconditioner is developed which includes information from the full Kohn-Sham matrix but avoids computationally demanding diagonalisation steps in its construction. This reduces the computational cost of its construction, eliminating a bottleneck in large scale simulations, while maintaining rapid convergence. In addition, a modified form of Hotelling's iterative inversion is introduced to replace the exact inversion of the preconditioner matrix. This method is highly effective during molecular dynamics (MD), as the solution obtained in earlier MD steps is a suitable initial guess. Filtering small elements during sparse matrix multiplication leads to linear scaling inversion, while retaining robustness, already for relatively small systems. For system sizes ranging from a few hundred to a few thousand atoms, which are typical for many practical applications, the improvements to the algorithm lead to a 2-5 fold speedup per MD step.

Redox potentials and acidity constants from density functional theory based molecular dynamics.

CONSPECTUS: All-atom methods treat solute and solvent at the same level of electronic structure theory and statistical mechanics. All-atom computation of acidity constants (pKa) and redox potentials is still a challenge. In this Account, we review such a method combining density functional theory based molecular dynamics (DFTMD) and free energy perturbation (FEP) methods. The key computational tool is a FEP based method for reversible insertion of a proton or electron in a periodic DFTMD model system. The free energy of insertion (work function) is computed by thermodynamic integration of vertical energy gaps obtained from total energy differences. The problem of the loss of a physical reference for ionization energies under periodic boundary conditions is solved by comparing with the proton work function computed for the same supercell. The scheme acts as a computational hydrogen electrode, and the DFTMD redox energies can be directly compared with experimental redox potentials. Consistent with the closed shell nature of acid dissociation, pKa estimates computed using the proton insertion/removal scheme are found to be significantly more accurate than the redox potential calculations. This enables us to separate the DFT error from other sources of uncertainty such as finite system size and sampling errors. Drawing an analogy with charged defects in solids, we trace the error in redox potentials back to underestimation of the energy gap of the extended states of the solvent. Accordingly the improvement in the redox potential as calculated by hybrid functionals is explained as a consequence of the opening up of the bandgap by the Hartree-Fock exchange component in hybrids. Test calculations for a number of small inorganic and organic molecules show that the hybrid functional implementation of our method can reproduce acidity constants with an uncertainty of 1-2 pKa units (0.1 eV). The error for redox potentials is in the order of 0.2 V.

The nature of excess electrons in anatase and rutile from hybrid DFT and RPA.

The behavior of excess electrons in undoped and defect free bulk anatase and rutile TiO2 has been investigated by state-of-the-art electronic structure methods including hybrid density functional theory (DFT) and the random phase approximation (RPA). Consistent with experiment, charge trapping and polaron formation is observed in both anatase and rutile. The difference in the anisotropic shape of the polarons is characterized, confirming for anatase the large polaron picture. For anatase, where polaron formation energies are small, charge trapping is observed also with standard hybrid functionals, provided the simulation cell is sufficiently large (864 atoms) to accommodate the lattice relaxation. Even though hybrid orbitals are required as a starting point for RPA in this system, the obtained polaron formation energies are relatively insensitive to the amount of Hartree-Fock exchange employed. The difference in trapping energy between rutile and anatase can be obtained accurately with both hybrid functionals and RPA. Computed activation energies for polaron hopping and delocalization clearly show that anatase and rutile might have different charge transport mechanisms. In rutile, only hopping is likely, whereas in anatase hopping and delocalization are competing. Delocalization will result in conduction-band-like and thus enhanced transport. Anisotropic conduction, in agreement with experimental data, is observed, and results from the tendency to delocalize in the [001] direction in rutile and the (001) plane in anatase. For future work, our calculations serve as a benchmark and suggest RPA on top on hybrid orbitals (PBE0 with 30% Hartree-Fock exchange), as a suitable method to study the rich chemistry and physics of TiO2.