Search

Rate theory and DFT calculations of hydrogen evolution reaction (HER) on MoS2 with Co, Ni and Pt impurities show the significance of dihydrogen (H2*) complex where both hydrogen atoms are interacting with the surface. Stabilization of such a complex affects the competing Volmer-Heyrovsky (direct H2 release) and Volmer-Tafel (H2* intermediate) pathways. The resulting evolution proceeds with a very small overpotential for all dopants ($\eta$ = 0.1 to 0.2 V) at 25% edge substitution, significantly reduced from the already low $\eta$ = 0.27 V for the undoped edge. At full edge substitution, Co-MoS2 remains highly active ($\eta$ = 0.18 V) while Ni- and Pt-MoS2 are deactivated ($\eta$ = 0.4 to 0.5 V) due to unfavorable interaction with H2*. Instead of the single S-vacancy, the site of intrinsic activity in the basal plane was found to be the undercoordinated central Mo-atom in threefold S-vacancy configurations, enabling hydrogen evolution with $\eta$ = 0.52 V via a H2* intermediate. The impurity atoms interact favorably with the intrinsic sulfur vacancies on the basal plane, stabilizing but simultaneously deactivating the triple vacancy configuration. The calculated shifts in overpotential are consistent with reported measurements, and the dependence on doping level may explain variations in experimental observations.

Electrochemical reduction is a promising way to make use of CO$_2$ as feedstock for generating renewable fuel and valuable chemicals. Several metals can be used in the electrocatalyst to generate CO and formic acid but hydrogen formation is an unwanted side reaction that can even be dominant. The lack of selectivity is in general a significant problem, but silver-based electrocatalysts have been shown to be highly selective for CO with over over 90% faradaic efficiency when the applied voltage is below -1 V vs. RHE. Hydrogen formation is then insignificant and little formate is formed even though it is thermodynamically favored. We present calculations of the activation free energy for the various elementary steps as a function of applied voltage at the three low index facets, Ag(111), Ag(100) and Ag(110), as well as experimental measurements on polycrystalline electrodes, to identify the reason for this high selectivity. The formation of formic acid is suppressed because of the low coverage of adsorbed hydrogen and kinetic hindrance to the formation of the HCOO* intermediate, while *COOH, a key intermediate in CO formation, is thermodynamically unstable until the applied voltage reaches -1 V vs. RHE, at which point the kinetics for its formation are more favorable than for hydrogen. The calculated results are consistent with experimental measurements carried out for acidic conditions and provide an atomic scale insight into the high CO selectivity of silver-based electrocatalysts.

We apply a recently proposed computational protocol for a neural-network-supported configuration interaction (NN CI) calculation to the paradigmatic N$_2$ molecule. By comparison of correlation energy, binding energy, and the full dissociation curve to experimental and full CI benchmarks, we demonstrate the applicability and robustness of our approach for the first time in the context of molecular systems, and offer thereby a new complementary tool in the family of machine-learning-based computation methods. The main advantage of the method lies in the efficiency of the neural-network-selected many-body basis set. Specifically, we approximate full CI results obtained on bases of $\approx 10^{10}$ Slater Determinants with only $\approx10^{5}$ determinants with good accuracy. The high efficiency of the NN CI approach underlines its potential for broader applications such as structural optimizations and even computation of spectroscopic observables in systems for which computational resources are a limiting factor., Comment: 8 pages, 4 figures

The relaxation time of several second generation molecular motors is analysed by calculating the minimum energy path between the metastable and stable states and evaluating the transition rate within harmonic transition state theory based on energetics obtained from density functional theory. Comparison with published experimental data shows remarkably good agreement and demonstrates the predictive capability of the theoretical approach. While previous measurements by Feringa and coworkers [Chem.\,Eur.\,J.\,(2017) 23, 6643] have shown that a replacement of the stereogenic hydrogen by a fluorine atom increases the relaxation time because of destabilization of the transition state for the thermal helix inversion, we find that a replacement of CH$_3$ by a CF$_3$ group at the same site shortens the relaxation time because of elevated energy of the metastable state without a significant shift in the transition state energy. Since these two fluorine substitutions have an opposite effect on the relaxation time, the two combined can provide a way to fine tune the rotational speed of a molecular motor.

Copper is a good CO2 electroreduction catalyst as products beyond CO form, but efficiency and selectivity is low. Experiments have shown that admixture of other elements can help, and computational screening studies have pointed out various promising candidates based on the adsorption of a single CO molecule as a descriptor. Our calculations of CO adsorption on surfaces where a first row transition metal atom replaces a Cu atom show that multiple CO molecules, not just one, bind to the substitutional atom. For Fe, Co, and Ni atoms, a decrease in binding energy is found, but the reverse trend, namely increasing bond strength, is found for V, Cr, and Mn and the first three CO molecules. Magnetic moment, charge, and position of the substitutional atom are also strongly affected by the CO adsorption in most cases. Magnetic moment is stepwise reduced to zero, and the outward displacement of the substitutional atom increased., Comment: 8 pages, 2 figures

A method for locating first order saddle points on the energy surface of magnetic systems is described and several applications presented. The starting point of the iterative algorithm involved in the method can be anywhere, even close to a local energy minimum representing an initial state of a magnetic system and, in contrast to chain-of-states methods, the final state need not be specified. Convergence on the saddle points is guided by a negative energy gradient whose component along the minimum mode of the system is inverted, effectively transforming the neighbourhood of the saddle point to that of a local minimum. The method requires only the lowest two eigenvalues and corresponding eigenvectors of the Hessian of the system's energy and they are found using a quasi-Newton limited-memory Broyden-Fletcher-Goldfarb-Shanno solver for the minimization of the Rayleigh quotient without evaluation of the Hessian itself. The efficient implementation of the method and its linear scaling with the system size make it applicable to large systems. Applications are presented to transitions in systems that reveal significant complexity of co-existing magnetic states, such as skyrmions, skyrmion bags, skyrmion tubes, chiral bobbers and globules., Comment: 16 Pages, 10 Figures

Excited electronic states of molecules and solids play a fundamental role in fields such as catalysis and electronics. In electronic structure calculations, excited states typically correspond to saddle points on the surface described by the variation of the energy as a function of the electronic degrees of freedom. A direct optimization algorithm based on generalized mode following is presented for density functional calculations of excited states. While conventional direct optimization methods based on quasi-Newton algorithms usually converge to the stationary point closest to the initial guess, even minima, the generalized mode following approach systematically targets a saddle point of a specific order l by following the l lowest eigenvectors of the electronic Hessian up in energy. This approach thereby recasts the challenging saddle point search as a minimization, enabling the use of efficient and robust minimization algorithms. The initial guess orbitals and the saddle point order of the target excited state solution are evaluated by performing an initial step of constrained optimization freezing the electronic degrees of freedom involved in the excitation. In the context of Kohn-Sham density functional calculations, typical approximations to the exchange-and-correlation functional suffer from a self-interaction error. The Perdew and Zunger self-interaction correction can alleviate this problem, but makes the energy variant to unitary transformations in the occupied orbital space, introducing a large amount of unphysical solutions that do not fully minimize the self-interaction error. An extension of the generalized mode following method is proposed that ensures convergence to the solution minimizing the self-interaction error., Comment: 11 pages, 3 figures

The structure of the solvation shell of aqueous Fe$^{+3}$ ion has been a subject of controversy due to discrepancies between experiments and different levels of theory. We address this issue by performing simulations for a wide range of ion concentrations, using various empirical potential energy functions, as well as density functional theory calculations of selected configurations. The solvation shell undergoes abrupt transitions between two states: an octahedral (OH) state with 6-fold coordination, and a capped trigonal prism (CTP) state with 7-fold coordination. The lifetime of these states is concentration dependent. In dilute $\mathrm{FeCl_3}$ solutions, the lifetime of the two states is similar ($\approx 1$ ns). However, the lifetime of the OH state increases with ion concentration, while that of the CTP state decreases slightly. When a uniform negative background charge is used instead of explicit counterions, the lifetime of the OH state is greatly overestimated. These findings underscore the need for further experimental measurements as well as high-level simulations over sufficiently long timescales and low concentration.

Density functional calculations of Rydberg excited states up to high energy are carried out for several molecules using an approach where the orbitals are variationally optimized by converging on saddle points on the electronic energy surface within a real space grid representation. Remarkably good agreement with experimental estimates of the excitation energy is obtained using the generalized gradient approximation (GGA) functional of Perdew, Burke and Ernzerhof (PBE) when Perdew-Zunger self-interaction correction is applied in combination with complex-valued orbitals. Even without the correction, the PBE functional gives quite good results despite the fact that corresponding Rydberg virtual orbitals have positive energy in the ground state calculation. Results obtained using the TPSS and r2SCAN meta-GGA functionals are also presented, but they do not provide a systematic improvement over the results from the uncorrected PBE functional. The grid representation combined with the projector augmented-wave approach gives a simpler and better representation of the diffuse Rydberg orbitals than a linear combination of atomic orbitals with commonly used basis sets, the latter leading to an overestimation of the excitation energy due to confinement of the excited states., Comment: 34 pages, 8 figures

We review the GPAW open-source Python package for electronic structure calculations. GPAW is based on the projector-augmented wave method and can solve the self-consistent density functional theory (DFT) equations using three different wave-function representations, namely real-space grids, plane waves, and numerical atomic orbitals. The three representations are complementary and mutually independent and can be connected by transformations via the real-space grid. This multi-basis feature renders GPAW highly versatile and unique among similar codes. By virtue of its modular structure, the GPAW code constitutes an ideal platform for implementation of new features and methodologies. Moreover, it is well integrated with the Atomic Simulation Environment (ASE) providing a flexible and dynamic user interface. In addition to ground-state DFT calculations, GPAW supports many-body GW band structures, optical excitations from the Bethe-Salpeter Equation (BSE), variational calculations of excited states in molecules and solids via direct optimization, and real-time propagation of the Kohn-Sham equations within time-dependent DFT. A range of more advanced methods to describe magnetic excitations and non-collinear magnetism in solids are also now available. In addition, GPAW can calculate non-linear optical tensors of solids, charged crystal point defects, and much more. Recently, support of GPU acceleration has been achieved with minor modifications of the GPAW code thanks to the CuPy library. We end the review with an outlook describing some future plans for GPAW.

In reaction path optimization, such as the calculation of a minimum energy path (MEP) between given reactant and product configurations of atoms, it is advantageous to start with an initial guess where close proximity of atoms is avoided and bonds are not unnecessarily broken only to be reformed later. When the configurations of the atoms are described with Cartesian coordinates, a linear interpolation between the endpoints can be problematic, and a better option is provided by the so-called image dependent pair potential (IDPP) approach where interpolated pairwise distances are generated to form an objective function that can be used to generate an improved initial path. When started with a linear interpolation, this method can, however, still lead to unnecessary bond breaking in, for example, reactions where a molecular subgroup undergoes significant rotation. In the method presented here, this problem is addressed by constructing the path gradually, introducing images sequentially starting from the vicinity of the endpoints while the distance between images in the central region is larger. The distribution of images is controlled by systematically scaling the tightness of springs acting between the images until the desired number of images is obtained and they are evenly spaced. This procedure generates an initial path on the IDPP surface, a task that requires insignificant computational effort, as no evaluation of the energy of the system is needed. The calculation of the MEP, typically using electronic structure calculations, is then subsequently carried out in a way that makes efficient use of parallel computing with the nudged elastic band method. Several examples of reactions are given where the linear interpolation IDPP (LI-IDPP) method yields problematic paths with unnecessary bond breaking in some of the intermediate images, while the sequential IDPP (S-IDPP) method yields pat, Comment: 36 pages, 6 figures

Elucidation of the mechanism for optical spin initialization of point defects in solids in the context of quantum applications requires an accurate description of the excited electronic states involved. While variational density functional calculations have been successful in describing the ground state of a great variety of systems, doubts have been expressed in the literature regarding the ability of such calculations to describe electronic excitations of point defects. A direct orbital optimization method is used here to perform time-independent, variational density functional calculations of a prototypical defect, the negatively charged nitrogen-vacancy center in diamond. The calculations include up to 511 atoms subject to periodic boundary conditions and the excited state calculations require similar computational effort as ground state calculations. Contrary to some previous reports, the use of local and semilocal density functionals gives the correct ordering of the low-lying triplet and singlet states, namely ${}^{3}A_2 < {}^{1}E < {}^{1}A_1 < {}^{3}E$. Furthermore, the more advanced meta generalized gradient approximation functionals give results that are in remarkably good agreement with high-level, many-body calculations as well as available experimental estimates, even for the excited singlet state which is often referred to as having multireference character. The lowering of the energy in the triplet excited state as the atom coordinates are optimized in accordance with analytical forces is also close to the experimental estimate and the resulting zero-phonon line triplet excitation energy is underestimated by only 0.15 eV. The approach used here is found to be a promising tool for studying electronic excitations of point defects in, for example, systems relevant for quantum technologies.

The effect of solvent on the free energy of reaction intermediates adsorbed on electrocatalyst surfaces can significantly change the thermochemical overpotential, but accurate calculations of this are challenging. Here, we present computational estimates of the solvation energy for reaction intermediates in oxygen reduction reaction (ORR) on a B-doped graphene (BG) model system where the overpotential is found to reduce by up to 0.6 V due to solvation. BG is experimentally reported to be an active ORR catalyst but recent computational estimates using state-of-the-art hybrid density functionals in the absence of solvation effects have indicated low activity. To test whether the inclusion of explicit solvation can bring the calculated activity estimates closer to the experimental reports, up to 4 layers of water molecules are included in the simulations reported here. The calculations are based on classical molecular dynamics and local minimization of energy using atomic forces evaluated from electron density functional theory. Data sets are obtained from regular and coarse-grained dynamics, as well as local minimization of structures resampled from dynamics simulations. The results differ greatly depending on the method used and the solvation energy estimates and are deemed untrustworthy. It is concluded that a significantly larger number of water molecules is required to obtain converged results for the solvation energy. As the present system includes up to 139 atoms, it already strains the limits of computational feasibility, so this points to the need for a hybrid simulation approach where efficient simulations of much larger number of solvent molecules is carried out using a lower level of theory while retaining the higher level of theory for the reacting molecules as well as their near neighbors and the catalyst.

Variational calculations of excited electronic states are carried out by finding saddle points on the surface that describes how the energy of the system varies as a function of the electronic degrees of freedom. This approach has several advantages over commonly used methods especially in the context of density functional calculations, as collapse to the ground state is avoided and yet, the orbitals are variationally optimized for the excited state. This optimization makes it possible to describe excitations with large charge transfer where calculations based on ground state orbitals are problematic, as in linear response time-dependent density functional theory. A generalized mode following method is presented where an $n^{\text{th}}$-order saddle point is found by inverting the components of the gradient in the direction of the eigenvectors of the $n$ lowest eigenvalues of the electronic Hessian matrix. This approach has the distinct advantage of following a chosen excited state through atomic configurations where the symmetry of the single determinant wave function is broken, as demonstrated in calculations of potential energy curves for nuclear motion in the ethylene and dihydrogen molecules. The method is implemented using a generalized Davidson algorithm and an exponential transformation for updating the orbitals within a generalized gradient approximation of the energy functional. Convergence is found to be more robust than for a direct optimization approach previously shown to outperform standard self-consistent field approaches, as illustrated here for charge transfer excitations in nitrobenzene and N-phenylpyrrole, involving calculations of $4^{\text{th}}$- and $6^{\text{th}}$-order saddle points, respectively. Finally, calculations of a diplatinum and silver complex are presented, illustrating the applicability of the method to excited state energy curves of large molecules., Comment: 57 pages, 12 figures, submitted to the Journal of Chemical Theory and Computation

A Gaussian approximation machine learning interatomic potential for platinum is presented. It has been trained on DFT data computed for bulk, surfaces and nanostructured platinum, in particular nanoparticles. Across the range of tested properties, which include bulk elasticity, surface energetics and nanoparticle stability, this potential shows excellent transferability and agreement with DFT, providing state-of-the-art accuracy at low computational cost. We showcase the possibilities for modeling of Pt systems enabled by this potential with two examples: the pressure-temperature phase diagram of Pt calculated using nested sampling and a study of the spontaneous crystallization of a large Pt nanoparticle based on classical dynamics simulations over several nanoseconds.

Magnetic hopfions are string-like three-dimensional topological solitons, characterised by the Hopf invariant. They serve as a fundamental prototype for three-dimensional magnetic quasi-particles and are an inspiration for novel device concepts in the field of spintronics. Based on a micromagnetic model and without considering temperature, the existence of such hopfions has been predicted in certain magnets with competing exchange interactions. However, physical realisation of freely moving hopfions in bulk magnets have so far been elusive. Here, we consider an effective Heisenberg model with competing exchange interactions and study the stability of small toroidal hopfions with Hopf number $Q_\text{H}=1$ by finding first-order saddle points on the energy surface representing the transition state for the decay of hopfions via the formation of two coupled Bloch points. We combine the geodesic nudged elastic band method and an adapted implementation of the dimer method to resolve the sharp energy profile of the reaction path near the saddle point. Our analysis reveals that the energy barrier can reach substantial height and is largely determined by the size of the hopfion relative to the lattice constant., Comment: 13 pages, 8 figures

Transition metal dichalcogenides are cheap and earth-abundant candidates for the replacement of precious metals as catalyst materials. Experimental measurements of the hydrogen evolution reaction (HER), for example, have demonstrated significant electrocatalytic activity of MoS2 but there is large variation depending on preparation method. In order to gain information about the mechanism and active sites for HER, we have carried out calculations of the reaction and activation energy for HER at the transition metal doped basal plane of MoS2 under electrochemical conditions, i.e. including applied electrode potential and solvent effects. The calculations are based on identifying the relevant saddle points on the energy surface obtained from density functional theory within the generalized gradient approximation, and the information on energetics is used to construct voltage-dependent volcano plots. Doping with 3d-metal atoms as well as Pt is found to enhance hydrogen adsorption onto the basal plane by introducing electronic states within the band gap, and in some cases (Co, Ni, Cu, Pt) significant local symmetry breaking. The Volmer-Heyrovsky mechanism is found to be most likely and the associated energetics show considerable dopant and voltage-dependence. While the binding free energy of hydrogen can be tuned to be seemingly favorable for HER, the calculated activation energy turns out to be significant, at least 0.7 eV at a voltage of -0.5 V vs. SHE, indicating low catalytic activity of the doped basal plane. This suggests that other sites are responsible for the experimental activity, possibly edges or basal plane defects.

Icosahedral Au clusters with three and four shells of atoms are found to deviate significantly from the commonly assumed Mackay structures. By introducing additional atoms in the surface shell and creating a vacancy in the center of the cluster, the calculated energy per atom can be lowered significantly, according to several different descriptions of the interatomic interaction. Analogous icosahedral structures with five and six shells of atoms are generated using the same structural motifs and are similarly found to be more stable than Mackay icosahedra. The lowest energy per atom is obtained with clusters containing 310, 564, 928 and 1426 atoms, as compared with the commonly assumed magic numbers of 309, 561, 923 and 1415. Some of the vertices in the optimized clusters have a hexagonal ring of atoms, rather than a pentagon, with the vertex atom missing. An inner shell atom in some cases moves outwards by more than an \AA{}ngstr\"om into the surface shell at the vertex site. This feature, as well as the wide distribution of nearest-neighbor distances in the surface layer, can strongly influence the catalytic properties of icosahedral clusters. The structural optimization is initially carried out using the GOUST method with atomic forces estimated with the EMT empirical potential function, but the atomic coordinates are then refined by minimization using electron density functional theory (DFT) or Gaussian approximation potential (GAP). A single energy barrier is found to separate the Mackay icosahedron from a lower energy structure where a string of atoms moves outwards in a concerted manner from the center so as to create a central vacancy while placing an additional atom in the surface shell.

Theoretical studies of photochemical processes require a description of the energy surfaces of excited electronic states, especially near degeneracies, where transitions between states are most likely. Systems relevant to photochemical applications are typically too large for high-level multireference methods, and while time-dependent density functional theory (TDDFT) is efficient, it can fail to provide the required accuracy. A variational, time-independent density functional approach is applied to the twisting of the double bond and pyramidal distortion in ethylene, the quintessential model for photochemical studies. By allowing for symmetry breaking, the calculated energy surfaces exhibit the correct topology around the twisted-pyramidalized conical intersection even when using a semilocal functional approximation, and by including explicit self-interaction correction, the torsional energy curves are in close agreement with published multireference results. The findings of the present work point to the possibility of using a single determinant time-independent density functional approach to simulate nonadiabatic dynamics, even for large systems where multireference methods are impractical and TDDFT is often not accurate enough., Comment: 53 pages, 14 figures, submitted to the Journal of Physical Chemistry Letters

Calculations of the diffusion of a Au adatom on the dimer reconstructed Si(100)-2x1 surface reveal an interesting mechanism that differs significantly from a direct path between optimal binding sites, which are located in between dimer rows. Instead, the active diffusion mechanism involves promotion of the adatom to higher energy sites on top of a dimer row and then fast migration along the row, visiting ca. a hundred sites at room temperature, before falling back down into an optimal binding site. This top-of-row mechanism becomes more important the lower the temperature is. The calculations are carried out by finding minimum energy paths on the energy surface obtained from density functional theory within the PBEsol functional approximation followed by kinetic Monte Carlo simulations of the diffusion over a range of temperature from 200 K to 900 K. While the activation energy for the direct diffusion mechanism is calculated to be 0.84 eV, the effective activation energy for the indirect mechanism is on average 0.56 eV.

The activity and stability of a platinum nanoparticle (NP) is not only affected by its size but additionally depends on its shape. To this end, simulations can identify structure-property relationships to make a priori decisions on the most promising structures. While activity is routinely probed by electronic structure calculations on simplified surface models, modeling the stability of NP model systems in electrochemical reactions is challenging due to the long timescale of relevant processes such as oxidation beyond the point of reversibility. In this work, a routine for simulating electrocatalyst stability is presented. The procedure is referred to as GREG after its main ingredients - a grand-canonical simulation approach using reactive force fields to model electrochemical reactions as a function of the galvanic cell potential. The GREG routine is applied to study the oxidation of 3 nm octahedral, cubic, dodecahedral, cuboctahedral, spherical, and tetrahexahedral platinum NPs. The oxidation process is analyzed using adsorption isobars as well as interaction energy heat maps that provide the basis for constructing electrochemical phase diagrams. Onset potentials for surface oxidation increase in the sequence cube ~= dodecahedron <= octahedron <= tetrahexahdron < sphere < cuboctahedron, establishing a relationship between oxidation behavior and surface facet structure. The electrochemical results are rationalized using structural and electronic analysis., Comment: Submitted to J. Phys. Chem. C

Experimental studies of the oxygen reduction reaction (ORR) at nitrogen doped graphene electrodes have reported a remarkably low overpotential, on the order of 0.5 V, similar to Pt based electrodes. Theoretical calculations using density functional theory have lent support for this claim. However, other measurements have indicated that transition metal impurities are actually responsible for the ORR activity, thereby raising questions about the reliability of both the experiments and the calculations. In order to assess the accuracy of the theoretical calculations, various generalized gradient approximation (GGA), meta-GGA and hybrid functionals are employed here and calibrated against high-level wave function based coupled cluster calculations (CCSD(T)) of the overpotential as well as self-interaction corrected density functional calculations and published quantum Monte Carlo calculations of O adatom binding to graphene. The PBE0 and HSE06 hybrid functionals are found to give more accurate results than the GGA and meta-GGA functionals, as would be expected, and for low dopant concentration, 3.1%, the overpotential is calculated to be 1.0 V. The GGA and meta-GGA functionals give a lower estimate by as much as 0.4 V. When the dopant concentration is doubled, the overpotential calculated with hybrid functionals drops, while it increases in GGA functional calculations. The opposite trends result from different potential determining steps, the *OOH species being of central importance in the hybrid functional calculations while the reduction of *O determines the overpotential obtained in GGA and meta-GGA calculations. The results presented here are mainly based on calculations of periodic representations of the system, but a comparison is also made with molecular flake models which are found to give erratic results., Comment: To be published in Journal of Chemical Theory and Computation

Aberration-corrected scanning electron microscopy (AC-STEM) can provide valuable information on the atomic structure of nanoclusters, an essential input for gaining an understanding of their physical and chemical properties. A systematic method is presented here for the extraction of atom coordinates from an AC-STEM image in a way that is general enough to be applicable to irregular structures. The two-dimensional information from the image is complemented with an approximate description of the atomic interactions so as to construct a three-dimensional structure and, at a final stage, the structure is refined using electron density functional theory (DFT) calculations. The method is applied to an AC-STEM image of Au55. Analysis of the local structure shows that the cluster is a combination of a part with icosahedral structure elements and a part with local atomic arrangement characteristic of crystal packing, including a segment of a flat surface facet. The energy landscape of the cluster is explored in calculations of minimum energy paths between the optimal fit structure and other candidates generated in the analysis. This reveals low energy barriers for conformational changes, showing that such transitions can occur on laboratory timescale even at room temperature and lead to large changes in the AC-STEM image. The paths furthermore reveal additional cluster configurations, some with lower DFT energy and providing nearly as good fit to the experimental image.

The scattering-adapted flexible inner region ensemble separator (SAFIRES) is a partitioning scheme designed to divide a simulation cell into two regions to be treated with different computational methodologies. SAFIRES prevents particles from crossing between regions and resolves boundary events through elastic collisions of the particles mediated by the boundary, conserving energy and momenta. A multiple-time-step propagation algorithm is introduced where the time step is scaled automatically to identify the moment a collision occurs. If the length of the time step is kept constant, the new propagator reduces to a regular algorithm for Langevin dynamics, and to the velocity Verlet algorithm for classical dynamics if the friction coefficient is set to zero. SAFIRES constitutes the exact limit of the premise behind boundary-based methods such as FIRES, BEST, and BCC which take advantage of the indistinguishability of molecules on opposite sides of the separator. It gives correct average ensemble statistics despite the introduction of an ensemble separator. SAFIRES is tested in simulations where the molecules on the two sides are treated in the same way, for a Lennard-Jones (LJ) liquid and a LJ liquid in contact with a surface, as well as for liquid modelling simulations using the TIP4P force field. Simulations using SAFIRES are shown to reproduce the unconstrained reference simulations without significant deviations., Comment: To be submitted to Journal of Chemical Theory and Computation

A direct orbital optimization method is presented for density functional calculations of excited electronic states using either a real space grid or a plane wave basis set. The method is variational, provides atomic forces in the excited states, and can be applied to Kohn-Sham (KS) functionals as well as orbital-density dependent functionals (ODD) including explicit self-interaction correction. The implementation for KS functionals involves two nested loops: (1) An inner loop for finding a stationary point in a subspace spanned by the occupied and a few virtual orbitals corresponding to the excited state; (2) an outer loop for minimizing the energy in a tangential direction in the space of the orbitals. For ODD functionals, a third loop is used to find the unitary transformation that minimizes the energy functional among occupied orbitals only. Combined with the maximum overlap method, the algorithm converges in challenging cases where conventional self-consistent field algorithms tend to fail. The benchmark tests presented include two charge-transfer excitations in nitrobenzene and an excitation of CO to degenerate $\pi^\ast$ orbitals where the importance of complex orbitals is illustrated. An application of the method to several metal-to-ligand charge-transfer and metal-centred excited states of an Fe$^{\rm II}$ photosensitizer complex is described and the results compared to reported experimental estimates. The method is also used to study the effect of Perdew-Zunger self-interaction correction on valence and Rydberg excited states of several molecules, both singlet and triplet states, and the performance compared to semilocal and hybrid functionals., Comment: 68 pages, 15 figures, including supporting information

Qualitatively incorrect results are obtained for the Mn dimer in density functional theory calculations using the generalized gradient approximation (GGA) and similar results are obtained from local density and meta-GGA functionals. The coupling is predicted to be ferromagnetic rather than antiferromagnetic and the bond between the atoms is predicted to be an order of magnitude too strong and about an {\AA}ngstr{\o}m too short. Explicit, self-interaction correction (SIC) applied to a commonly used GGA energy functional, however, provides close agreement with both experimental data and high-level, multi-reference wave function calculations. These results show that the failure is not due to strong correlation but rather the single electron self-interaction that is necessarily introduced in estimates of the classical Coulomb and exchange-correlation energy when only the total electron density is used as input. The corrected functional depends explicitly on the orbital densities and can, therefore, avoid the introduction of self-Coulomb interaction. The error arises because of over-stabilization of bonding $d$-states in the minority spin channel resulting from an overestimate of the $d$-electron self-interaction in the semi-local exchange-correlation functionals. Since the computational effort in the self-interaction corrected calculations scales with system size in the same way as for regular semi-local functional calculations, this approach provides a way to calculate properties of Mn nanoclusters as well as biomolecules and extended solids where Mn dimers and larger cluster are present, while multi-reference wave function calculations can only be applied to small systems., Comment: 20 pages, 4 figures

The energy minimization involved in density functional calculations of electronic systems can be carried out using an exponential transformation that preserves the orthonormality of the orbitals. The energy of the system is then represented as a function of the elements of a skew-Hermitian matrix that can be optimized directly using unconstrained minimization methods. An implementation based on the limited memory Broyden-Fletcher-Goldfarb-Shanno approach with inexact line search and a preconditioner is presented and the performance compared with that of the commonly used self-consistent field approach. Results are presented for the G2 set of 148 molecules, liquid water configurations with up to 576 molecules and some insulating crystals. A general preconditioner is presented that is applicable to systems with fractional orbital occupation as is, for example, needed in the k-point sampling for periodic systems. This exponential transformation direct minimization approach is found to outperform the standard implementation of the self-consistent field approach in that all the calculations converge with the same set of parameter values and it requires less computational effort on average. The formulation of the exponential transformation and the gradients of the energy presented here are quite general and can be applied to energy functionals that are not unitary invariant such as self-interaction corrected functionals.

A potential function is presented for describing a system of flexible H$_2$O molecules based on the single center multipole expansion (SCME) of the electrostatic interaction. The model, referred to as SCME/f, includes the variation of the molecular quadrupole moment as well as the dipole moment with changes in bond length and angle so as to reproduce results of high level electronic structure calculations. The multipole expansion also includes fixed octupole and hexadecapole moments, as well as anisotropic dipole-dipole, dipole-quadrupole and quadrupole-quadrupole polarizability tensors. The model contains five adjustable parameters related to the repulsive interaction and damping functions in the electrostatic and dispersion interactions. Their values are adjusted to reproduce the lowest energy isomers of small clusters, (H$_2$O)$_n$ with $n=2-6$, as well as measured properties of the ice Ih crystal. Subsequent calculations of the energy difference between the various isomer configurations of the clusters show that SCME/f gives good agreement with results of electronic structure calculations and represents a significant improvement over the previously presented rigid SCME potential function. Analysis of the vibrational frequencies of the clusters and structural properties of ice Ih crystal show the importance of accurately describing the variation of the quadrupole moment with molecular structure.

Olivine is involved in many natural reactions and industrial reactions as a catalyst. The catalytic ability is highly possible rely on the Fe$^{2+}$ in olivine. We use density functional theory calculation and thermodynamics to investigate the site preference of Fe atom in olivine which composition from iron-rich to iron-poor and its surfaces. The Fe$^{2+}$ always shows its high spin (quintet) state which has larger ion radius than Mg$^{2+}$ in olivine crystal and surfaces. The Fe$^{2+}$ inside the surface slab prefers the smaller M1 site than M2 site by enlarging the metal-oxygen octahedra when occupied the metal site as in the bulk system. Energy contribution of entropies accumulation caused temperature raise stops this preference at the temperature where a cation order-disorder distribution energy crossover happen in olivine. Surface exposed site provide Fe$^{2+}$ large space due its unsaturated nature. This lead a higher level of preference of Fe$^{2+}$ to the surface site than any metal site inside the crystal no matter M1 or M2 site is exposed. This indicate the Fe$^{2+}$ in the bulk system can diffuse to a metal site exposed on the surface driven by the energy difference. Many reactions can use the on surface Fe$^{2+}$ as a catalyst because of the active chemical behavior of Fe. Meanwhile this energetics preference should be considered in the future model to explain the natural observed zoning olivine have a high Fe edge and low Fe center. These microscopic understanding can be essential to many olivine related geochemical and astrochemical reactions., Comment: 8 figures

The development of variational density functional theory approaches to excited electronic states is impeded by limitations of the commonly used self-consistent field (SCF) procedure. A method based on a direct optimization approach as well as the maximum overlap method is presented and the performance compared with previously proposed SCF strategies. Excited-state solutions correspond to saddle points of the energy as a function of the electronic degrees of freedom. The approach presented here makes use of a preconditioner determined with the help of the maximum overlap method to guide the convergence on a target nth-order saddle point. The method is found to be more robust and to converge faster than previously proposed SCF approaches for a set of 89 excited states of molecules. A limited-memory formulation of the symmetric rank-one method for updating the inverse Hessian is found to give the best performance. A conical intersection for the carbon monoxide molecule is calculated without resorting to fractional occupation numbers. Calculations on excited states of the hydrogen atom and a doubly excited state of the dihydrogen molecule using a self-interaction corrected functional are presented. For these systems, the self-interaction correction is found to improve the accuracy of density functional calculations of excited states.

Calculations of skyrmions in antiferromagnets (AFMs) are presented, and their properties compared with skyrmions in corresponding ferromagnets (FMs). The rates of skyrmion collapse and escape through the boundary of a track, as well as the binding to and collapse at a non-magnetic impurity, are calculated as a function of applied magnetic field. The activation energy for skyrmion annihilation is the same in AFMs and corresponding FMs in the absence of an applied magnetic field. The pre-exponential factor in the Arrhenius rate law is, however, different because skyrmion dynamics is different in the two systems. An applied magnetic field has opposite effects on skyrmions in the two types of materials. In AFMs the rate of collapse of skyrmions as well as the rate of escape through the edge of a magnetic strip decreases slightly with increasing field, while these rates increase strongly for a skyrmion in the corresponding FMs when the field is directed antiparallel to the magnetization in the center of the skyrmion. A non-magnetic impurity is less likely to trap a skyrmion in AFMs especially in the presence of a magnetic field. This, together with the established fact that a spin polarized current moves skyrmions in AFMs in the direction of the current, while in FMs skyrmions move at an angle to the current, demonstrates that skyrmions in AFMs have several advantageous properties over skyrmions in FMs for memory and spintronic devices.

Efficient algorithms for the calculation of minimum energy paths of magnetic transitions are implemented within the geodesic nudged elastic band (GNEB) approach. While an objective function is not available for GNEB and a traditional line search can, therefore, not be performed, the use of limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) and conjugate gradient algorithms in conjunction with orthogonal spin optimization (OSO) approach is shown to greatly outperform the previously used velocity projection and dissipative Landau-Lifschitz dynamics optimization methods. The implementation makes use of energy weighted springs for the distribution of the discretization points along the path and this is found to improve performance significantly. The various methods are applied to several test problems using a Heisenberg-type Hamiltonian, extended in some cases to include Dzyaloshinskii-Moriya and exchange interactions beyond nearest neighbours. Minimum energy paths are found for magnetization reversals in a nano-island, collapse of skyrmions in two-dimensional layers and annihilation of a chiral bobber near the surface of a three-dimensional magnet. The LBFGS-OSO method is found to outperform the dynamics based approaches by up to a factor of 8 in some cases.

Locally, the atomic structure in well annealed amorphous silicon appears similar to that of crystalline silicon. We address here the question whether a point defect, specifically a vacancy, in amorphous silicon also resembles that in the crystal. From density functional theory calculations of a large number of nearly defect free configurations, relaxed after an atom has been removed, we conclude that there is little similarity. The analysis is based on formation energy, relaxation energy, bond lengths, bond angles, Vorono\"i volume, coordination, atomic charge and electronic gap states. All these quantities span a large and continuous range in amorphous silicon and while the removal of an atom leads to the formation of one to two bond defects and to a lowering of the local atomic density, the relaxation of the bonding network is highly effective, and the signature of the vacancy generally unlike that of a vacancy in the crystal., Comment: 9 pages, 8 figures

Magnetic skyrmions are nano-scale magnetic states that could be used in various spintronics devices. A central issue is the mechanism and rate of various possible annihilation processes and the lifetime of metastable skyrmions. While most studies have focused on classical over-the-barrier mechanism for annihilation, it is also possible that quantum mechanical tunneling through the energy barrier takes place. Calculations of the lifetime of magnetic skyrmions in a two-dimensional lattice are presented and the rate of tunneling compared with the classical annihilation rate. A remarkably strong variation in the crossover temperature and the lifetime of the skyrmion is found as a function of the values of parameters in the extended Heisenberg Hamiltonian, i.e. the out-of-plane anisotropy, Dzyaloshinskii-Moriya interaction (DMI) and applied magnetic field. Materials parameters and conditions are identified where the onset of tunneling could be observed on a laboratory time scale. In particular, it is predicted that skyrmion tunneling could be observed in the PdFe/Ir(111) system when an external magnetic field on the order of 6 T is applied.

An algorithm for the minimization of the energy of magnetic systems is presented and applied to the analysis of thermal configurations of a ferromagnet to identify inherent structures, i.e. the nearest local energy minima, as a function of temperature. Over a rather narrow temperature interval, skyrmions appear and reach a high temperature limit for the skyrmion density. In addition, the performance of the algorithm is further demonstrated in a self-consistent field calculation of a skyrmion in an itinerant magnet. The algorithm is based on a geometric approach in which the curvature of the spherical domain is taken into account and as a result the length of the magnetic moments is preserved in every iteration. In the limit of infinitesimal rotations, the minimization path coincides with that obtained using damped spin dynamics while the use of limited-memory quasi-newton minimization algorithms, such as the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm, significantly accelerates the convergence.

The \textit{Spirit} framework is designed for atomic scale spin simulations of magnetic systems of arbitrary geometry and magnetic structure, providing a graphical user interface with powerful visualizations and an easy to use scripting interface. An extended Heisenberg type spin-lattice Hamiltonian including competing exchange interactions between neighbors at arbitrary distance, higher-order exchange, Dzyaloshinskii-Moriya and dipole-dipole interactions is used to describe the energetics of a system of classical spins localised at atom positions. A variety of common simulations methods are implemented including Monte Carlo and various time evolution algorithms based on the Landau-Lifshitz-Gilbert equation of motion, which can be used to determine static ground state and metastable spin configurations, sample equilibrium and finite temperature thermodynamical properties of magnetic materials and nanostructures or calculate dynamical trajectories including spin torques induced by stochastic temperature or electric current. Methods for finding the mechanism and rate of thermally assisted transitions include the geodesic nudged elastic band method, which can be applied when both initial and final states are specified, and the minimum mode following method when only the initial state is given. The lifetime of magnetic states and rate of transitions can be evaluated within the harmonic approximation of transition-state theory. The framework offers performant CPU and GPU parallelizations. All methods are verified and applications to several systems, such as vortices, domain walls, skyrmions and bobbers are described.

Bridgmanite, a high temperature and pressure form of $MgSiO_3$, is believed to be Earth's most abundant mineral and responsible for the observed seismic anisotropy in the mantle. Little is known about surfaces of bridgmanite but knowledge of the most stable surface terminations is important for understanding various geochemical processes as well as likely slip planes. A density functional theory based thermodynamic approach is used here to establish the range of stability of bridgmanite as well as possible termination structures of the (001), (010), (100) and (011) surfaces as a function of the chemical potential of oxygen and magnesium. The results presented provide a basis for further theoretical studies of the chemical processes on bridgmanite surfaces in the Earth's mantle and slip plane analysis., Comment: 4 Pages,4 figures

The stability of possible termination structures for the (010) surface of forsterite, $ Mg_2SiO_4 $, is studied using a density functional theory (DFT) based thermodynamic approach. The DFT calculations are used to estimate the surface Gibbs free energy of various surface structures and compare their stability as a function of the chemical environment. Among 9 possible terminations, the SiO-II, M2, O-II terminations are found to be most stable as conditions range from Mg-poor to Mg-rich. This relative stability order remains the same at an elevated temperature. The surface phase diagram obtained provides ground for further theoretical studies of chemical processes on forsterite surfaces in terrestrial planets.

Various transitions that a magnetic skyrmion can undergo are found in calculations using a method for climbing up the energy surface and converging onto first order saddle points. In addition to collapse and escape through a boundary, the method identifies a transition where the skyrmion divides and forms two skyrmions. The activation energy for this duplication process can be similar to that of collapse and escape. A tilting of the external magnetic field for a certain time interval is found to induce the duplication process in a dynamical simulation. Such a process could turn out to be an important avenue for the creation of skyrmions in future magnetic devices.

We present an efficient implementation of a surface Green's-function method for atomistic modeling of surfaces within the framework of density functional theory using a pseudopotential localized basis set approach. In this method, the system is described as a truly semi-infinite solid with a surface region coupled to an electron reservoir, thereby overcoming several fundamental drawbacks of the traditional slab approach. The versatility of the method is demonstrated with several applications to surface physics and chemistry problems that are inherently difficult to address properly with the slab method, including metal work function calculations, band alignment in thin-film semiconductor heterostructures, surface states in metals and topological insulators, and surfaces in external electrical fields. Results obtained with the surface Green's-function method are compared to experimental measurements and slab calculations to demonstrate the accuracy of the approach.

The stability of skyrmions in various environments is estimated by analyzing the multidimensional surface describing the energy of the system as a function of the directions of the magnetic moments in the system. The energy is given by a Heisenberg-like Hamiltonian that includes Dzyaloshinskii-Moriya interaction, anisotropy and external magnetic field. Local minima on this surface correspond to the ferromagnetic and skyrmion states. Minimum energy paths (MEP) between the minima are calculated using the geodesic nudged elastic band method. The maximum energy along an MEP corresponds to a first order saddle point on the energy surface and gives an estimate of the activation energy for the magnetic transition, such as creation and annihilation of a skyrmion. The pre-exponential factor in the Arrhenius law for the rate, the so-called attempt frequency, is estimated within harmonic transition state theory where the eigenvalues of the Hessian at the saddle point and the local minima are used to characterize the shape of the energy surface. For some degrees of freedom, so-called 'zero modes', the energy of the system remains invariant. They need to be treated separately and give rise to temperature dependence of the attempt frequency. As an example application of this general theory, the lifetime of a skyrmion in a track of finite width for a PdFe overlayer on a Ir(111) substrate is calculated as a function of track width and external magnetic field. Also, the effect of non-magnetic impurities is studied. Various MEPs for annihilation inside a track, via the boundary of a track and at an impurity are presented. The attempt frequency as well as the activation energy has been calculated for each mechanism to estimate the transition rate as a function of temperature.

Minimum energy paths for transitions such as atomic and/or spin rearrangements in thermalized systems are the transition paths of largest statistical weight. Such paths are frequently calculated using the nudged elastic band method, where an initial path is iteratively shifted to the nearest minimum energy path. The computational effort can be large, especially when ab initio or electron density functional calculations are used to evaluate the energy and atomic forces. Here, we show how the number of such evaluations can be reduced by an order of magnitude using a Gaussian process regression approach where an approximate energy surface is generated and refined in each iteration. When the goal is to evaluate the transition rate within harmonic transition state theory, the evaluation of the Hessian matrix at the initial and final state minima can be carried out beforehand and used as input in the minimum energy path calculation, thereby improving stability and reducing the number of iterations needed for convergence. A Gaussian process model also provides an uncertainty estimate for the approximate energy surface, and this can be used to focus the calculations on the lesser-known part of the path, thereby reducing the number of needed energy and force evaluations to a half in the present calculations. The methodology is illustrated using the two-dimensional M\"uller-Brown potential surface and performance assessed on an established benchmark involving 13 rearrangement transitions of a heptamer island on a solid surface.

The stability of magnetic skyrmions against thermal fluctuations and external perturbations is investigated within the framework of harmonic transition state theory for magnetic degrees of freedom. The influence of confined geometry and atomic scale non-magnetic defects on the skyrmion lifetime is estimated. It is shown that a skyrmion on a track has lower activation energy for annihilation and higher energy for nucleation if the size of the skyrmion is comparable with the width of the track. Two mechanisms of skyrmion annihilation are considered: inside the track and escape through the boundary. For both mechanisms, the dependence of activation energy on the track width is calculated. Non-magnetic defects are found to localize skyrmions in their neighborhood and strongly decrease the activation energy for creation and annihilation. This is in agreement with experimental measurements that have found nucleation of skyrmions in presence of spin-polarized current preferably occurring near structural defects.

The calculation of minimum energy paths for transitions such as atomic and/or spin re-arrangements is an important task in many contexts and can often be used to determine the mechanism and rate of transitions. An important challenge is to reduce the computational effort in such calculations, especially when ab initio or electron density functional calculations are used to evaluate the energy since they can require large computational effort. Gaussian process regression is used here to reduce significantly the number of energy evaluations needed to find minimum energy paths of atomic rearrangements. By using results of previous calculations to construct an approximate energy surface and then converge to the minimum energy path on that surface in each Gaussian process iteration, the number of energy evaluations is reduced significantly as compared with regular nudged elastic band calculations. For a test problem involving rearrangements of a heptamer island on a crystal surface, the number of energy evaluations is reduced to less than a fifth. The scaling of the computational effort with the number of degrees of freedom as well as various possible further improvements to this approach are discussed.