• 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.
  • 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 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.
  • 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.
  • The theory for the generation of Wannier functions within the generalized Pipek--Mezey approach [Lehtola, S.; J\'onsson, H. J. Chem. Theory Comput. 2014, 10, 642] is presented and an implementation thereof is described. Results are presented for systems with periodicity in one, two and three dimensions as well as isolated molecules. The generalized Pipek--Mezey Wannier functions (PMWF) are highly localized orbitals consistent with chemical intuition where a distinction is maintained between {\sigma}- and {\pi}-orbitals. The PMWF method is compared with the so-called maximally localized Wannier functions (MLWF) that are frequently used for the analysis of condensed matter calculations. Whereas PMWFs maximize the localization criterion of Pipek and Mezey, MLWFs maximize that of Foster and Boys and have the disadvantage of mixing {\sigma}- and {\pi}-orbitals in many cases. The PMWF orbitals turn out to be as localized as the MLWF orbitals as evidenced by cross-comparison of the values of the PMWF and MLWF objective functions for the two types of orbitals. Our implementation in the atomic simulation environment (ASE) is compatible with various representations of the wave function, including real-space grids, plane waves and linear combinations of atomic orbitals. The projector augmented wave formalism for the representation of atomic core electrons is also supported. Results of calculations with the GPAW software are described here, but our implementation can also use output from other electronic structure software such as ABINIT, NWChem and VASP.
  • The recrossing correction to the transition state theory estimate of a thermal rate can be difficult to calculate when the energy barrier is flat. This problem arises, for example, in polymer escape if the polymer is long enough to stretch between the initial and final state energy wells while the polymer beads undergo diffusive motion back and forth over the barrier. We present an efficient method for evaluating the correction factor by constructing a sequence of hyperplanes starting at the transition state and calculating the probability that the system advances from one hyperplane to another towards the product. This is analogous to what is done in forward flux sampling except that there the hyperplane sequence starts at the initial state. The method is applied to the escape of polymers with up to 64 beads from a potential well. For high temperature, the results are compared with direct Langevin dynamics simulations as well as forward flux sampling and excellent agreement between the three rate estimates is found. The use of a sequence of hyperplanes in the evaluation of the recrossing correction speeds up the calculation by an order of magnitude as compared with the traditional approach. As the temperature is lowered, the direct Langevin dynamics simulations as well as the forward flux simulations become computationally too demanding, while the harmonic transition state theory estimate corrected for recrossings can be calculated without significant increase in the computational effort.
  • Global optimization of transition paths in complex atomic scale systems is addressed in the context of misfit dislocation formation in a strained Ge film on Si(001). Such paths contain multiple intermediate minima connected by minimum energy paths on the energy surface emerging from the atomic interactions in the system. The challenge is to find which intermediate states to include and to construct a path going through these intermediates in such a way that the overall activation energy for the transition is minimal. In the numerical approach presented here, intermediate minima are constructed by heredity transformations of known minimum energy structures and by identifying local minima in minimum energy paths calculated using a modified version of the nudged elastic band method. Several mechanisms for the formation of a 90{\deg} misfit dislocation at the Ge-Si interface are identified when this method is used to construct transition paths connecting a homogeneously strained Ge film and a film containing a misfit dislocation. One of these mechanisms which has not been reported in the literature is detailed. The activation energy for this path is calculated to be 26% smaller than the activation energy for half loop formation of a full, isolated 60{\deg} dislocation. An extension of the common neighbor analysis method involving characterization of the geometrical arrangement of second nearest neighbors is used to identify and visualize the dislocations and stacking faults.
  • A possible mechanism for the formation of a 90{\deg} misfit dislocation at the Ge/Si(001) interface through homogeneous nucleation is identified from atomic scale calculations where a minimum energy path connecting the coherent epitaxial state and a final state with a 90{\deg} misfit dislocation is found using the nudged elastic band method. The initial path is generated using a repulsive bias activation procedure in a model system including 75000 atoms. The energy along the path exhibits two maxima in the energy. The first maximum occurs as a 60{\deg} dislocation nucleates. The intermediate minimum corresponds to an extended 60{\deg} dislocation. The subsequent energy maximum occurs as a second 60{\deg} dislocation nucleates in a complementary, mirror glide plane, simultaneously starting from the surface and from the first 60{\deg} dislocation. The activation energy of the nucleation of the second dislocation is 30% lower than that of the first one showing that the formation of the second 60{\deg} dislocation is aided by the presence of the first one. The simulations represent a step towards unraveling the formation mechanism of 90{\deg} dislocations, an important issue in the design of growth procedures for strain released Ge overlayers on Si(100) surfaces, and more generally illustrate an approach that can be used to gain insight into the mechanism of complex nucleation paths of extended defects in solids.
  • The temperature dependence of the response of a magnetic system to an applied field can be understood qualitatively by considering variations in the energy surface characterizing the system and estimated quantitatively with rate theory. In the system analysed here, Fe/Sm-Co spring magnet, the width of the hysteresis loop is reduced to a half when temperature is raised from 25~K to 300~K. This narrowing can be explained and reproduced quantitatively without invoking temperature dependence of model parameters as has typically been done in previous data analysis. The applied magnetic field lowers the energy barrier for reorientation of the magnetization but thermal activation brings the system over the barrier. A 2-dimensional representation of the energy surface is developed and used to gain insight into the transition mechanism and to demonstrate how the applied field alters the transition path. Our results show the importance of explicitly including the effect of thermal activation when interpreting experiments involving the manipulation of magnetic systems at finite temperature.
  • In finite systems, such as nanoparticles and gas-phase molecules, calculations of minimum energy paths (MEPs) connecting initial and final states of transitions as well as searches for saddle points are complicated by the presence of external degrees of freedom, such as overall translation and rotation. A method based on quaternion algebra for removing the external degrees of freedom is described here and applied in calculations using two commonly used methods: the nudged elastic band (NEB) method for MEPs and the DIMER method for finding the minimum mode in minimum mode following searches of first-order saddle points. With the quaternion approach, fewer images in the NEB are needed to represent MEPs accurately. In both NEB and DIMER calculations of finite systems, the number of iterations required to reach convergence is significantly reduced. The algorithms have been implemented in the Atomic Simulation Environment (ASE) open source software.
  • Long-timescale simulations of the diffusion of a H$_2$O admolecule on the (0001) basal plane of ice Ih were carried out over a temperature range of 100 to 200 K using the adaptive kinetic Monte Carlo method and TIP4P/2005f interaction potential function. The arrangement of dangling H atoms was varied from the proton-disordered surface to the perfectly ordered Fletcher surface. A large variety of sites was found leading to a broad distribution in adsorption energy at both types of surfaces. Up to 4 % of the sites on the proton-disordered surface have an adsorption energy exceeding the cohesive energy of ice Ih. The mean squared displacement of a simulated trajectory at 175 K for the proton-disordered surface gave a diffusion constant of 6$\cdot$10$^{-10}$ cm$^2$/s, consistent with an upper bound previously reported from experimental measurements. During the simulation, dangling H atoms were found to rearrange so as to reduce clustering, thereby approaching a linear Fletcher type arrangement. Diffusion on the perfectly ordered Fletcher surface was estimated to be significantly faster, especially in the direction along the rows of dangling hydrogen atoms. From simulations over the range in temperature, an effective activation energy of diffusion was estimated to be 0.16 eV and 0.22 eV for diffusion parallel and perpendicular to the rows, respectively. Even a slight disruption of the rows of the Fletcher surface made the diffusion isotropic.
  • The rate of escape of an ideal bead-spring polymer in a symmetric double-well potential is calculated using transition state theory (TST) and the results compared with direct dynamical simulations. The minimum energy path of the transitions becomes flat and the dynamics diffusive for long polymers making the Kramers-Langer estimate poor. However, TST with dynamical corrections based on short time trajectories started at the transition state gives rate constant estimates that agree within a factor of two with the molecular dynamics simulations over a wide range of bead coupling constants and polymer lengths. The computational effort required by the TST approach does not depend on the escape rate and is much smaller than that required by molecular dynamics simulations.
  • The rate of escape of polymers from a two-dimensionally confining potential well has been evaluated using self-avoiding as well as ideal chain representations of varying length, up to 80 beads. Long timescale Langevin trajectories were calculated using the path integral hyperdynamics method to evaluate the escape rate. A minimum is found in the rate for self-avoiding polymers of intermediate length while the escape rate decreases monotonically with polymer length for ideal polymers. The increase in the rate for long, self-avoiding polymers is ascribed to crowding in the potential well which reduces the free energy escape barrier. An effective potential curve obtained using the centroid as an independent variable was evaluated by thermodynamic averaging and Kramers rate theory then applied to estimate the escape rate. While the qualitative features are well reproduced by this approach, it significantly overestimates the rate, especially for the longer polymers. The reason for this is illustrated by constructing a two-dimensional effective energy surface using the radius of gyration as well as the centroid as controlled variables. This shows that the description of a transition state dividing surface using only the centroid fails to confine the system to the region corresponding to the free energy barrier and this problem becomes more pronounced the longer the polymer is. A proper definition of a transition state for polymer escape needs to take into account the shape as well as the location of the polymer.
  • A method for finding minimum energy paths of transitions in magnetic systems is presented and used to determine the mechanism and estimate the activation energy of skyrmion and antivortex annihiliation in nano-systems. The path is optimized with respect to orientation of the magnetic vectors while their magnitudes are fixed or obtained from separate calculations. The curvature of the configuration space is taken into account by: (1) using geodesics to evaluate distances and displacements of the system during the optimization, and (2) projecting the path tangent and the magnetic force on the tangent space of the manifold defined by all possible orientations of the magnetic vectors. The method, named geodesic nudged elastic band (GNEB), and its implementation are illustrated with calculations of complex transitions involving annihilation and creation of skyrmion and antivortex states. The lifetime of the latter was determined within harmonic transition state theory using a noncollinear extension of the Alexander-Anderson model.
  • We report results of long timescale adaptive kinetic Monte Carlo simulations aimed at identifying possible molecular reordering processes on both proton-disordered and ordered (Fletcher) basal plane (0001) surfaces of hexagonal ice. The simulations are based on a force field for flexible molecules and span a time interval of up to 50 {\mu}s at a temperature of 100 K, which represents a lower bound to the temperature range of Earth's atmosphere. Additional calculations using both density functional theory and an ab initio based polarizable potential function are performed to test and refine the force field predictions. Several distinct processes are found to occur readily even at this low temperature, including concerted reorientation (flipping) of neighboring surface molecules, which changes the pattern of dangling H-atoms, and the formation of interstitial defects by the downwards motion of upper-bilayer molecules. On the proton-disordered surface, one major surface roughening process is observed that significantly disrupts the crystalline structure. Despite much longer simulation time, such roughening processes are not observed on the highly ordered Fletcher surface which is energetically more stable because of smaller repulsive interaction between neighboring dangling H-atoms. However, a more localized process takes place on the Fletcher surface involving a surface molecule transiently leaving its lattice site. The flipping process provides a facile pathway of increasing proton-order and stabilizing the surface, supporting a predominantly Fletcher-like ordering of low-temperature ice surfaces, but our simulations also show that proton- disordered patches on the surface may induce significant local reconstructions. Further, a subset of the molecules on the Fletcher surface are susceptible to forming interstitial defects.
  • Calculations of stable and metastable magnetic states as well as minimum energy paths for transitions between states are carried out using a noncollinear extension of the multiple-impurity Alexander-Anderson model and a magnetic force theorem which is derived and used to evaluate the total energy gradient with respect to orientation of magnetic moments -- an important tool for efficient navigation on the energy surface. By using this force theorem, the search for stable and metastable magnetic states as well as minimum energy paths revealing the mechanism and activation energy of transitions can be carried out efficiently. For Fe monolayer on W(110) surface, the model gives magnetic moment as well as exchange coupling between nearest and next-nearest neighbors that are in good agreement with previous density functional theory calculations. When applied to nanoscale Fe islands on this surface, the magnetic moment is predicted to be 10\% larger for atoms at the island rim, explaining in part an experimentally observed trend in the energy barrier for magnetization reversal in small islands. Surprisingly, the magnetic moment of the atoms does not change much along the minimum energy path for the transitions, which for islands containing more than 15 atom rows along either $[001]$ or $[1\bar{1}0]$ directions involves the formation of a thin, temporary domain wall. A noncollinear magnetic state is identified in a $7\times 7$ atomic row Fe island where the magnetic moments are arranged in an antivortex configuration with the central ones pointing out of the $(110)$ plane. This illustrates how the model can describe complicated exchange interactions even though it contains only a few parameters. The minimum energy path between this antivortex state and the collinear ground state is also calculated and the thermal stability of the antivortex state estimated.
  • A method is presented for generating a good initial guess of a transition path between given initial and final states of a system without evaluation of the energy. An objective function surface is constructed using an interpolation of pairwise distances at each discretization point along the path and the nudged elastic band method then used to find an optimal path on this image dependent pair potential (IDPP) surface. This provides an initial path for the more computationally intensive calculations of the true minimum energy path using some method of choice for evaluating the energy and atomic forces, for example by ab initio or density functional theory. The optimal path on the IDPP surface is significantly closer to the true minimum energy path than a linear interpolation of the Cartesian coordinates and, therefore, reduces the number of iterations needed to reach convergence and averts divergence in the electronic structure calculations when atoms are brought too close to each other in the initial path. The method is illustrated with three examples: (1) rotation of a methyl group in an ethane molecule, (2) an exchange of atoms in an island on a crystal surface, and (3) an exchange of two Si-atoms in amorphous silicon. In all three cases, the computational effort in finding the minimum energy path with DFT was reduced by a factor ranging from 50 % to an order of magnitude by using an IDPP path as the initial path. The time required for parallel computations was reduced even more because of load imbalance when linear interpolation of Cartesian coordinates was used.
  • Self-consistent calculations using the Perdew-Zunger self-interaction correction (PZ-SIC) to local density and gradient dependent energy functionals are presented for the binding energy and equilibrium geometry of small molecules as well as energy barriers of reactions. The effect of the correction is to reduce binding energy and bond lengths and increase activation energy barriers when bond breaking is involved. The accuracy of the corrected functionals varies strongly, the correction to the binding energy being too weak for the local density approximation but too strong for the gradient dependent functionals considered. For the Perdew, Burke, and Ernzerhof (PBE) functional, a scaling of the PZ-SIC by one half gives improved results on average for both binding energy and bond lengths. The PZ-SIC does not necessarily give more accurate total energy, but it can result in a better cancellation of errors. An essential aspect of these calculations is the use of complex orbitals. A restriction to real orbitals leads to less accurate results as was recently shown for atoms [S. Kl\"upfel, P. Kl\"upfel, and H. J\'onsson, Phys. Rev. A {84}, 050501 (2011)]. The molecular geometry of radicals can be strongly affected by PZ-SIC. An incorrect, non-linear structure of the \ce{C2H} radical predicted by PBE is corrected by PZ-SIC. The \ce{CH3} radical is correctly predicted to be planar when complex orbitals are used, while it is non-planar when the PZ-SIC calculation is restricted to real orbitals.
  • The ground state of atoms from H to Ar was calculated using a self-interaction correction to local and gradient dependent density functionals. The correction can significantly improve the total energy and makes the orbital energies consistent with ionization energies. However, when the calculation is restricted to real orbitals, application of the self-interaction correction can give significantly higher total energy and worse results, as illustrated by the case of the Perdew-Burke-Ernzerhof gradient dependent functional. This illustrates the importance of using complex orbitals for systems described by orbital density dependent energy functionals.
  • The effect of hydrogen adsorption on the magnetic properties of an Fe$_3$ cluster immersed in a Cu(111) surface has been calculated using densifty functional theory and the results used to parametrize an Alexander-Anderson model which takes into account the interaction of d-electrons with itinerant electrons. A number of adatom configurations containing one to seven H-atoms were analyzed. The sequential addition of hydrogen atoms is found to monotonically reduce the total magnetic moment of the cluster with the effect being strongest when the H-atoms sit at low coordinated sites. Decomposition of the charge density indicates a transfer of 0.4 electrons to each of the H-atoms from both the Fe-atoms and from the copper substrate, irrespective of adsorption site and coverage. The magnetic moment of only the nearest neighbor Fe-atoms is reduced and mainly due to increased population of minority spin d-states. This can be modeled by increased indirect coupling of d-states via the conduction s-band in the Alexander-Anderson model.
  • Theoretical calculations of the structure, formation and migration of kinks on a non-dissociated screw dislocation in silicon have been carried out using density functional theory calculations as well as calculations based on interatomic potential functions. The results show that the structure of a single kink is characterized by a narrow core and highly stretched bonds between some of the atoms. The formation energy of a single kink ranges from 0.9 to 1.36 eV, and is of the same order as that for kinks on partial dislocations. However, the kinks migrate almost freely along the line of an undissociated dislocation unlike what is found for partial dislocations. The effect of stress has also been investigated in order to compare with previous silicon deformation experiments which have been carried out at low temperature and high stress. The energy barrier associated with the formation of a stable kink pair becomes as low as 0.65 eV for an applied stress on the order of 1 GPa, indicating that displacements of screw dislocations likely occur via thermally activated formation of kink pairs at room temperature.
  • A practical method for finding free energy barriers for transitions in high-dimensional classical and quantum systems is presented and used to calculate the dissociative sticking probability of H2 on a metal surface within transition state theory (TST). The reversible work involved in shifting the system confined to a hyperplane from the reactant region towards products is evaluated directly. Quantum mechanical degrees of freedom are included by using Feynman Path Integrals with the hyperplane constraint applied to the centroid of the cyclic paths. An optimal dividing surface for the rate estimated by TST is identified naturally in the course of the reversible work evaluation. The free energy barrier is determined relative to the reactant state directly, so an estimate of the transition rate can be obtained without requiring a solvable reference model for the transition state. The method has been applied to calculations of the sticking probability of a thermalized hydrogen gas on a Cu(110) surface. The two hydrogen atoms were included quantum mechanically, and over two hundred atoms in the Cu crystal where included classically. The activation energy for adsorption and desorption was determined and found to be significantly lowered by tunneling at low temperature. The calculated values agree quite well with experimental estimates. Dynamical corrections to the classical TST rate estimate were evaluated and found to be small.
  • We have developed a flexible hybrid decomposition parallel implementation of the first-principles molecular dynamics algorithm of Car and Parrinello. The code allows the problem to be decomposed either spatially, over the electronic orbitals, or any combination of the two. Performance statistics for 32, 64, 128 and 512 Si atom runs on the Touchstone Delta and Intel Paragon parallel supercomputers and comparison with the performance of an optimized code running the smaller systems on the Cray Y-MP and C90 are presented.