
We present an efficient implementation of a surface Green'sfunction 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 semiinfinite 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 thinfilm
semiconductor heterostructures, surface states in metals and topological
insulators, and surfaces in external electrical fields. Results obtained with
the surface Green'sfunction 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 lesserknown 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
twodimensional M\"ullerBrown 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 Heisenberglike Hamiltonian that includes DzyaloshinskiiMoriya 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
preexponential factor in the Arrhenius law for the rate, the socalled 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,
socalled '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 nonmagnetic 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 nonmagnetic 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. Nonmagnetic
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 spinpolarized current preferably occurring near
structural defects.

The calculation of minimum energy paths for transitions such as atomic and/or
spin rearrangements 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
PipekMezey 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 PipekMezey 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 socalled 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 crosscomparison 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 realspace 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 GeSi 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/SmCo 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 2dimensional 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 gasphase 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
firstorder 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.

Longtimescale 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 protondisordered 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
protondisordered 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 protondisordered 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 beadspring polymer in a symmetric doublewell
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 KramersLanger 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 twodimensionally confining potential
well has been evaluated using selfavoiding 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 selfavoiding
polymers of intermediate length while the escape rate decreases monotonically
with polymer length for ideal polymers. The increase in the rate for long,
selfavoiding 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 twodimensional 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 nanosystems. 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 AlexanderAnderson model.

We report results of long timescale adaptive kinetic Monte Carlo simulations
aimed at identifying possible molecular reordering processes on both
protondisordered 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 Hatoms, and the formation of interstitial defects by the downwards
motion of upperbilayer molecules. On the protondisordered 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 Hatoms. 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
protonorder and stabilizing the surface, supporting a predominantly
Fletcherlike ordering of lowtemperature 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 multipleimpurity AlexanderAnderson 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 nextnearest
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 Siatoms 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.

Selfconsistent calculations using the PerdewZunger selfinteraction
correction (PZSIC) 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 PZSIC by one half gives improved results on
average for both binding energy and bond lengths. The PZSIC 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 PZSIC. An incorrect, nonlinear structure of the \ce{C2H}
radical predicted by PBE is corrected by PZSIC. The \ce{CH3} radical is
correctly predicted to be planar when complex orbitals are used, while it is
nonplanar when the PZSIC calculation is restricted to real orbitals.

The ground state of atoms from H to Ar was calculated using a
selfinteraction 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
selfinteraction correction can give significantly higher total energy and
worse results, as illustrated by the case of the PerdewBurkeErnzerhof
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 AlexanderAnderson
model which takes into account the interaction of delectrons with itinerant
electrons. A number of adatom configurations containing one to seven Hatoms
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 Hatoms sit at low coordinated sites. Decomposition of
the charge density indicates a transfer of 0.4 electrons to each of the Hatoms
from both the Featoms and from the copper substrate, irrespective of
adsorption site and coverage. The magnetic moment of only the nearest neighbor
Featoms is reduced and mainly due to increased population of minority spin
dstates. This can be modeled by increased indirect coupling of dstates via
the conduction sband in the AlexanderAnderson model.

Theoretical calculations of the structure, formation and migration of kinks
on a nondissociated 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
highdimensional 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 firstprinciples 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 YMP and C90 are presented.