Search

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., Competing Interests: Conflict of interestThe authors declare that they have no conflict of interest., (© The Author(s) 2021.)

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 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.

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., (© 2019 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim.)

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 SiO 2 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.

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.

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.

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 (NQEs) on the reduction and oxidation properties of small aqueous species (CO 2 , HO 2 , and O 2 ) 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 CO 2 , 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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Ice XV is the hydrogen-ordered form of the ice VI phase whose structure was predicted to be Cc and ferroelectric using periodic DFT approaches. However, neutron diffraction and Raman spectroscopy data show the structure to have P1̅ symmetry and to be antiferroelectric. Recent work1 using fragment-based MP2 and CCSD(T) approaches predicts the experimental structure as the ground state. We have reconsidered this problem using fully periodic MP2 and RPA approaches and find that the ferroelectric Cc structure is the lowest energy configuration. However, ubiquitously employed tinfoil boundary conditions stabilize polar structures. We suggest that ferroelectric Cc crystals can grow within a paraelectric ice VI matrix but may become unstable once a fraction of the matrix has become hydrogen-ordered. The reduction in dielectric constant causes P1̅ and other structures with small polarization to become favored, providing a possible resolution between observation and theoretical predictions.

The high overpotential in water oxidation on anodes is a limiting factor for the large-scale application of photoelectrochemical cells. To overcome this limitation, it is essential to understand the four proton-coupled electron transfer (PCET) steps in the reaction mechanism and their implications to the overpotential. Herein, a simple scheme to compute the energies of the PCET steps in water oxidation on the aqueous TiO2 surface using a hybrid density functional is described. An energy level diagram for fully decoupled electron- and proton-transfer reactions in which both electronic and protonic levels are placed on the same potential scale is also described. The level diagram helps to visualize the electronic and protonic components of the overpotential, and points out what are needed to improve. For TiO2, it is found that its catalytic activity is due to aligning the protonic energy levels in the PCET steps, while improving the activity requires also aligning the electronic levels., (© 2014 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.)

Monolayer sheets have gained attention due to the unique properties derived from their two-dimensional structure. One of the key challenges in sheet modification/synthesis is to exchange integral parts while keeping them intact. We describe site-to-site transmetalation of Zn(2+) in the netpoints of cm(2)-sized, metal-organic sheets by Fe(2+), Co(2+), and Pb(2+). This novel transformation was done both randomly and at predetermined patterns defined by photolithography to create monolayer sheets composed of different netpoints. All transmetalated sheets are mechanically strong enough to be spanned over 20 × 20 μm(2) sized holes. Density functional theory calculations provide both a model for the molecular structure of an Fe(2+)-based sheet and first insights into how transmetalation proceeds. Such transmetalated sheets with random and patterned netpoints can be considered as two-dimensional analogues of linear copolymers. Their nanoscale synthesis presents an advance in monolayer/polymer chemistry with applications in fields such as surface coating, molecular electronics, device fabrication, imaging, and sensing.

Covalent monolayer sheets in 2 hours: spreading of threefold anthracene-equipped shape-persistent and amphiphilic monomers at the air/water interface followed by a short photochemical treatment provides access to infinitely sized, strictly monolayered, covalent sheets with in-plane elastic modulus in the range of 19 N/m., (© 2013 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.)

The dielectric properties of the hydrogen disordered hexagonal phase (Ih) of water ice have been computed using density functional theory (DFT) based Monte Carlo simulations in the isobaric-isothermal ensemble. Temperature dependent data yield a fit for the Curie-Weiss law of the system and hence a prediction of the temperature of the phase transition from the Ih phase to the hydrogen ordered ice XI phase. Direct simulations around the phase transition temperature confirm and refine the predicted phase transition temperatures and provide data for further properties, such as the linear thermal expansion coefficient. Results have been obtained with both hybrid and semilocal density functionals, which yields insight in the performance of the electronic structure method. In particular, the hybrid functional yields significantly more realistic dielectric constants than the semilocal variant, namely ε ≈ 116 as opposed to ε ≈ 151 at 273 K (εexperiment = 95). This can be attributed to the tendency of semilocal functionals to be biased to configurations with a large dipole moment, and their overestimation of the dipole moments of these configurations. This is also reflected in the estimates of the Ih/XI transition temperature, which is 70-80 and 90-100 K for the hybrid and semilocal functional respectively. DFT based sampling of the millions of configurations necessary for this work has been enabled by a Tree Monte Carlo algorithm, designed for massively parallel computers.

Despite recent progress in linear scaling (LS) density function theory (DFT), the computational cost of the existing LS methods remains too high for a widespread adoption at present. In this work, we exploit nonorthogonal localized molecular orbitals to develop a series of LS methods for molecular systems with a low computational overhead. High efficiency of the proposed methods is achieved with a new robust two-stage variational procedure or by replacing the optimization altogether with an accurate nonself-consistent approach. We demonstrate that, even for challenging condensed-phase systems, the implemented LS methods are capable of extending the range of accurate DFT simulations to molecular systems that are an order of magnitude larger than those previously treated.

The second-order Møller-Plesset perturbation energy (MP2) and the Random Phase Approximation (RPA) correlation energy are increasingly popular post-Kohn-Sham correlation methods. Here, a novel algorithm based on a hybrid Gaussian and Plane Waves (GPW) approach with the resolution-of-identity (RI) approximation is developed for MP2, scaled opposite-spin MP2 (SOS-MP2), and direct-RPA (dRPA) correlation energies of finite and extended system. The key feature of the method is that the three center electron repulsion integrals (μν|P) necessary for the RI approximation are computed by direct integration between the products of Gaussian basis functions μν and the electrostatic potential arising from the RI fitting densities P. The electrostatic potential is obtained in a plane waves basis set after solving the Poisson equation in Fourier space. This scheme is highly efficient for condensed phase systems and offers a particularly easy way for parallel implementation. The RI approximation allows to speed up the MP2 energy calculations by a factor 10 to 15 compared to the canonical implementation but still requires O(N(5)) operations. On the other hand, the combination of RI with a Laplace approach in SOS-MP2 and an imaginary frequency integration in dRPA reduces the computational effort to O(N(4)) in both cases. In addition to that, our implementations have low memory requirements and display excellent parallel scalability up to tens of thousands of processes. Furthermore, exploiting graphics processing units (GPU), a further speedup by a factor ∼2 is observed compared to the standard only CPU implementations. In this way, RI-MP2, RI-SOS-MP2, and RI-dRPA calculations for condensed phase systems containing hundreds of atoms and thousands of basis functions can be performed within minutes employing a few hundred hybrid nodes. In order to validate the presented methods, various molecular crystals have been employed as benchmark systems to assess the performance, while solid LiH has been used to study the convergence with respect to the basis set and system size in the case of RI-MP2 and RI-dRPA.

The electronic states of aqueous species can mix with the extended states of the solvent if they are close in energy to the band edges of water. Using density functional theory-based molecular dynamics simulation, we show that this is the case for OH(-) and Cl(-). The effect is, however, badly exaggerated by the generalized gradient approximation leading to systematic underestimation of redox potentials and spurious nonlinearity in the solvent reorganization. Drawing a parallel to charged defects in wide gap solid oxides, we conclude that misalignment of the valence band of water is the main source of error turning the redox levels of OH(-) and Cl(-) in resonant impurity states. On the other hand, the accuracy of energies of levels corresponding to strongly negative redox potentials is acceptable. We therefore predict that mixing of the vertical attachment level of CO2 and the unoccupied states of water is a real effect.

A novel algorithm, based on a hybrid Gaussian and plane waves (GPW) approach, is developed for the canonical second-order Møller-Plesset perturbation energy (MP2) of finite and extended systems. The key aspect of the method is that the electron repulsion integrals (ia|λσ) are computed by direct integration between the products of Gaussian basis functions λσ and the electrostatic potential arising from a given occupied-virtual pair density ia. The electrostatic potential is obtained in a plane waves basis set after solving the Poisson equation in Fourier space. In particular, for condensed phase systems, this scheme is highly efficient. Furthermore, our implementation has low memory requirements and displays excellent parallel scalability up to 100 000 processes. In this way, canonical MP2 calculations for condensed phase systems containing hundreds of atoms or more than 5000 basis functions can be performed within minutes, while systems up to 1000 atoms and 10 000 basis functions remain feasible. Solid LiH has been employed as a benchmark to study basis set and system size convergence. Lattice constants and cohesive energies of various molecular crystals have been studied with MP2 and double-hybrid functionals.

In this work, the applicability and performance of a linear scaling algorithm is investigated for three-dimensional condensed phase systems. A simple but robust approach based on the matrix sign function is employed together with a thresholding matrix multiplication that does not require a prescribed sparsity pattern. Semiempirical methods and density functional theory have been tested. We demonstrate that self-consistent calculations with 1 million atoms are feasible for simple systems. With this approach, the computational cost of the calculation depends strongly on basis set quality. In the current implementation, high quality calculations for dense systems are limited to a few hundred thousand atoms. We report on the sparsities of the involved matrices as obtained at convergence and for intermediate iterations. We investigate how determining the chemical potential impacts the computational cost for very large systems.

We have investigated the structure and dynamics of thin water films adsorbed on TiO2(110) using infrared reflection-absorption spectroscopy (IRAS) and ab initio molecular dynamics. Infrared spectra were obtained for s- and p-polarized light with the plane of incidence parallel to the [001] and [11̅0] azimuths for water coverages ≤ 4 monolayers. The spectra indicate strong anisotropy in the water films. The vibrational densities of states predicted by the ab initio simulations for 1 and 2 monolayer coverages agree well with the observations. The results provide new insight into the structure of water on TiO2(110) and resolve a long-standing puzzle regarding the hydrogen bonding between molecules in the first and second monolayers on this surface. The results also demonstrate the capabilities of polarization- and azimuth-resolved IRAS for investigating the structure and dynamics of adsorbates on dielectric substrates.

We have carried out "first-principles" Born-Oppenheimer molecular dynamics (BOMD) simulations of the phosphate ions H₂PO₄⁻ and HPO₄²⁻ in liquid water and have calculated their IR spectra by Fourier transform techniques from the trajectories. IR bands were assigned by a so-called "generalized normal coordinate analysis". The effects of including Hartree-Fock (HF) exchange into the density functional theory (DFT) computation of forces were studied by comparing results obtained with the well-known BP, BLYP, and B3LYP functionals. The neglect of dispersion in the functionals was empirically corrected. The inclusion of HF exchange turned out to yield dramatically improved and, thus, quite accurate descriptions of the IR spectra observed for H₂PO₄⁻ and HPO₄²⁻ in aqueous solution. An analysis of earlier computational results (Klähn, M. et al. J. Phys. Chem. A 2004, 108, 6186-6194) on these vibrational spectra, which had been obtained in a hybrid setting combining a BP description of the respective phosphate with a simple molecular mechanics (MM) model of its aqueous environment, revealed three different sources of error, (i) the BP force field of the phosphates is much too soft and would have required a substantial scaling of frequencies, (ii) the oversimplified water force field entailed incorrect solvation structures and, thus, qualitatively wrong patterns of solvatochromic band shifts, and (iii) quantitative frequency computations additionally required the inclusion of HF exchange. Thus, the results of the B3LYP BOMD simulations do not only characterize physical properties like the IR spectra or the solvation structures of the phosphate systems but also provide clues for the future design of simplified but nevertheless reasonably accurate DFT/MM methods applicable to phosphates.

Understanding the properties of hydrated electrons, which were first observed using pulse radiolysis of water in 1962, is crucial because they are key species in many radiation chemistry processes. Although time-resolved spectroscopic studies and molecular simulations have shown that an electron in water (prepared, for example, by water photoionization) relaxes quickly to a localized, cavity-like structure ∼2.5 Å in radius, this picture has recently been questioned. In another experimental approach, negatively charged water clusters of increasing size were studied with photoelectron and IR spectroscopies. Although small water clusters can bind an excess electron, their character is very different from bulk hydrated species. As data on electron binding in liquid water have become directly accessible experimentally, the cluster-to-bulk extrapolations have become a topic of lively debate. Quantum electronic structure calculations addressing experimental measurables have, until recently, been largely limited to small clusters; extended systems were approached mainly with pseudopotential calculations combining a classical description of water with a quantum mechanical treatment of the excess electron. In this Account, we discuss our investigations of electrons solvated in water by means of ab initio molecular dynamics simulations. This approach, applied to a model system of a negatively charged cluster of 32 water molecules, allows us to characterize structural, dynamical, and reactive aspects of the hydrated electron using all of the system's valence electrons. We show that under ambient conditions, the electron localizes into a cavity close to the surface of the liquid cluster. This cavity is, however, more flexible and accessible to water molecules than an analogous area around negatively charged ions. The dynamical process of electron attachment to a neutral water cluster is strongly temperature dependent. Under ambient conditions, the electron relaxes in the liquid cluster and becomes indistinguishable from an equilibrated, solvated electron on a picosecond time scale. In contrast, for solid, cryogenic systems, the electron only partially localizes outside of the cluster, being trapped in a metastable, weakly bound "cushion-like" state. Strongly bound states under cryogenic conditions could only be prepared by cooling equilibrated, liquid, negatively charged clusters. These calculations allow us to rationalize how different isomers of electrons in cryogenic clusters can be observed experimentally. Our results also bring into question the direct extrapolation of properties of cryogenic, negatively charged water clusters to those of electrons in the bulk liquid. Ab initio molecular dynamics represents a unique computational tool for investigating the reactivity of the solvated electron in water. As a prototype, the electron-proton reaction was followed in the 32-water cluster. In accord with experiment, the molecular mechanism is a proton transfer process that is not diffusion limited, but rather controlled by a proton-induced deformation of the excess electron's solvent shell. We demonstrate the necessary ingredients of a successful density functional methodology for the hydrated electron that avoids potential pitfalls, such as self-interaction error, insufficient basis set, or lack of dispersion interactions. We also benchmark the density functional theory methods and outline the path to faithful ab initio simulations of dynamics and reactivity of electrons solvated in extended aqueous systems.

Proton transfer reactions, including acid-base recombination, are commonly considered to occur 'nearly instantaneously'. However, their actual time scales may stretch far into the microsecond range, as acid-base reactions are diffusion controlled and the concentrations are low near neutral pH. The interplay of competing bases in the pH relaxation is illustrated using a model acid-base system consisting of o-nitrobenzaldehyde (oNBA) as a proton cage and acetate ions and hydroxyl ions as bases. The kinetically controlled behavior leads to highly counterintuitive states, i.e. acetate ions are transiently protonated for hundreds of nanoseconds despite the presence of a much stronger base OH-.

The ultrafast dynamics of the cationic hole formed in bulk liquid water following ionization is investigated by ab initio molecular dynamics simulations and an experimentally accessible signature is suggested that might be tracked by femtosecond pump-probe spectroscopy. This is one of the fastest fundamental processes occurring in radiation-induced chemistry in aqueous systems and biological tissue. However, unlike the excess electron formed in the same process, the nature and time evolution of the cationic hole has been hitherto little studied. Simulations show that an initially partially delocalized cationic hole localizes within ~30 fs after which proton transfer to a neighboring water molecule proceeds practically immediately, leading to the formation of the OH radical and the hydronium cation in a reaction which can be formally written as H(2)O(+) + H(2)O → OH + H(3)O(+). The exact amount of initial spin delocalization is, however, somewhat method dependent, being realistically described by approximate density functional theory methods corrected for the self-interaction error. Localization, and then the evolving separation of spin and charge, changes the electronic structure of the radical center. This is manifested in the spectrum of electronic excitations which is calculated for the ensemble of ab initio molecular dynamics trajectories using a quantum mechanics/molecular mechanics (QM∕MM) formalism applying the equation of motion coupled-clusters method to the radical core. A clear spectroscopic signature is predicted by the theoretical model: as the hole transforms into a hydroxyl radical, a transient electronic absorption in the visible shifts to the blue, growing toward the near ultraviolet. Experimental evidence for this primary radiation-induced process is sought using femtosecond photoionization of liquid water excited with two photons at 11 eV. Transient absorption measurements carried out with ~40 fs time resolution and broadband spectral probing across the near-UV and visible are presented and direct comparisons with the theoretical simulations are made. Within the sensitivity and time resolution of the current measurement, a matching spectral signature is not detected. This result is used to place an upper limit on the absorption strength and/or lifetime of the localized H(2)O(+) ((aq)) species., (© 2011 American Institute of Physics)

We compare three different methods for direct energy minimization in electronic structure calculations where the gradient of the energy functional with respect to the molecular orbitals is available. These methods make use of the preconditioned gradient to increase robustness. An orbital transformation is used to ensure that the orthogonality constraint on the orbitals remains satisfied when using standard minimization methods. In addition, we propose an adaptive scheme for estimating the curvature of the energy functional to increase the performance of a line search free quasi-Newton method. We show that the performance of all methods is similar when robustness of the methods is ensured.

With a combination of transient pump-probe IR spectroscopy and ab initio molecular dynamics, the controversial pico- and nanosecond steps of the o-nitrobenzaldehyde (oNBA) photoreaction have been investigated in aqueous solution. In this way, the measured reaction kinetics have been complemented with an atomistic picture of the reactive events as obtained with unbiased simulations in explicit solvent. Our results allow for a detailed description of the oNBA proton photorelease, a process of fundamental importance and relevant to the use of oNBA as a proton cage in many experiments. In a first step, a stable ketene intermediate is formed on a subpicosecond time scale. This intermediate reacts in a solvent assisted way with an OH transfer to produce nitrosobenzoic acid with a characteristic time of 7 ps. Finally, in permitting pH conditions, this product molecule dissociates a carboxyl proton with a 21 ns time constant. The particular combination of theory and experiment employed in this work appears to be sufficiently general and powerful to find widespread application in the study of ultrafast reactive systems.

The calculation of Hartree-Fock exchange (HFX) is computationally demanding for large systems described with high-quality basis sets. In this work, we show that excellent performance and good accuracy can nevertheless be obtained if an auxiliary density matrix is employed for the HFX calculation. Several schemes to derive an auxiliary density matrix from a high-quality density matrix are discussed. Key to the accuracy of the auxiliary density matrix methods (ADMM) is the use of a correction based on standard generalized gradient approximations for HFX. ADMM integrates seamlessly in existing HFX codes and, in particular, can be employed in linear scaling implementations. Demonstrating the performance of the method, the effect of HFX on the structure of liquid water is investigated in detail using Born-Oppenheimer molecular dynamics simulations (300 ps) of a system of 64 molecules. Representative for large systems are calculations on a solvated protein (Rubredoxin), for which ADMM outperforms the corresponding standard HFX implementation by approximately a factor 20.

Using density functional theory we investigate whether intrinsic defects in ice surface segregate. We predict that hydronium, hydroxide, and the Bjerrum L- and D-defects are all more stable at the surface. However, the energetic cost to create a D-defect at the surface and migrate it into the bulk crystal is smaller than its bulk formation energy. Absolute and relative segregation energies are sensitive to the surface structure of ice, especially the spatial distribution of protons associated with dangling hydrogen bonds. It is found that the basal plane surface of hexagonal ice increases the bulk concentration of Bjerrum defects, strongly favoring D-defects over L-defects. Dangling protons associated with undercoordinated water molecules are preferentially injected into the crystal bulk as Bjerrum D-defects, leading to a surface dipole that attracts hydronium ions. Aside from the disparity in segregation energies for the Bjerrum defects, we find the interactions between defect species to be very finely balanced; surface segregation energies for hydronium and hydroxide species and trapping energies of these ionic species with Bjerrum defects are equal within the accuracy of our calculations. The mobility of the ionic hydronium and hydroxide species is greatly reduced at the surface in comparison to the bulk due to surface sites with high trapping affinities. We suggest that, in pure ice samples, the surface of ice will have an acidic character due to the presence of hydronium ions. This may be important in understanding the reactivity of ice particulates in the upper atmosphere and at the boundary layer.

A highly efficient mechanism for the regeneration of the cis-bis(isothiocyanato)bis(2,2'-bipyridyl-4,4'-dicarboxylato)-ruthenium(II) sensitizing dye (N3) by I(-) in acetonitrile has been identified by using molecular dynamics simulation based on density functional theory. Barrier-free complex formation of the oxidized dye with both I(-) and , and facile dissociation of and from the reduced dye are key steps in this process. In situ vibrational spectroscopy confirms the reversible binding of I(2) to the thiocyanate group. Additionally, simulations of the electrolyte near the interface suggest that acetonitrile is able to cover the (101) surface of anatase with a passivating layer that inhibits direct contact of the redox mediator with the oxide, and that the solvent structure specifically enhances the concentration of I(-) at a distance which further favors rapid dye regeneration.

Solvated electrons are ubiquitous intermediates in radiation-induced processes, with their lifetime being determined by quenching processes, such as the direct reaction with protons under acidic conditions. Ab initio molecular dynamics simulations allow us to unravel with molecular resolution the ultrafast reaction mechanism by which the electron and proton react in water. The path to a successful reaction involves a distortion and contraction of the hydrated electron and a rapid proton motion along a chain of hydrogen bonds, terminating on the water molecule most protruding into the electron cloud. This fundamental reaction is thus decidedly shown to be of a proton-transfer rather than electron-transfer character. Due to the desolvation penalty connected with breaking of the hydration shells of these charged particles, the reaction is, however, not diffusion-limited, in agreement with the interpretation of kinetics measurements.

Hartree-Fock exchange with a truncated Coulomb operator has recently been discussed in the context of periodic plane-waves calculations [Spencer, J.; Alavi, A. Phys. Rev. B: Solid State, 2008, 77, 193110]. In this work, this approach is extended to Gaussian basis sets, leading to a stable and accurate procedure for evaluating Hartree-Fock exchange at the Γ-point. Furthermore, it has been found that standard hybrid functionals can be transformed into short-range functionals without loss of accuracy. The well-defined short-range nature of the truncated exchange operator can naturally be exploited in integral screening procedures and makes this approach interesting for both condensed phase and gas phase systems. The presented Hartree-Fock implementation is massively parallel and scales up to ten thousands of cores. This makes it feasible to perform highly accurate calculations on systems containing thousands of atoms or ten thousands of basis functions. The applicability of this scheme is demonstrated by calculating the cohesive energy of a LiH crystal close to the Hartree-Fock basis set limit and by performing an electronic structure calculation of a complete protein (rubredoxin) in solution with a large and flexible basis set.

We present herein a comprehensive density functional theory study toward assessing the accuracy of two popular gradient-corrected exchange correlation functionals on the structure and density of liquid water at near ambient conditions in the isobaric-isothermal ensemble. Our results indicate that both PBE and BLYP functionals under predict the density and over structure the liquid. Adding the dispersion correction due to Grimme (1, 2) improves the predicted densities for both BLYP and PBE in a significant manner. Moreover, the addition of the dispersion correction for BLYP yields an oxygen-oxygen radial distribution function in excellent agreement with experiment. Thus, we conclude that one can obtain a very satisfactory model for water using BLYP and a correction for dispersion.

The silicate oligomerization reaction is key to sol-gel chemistry and zeolite synthesis. Numerous experimental and theoretical studies have addressed the physical chemistry of silicate oligomers in the prenucleation stage of siliceous zeolite formation. Here we report a study of a silica condensation reaction in aqueous solution in the presence of counter ions (Li(+) and NH(4)(+)). Ab-initio molecular dynamics simulations have been used to construct reaction energy diagrams including transition state free energies. Contact with Li+ as well as NH(4)(+) increases the activation energies of the dimerization step compared to the situation in the absence of counterions. The presence of NH(4)(+) has no effect on consecutive oligomerization steps. Hence NH(4)(+) will increase the relative formation rate of larger oligomers.

Combining photoemission and electrochemical data from the literature we argue that the difference between the vertical and adiabatic ionization energy of the aqueous hydroxide anion is 2.9 eV. We then use density functional theory based molecular dynamics to show that the solvent response to ionization is nonlinear. Adding this to the experimental data we predict a 4.1 eV difference between the energy for vertical attachment of an electron to the aqueous hydroxyl radical and the corresponding adiabatic electron affinity. This places the state accepting the electron only 2.2 eV below vacuum or 7.7 eV above the edge of the valence band of water.

Sulfuric acid and water clusters are important for new particle formation in the atmosphere. Recent experimental studies demonstrate that critical clusters in diverse atmospheric environments contain two acid molecules and may also include additional N-containing molecules (i.e., a base). We use first-principles molecular dynamics simulations to show that the presence of two sulfuric acid molecules in (H2SO4)m x base x (H2O)6 clusters is always sufficient to form a double ion, whereas a single acid molecule, even in the presence of a base, is not.

We present a computational study of the structure and dynamics of an excess electron in a medium-sized water cluster aimed at addressing the question of interior vs exterior solvation. Ab initio Born-Oppenheimer molecular dynamics simulations were performed within the DFT framework, employing a hybrid Gaussian and plane-wave formalism together with the PBE exchange-correlation functional and norm-conserving pseudopotentials. Analysis of a 15-ps trajectory allowed us to reach the following conclusions: (i) the excess electron is predominantly located at the cluster surface (even if it is initially placed in the interior), (ii) the computed electron binding energies correlate with the electron localization rather than with its bulk vs surface location, and (iii) a dynamical interconversion between two different H-bond patterns around the electron occurs. The computed electron binding energies and the most relevant features of the IR spectrum are in a very good agreement with results of previous experimental studies.

The spectroscopic signatures of proton transfer in the water dimer cation were investigated. The six lowest electronic states were characterized along the reaction coordinate using the equation-of-motion coupled-cluster with single and double substitutions method for ionized systems. The nature of the dimer states was explained in terms of the monomer states using a qualitative molecular orbital framework. We found that proton transfer induces significant changes in the electronic spectrum, thus suggesting that time-resolved electronic femtosecond spectroscopy is an effective strategy to monitor the dynamics following ionization. The electronic spectra at vertical and proton-transferred configurations include both local excitations (features similar to those of the monomers) and charge-transfer bands. Ab initio calculations were used to test the performance of a self-interaction correction for density functional theory (DFT). The corrected DFT/BLYP method is capable of quantitatively reproducing the proper energetic ordering of the (H2O)2(+) isomers and thus is a reasonable approach for calculations of larger systems.

Ab initio molecular dynamics simulations with hybrid density functionals have so far found little application due to their computational cost. In this work, an implementation of the Hartree-Fock exchange is presented that is specifically targeted at ab initio molecular dynamics simulations of medium sized systems. We demonstrate that our implementation, which is available as part of the CP2K/Quickstep program, is robust and efficient. Several prescreening techniques lead to a linear scaling cost for integral evaluation and storage. Integral compression techniques allow for in-core calculations on systems containing several thousand basis functions. The massively parallel implementation respects integral symmetry and scales up to hundreds of CPUs using a dynamic load balancing scheme. A time-reversible multiple time step scheme, exploiting the difference in computational efficiency between hybrid and local functionals, brings further time savings. With extensive simulations of liquid water, we demonstrate the ability to perform, for several tens of picoseconds, ab initio molecular dynamics based on hybrid functionals of systems in the condensed phase containing a few thousand Gaussian basis functions.