Nonadiabatic potential-energy surfaces by constrained density-functional theory

Nonadiabatic effects play an important role in many chemical processes. In order to study the underlying nonadiabatic potential-energy surfaces PESs , we present a locally constrained density-functional theory approach, which enables us to confine electrons to subspaces of the Hilbert space, e.g., to selected atoms or groups of atoms. This allows one to calculate nonadiabatic PESs for defined charge and spin states of the chosen subsystems. The capability of the method is demonstrated by calculating nonadiabatic PESs for the scattering of a sodium and a chlorine atom, for the interaction of a chlorine molecule with a small metal cluster, and for the dissociation of an oxygen molecule at the Al 111 surface.


I. INTRODUCTION
Assuming that electrons can react much faster to external perturbations than nuclei, the Born-Oppenheimer approximation 1 ͑BOA͒ is a widely employed approach to separate electronic and nuclear motion in dynamic processes. The nuclei then appear static to the electrons, which in turn set up a potential governing the motion of the nuclei. Formally the approach starts with the general many-body Schrödinger equation where T nuc is the kinetic energy operator of the nuclei, V nuc-nuc the interaction potential of the nuclei, and H e the electronic Hamiltonian containing the electron kinetic energy, as well as the electron-electron and electron-nuclei interactions. In an adiabatic representation, 2 where the ⌽ adia are chosen to be the eigenfunctions of the electronic Hamiltonian at the actual position of the nuclei Inserting Eq. ͑2͒ into Eq. ͑1͒ leads to a set of equations for the wave functions of the nuclei E⌳ adia = ͑T nuc + V nuc-nuc + E e ͒⌳ adia + nonadiabaticity terms,

͑4͒
in which the "nonadiabaticity terms" summarize matrix elements of the momentum and kinetic-energy operators of the nuclear motion. The Born-Oppenheimer approximation corresponds to setting these terms to zero, which implies that the electrons assume their electronic state instantaneously for any position of the nuclei, unaffected by the nuclear dynamics. The nuclei are then moving on the potential-energy surface ͑PES͒ V adia = V nuc-nuc + E e , also called the Born-Oppenheimer surface of the electronic state .
If, in the BOA, the system is initially in the electronic ground state, it will remain there irrespective of the dynamics of the nuclei. However, in real life electrons may not be able to follow the motion of the nuclei instantaneously, and, e.g., when selection rules apply, they may find themselves in an excited state. For example, chemical reactions forming singlet molecules from triplet and singlet reactants are forbidden by Wigner's spin-selection rule. The triplet multiplicity is the actual reason why most reactions of O 2 with other molecules or substances, although being exothermic, do not proceed at room temperature; they are kinetically hindered. In other words, in a chemical reaction the spin of the reactants must be conserved or transferred to some other entity. The transition from the O 2 ground state ͑ 3 ⌺ g − ͒ to the first excited state ͑ 1 ⌬ g ͒ is strictly forbidden for the isolated molecule, as is the reverse ͑deexcitation͒ process, once the molecule has been excited to the 1 ⌬ g state. Thus, when probabilities for transitions between different electronic states are low, e.g., due to selection rules, the assumption that the system will always remain in the electronic ground state may become incorrect. For the mentioned O 2 example this implies that when an external field shifts the 1 ⌬ g energy below the 3 ⌺ g − energy, the probability for an electronic transition toward the 1 ⌬ g state will be low. Indeed, this is an important aspect of the O 2 /Al͑111͒ interaction, one of our examples discussed below.
For such, or alike situations it is necessary to go beyond the BOA and to consider the coupling between different states. 5,6 For a full description the nonadiabaticity terms must be calculated, which in practice means to evaluate the matrix elements for the derivative couplings of the nuclear and electronic motion. For this, the use of another set of electronic functions ͉⌽ dia ͘ in the wave function expansion in Eq. ͑2͒ may be suitable. The idea behind such "diabatic" states is that in a dynamic ͑scattering͒ event, the electrons lag behind the nuclear motion and thus tend to conserve a given ͑initial͒ character of the electronic structure. Formally, the ͉⌽ dia ͘ are then constructed to minimize the nuclear derivative coupling terms, such that the matrix representation of the nuclear momentum operator becomes diagonal in the diabatic basis. [6][7][8] The equivalent to Eq. ͑4͒ in a diabatic representation is thus a matrix equation where the diagonal elements of the potential V dia = V nuc-nuc + H e are the diabatic PESs, and the off-diagonal elements H e describe the coupling between PES and PES . For the above example of an O 2 molecule, a suitable diabatic state would, e.g., be given by the 3 ⌺ g − state of the isolated molecule, and the corresponding diabatic PES can be used to describe the motion of a molecule that tries to conserve this electronic configuration, even in regions where external fields ͑or interactions with other species͒ bring the 1 ⌬ g below the 3 ⌺ g − energy. The aim of this work is to develop and apply a nonadiabatic approach, with a similar physical motivation as the diabaticity concept. We account for the tendency to maintain a given character of the electronic structure, e.g., due to selection rules, by focusing on the dynamics of a system with defined charge and spin in suitably chosen subsystems, typically the two reactants in a scattering event. Using again the example of an O 2 molecule, this could, e.g., be a state with fixed triplet spin on the molecule, even when it interacts with another reactant. Instead of using the properties of the manybody wave function as the formal basis, we thus build our approach on a defined total spin ͑or total charge, or other well-defined quantity͒ in the local Hilbert space of a chosen subsystem, which can be an atom or a group of atoms. The thus defined nonadiabatic PESs are computed with densityfunctional theory ͑DFT͒, 9,10 and, in particular, using the idea of constrained DFT formulated by Dederichs et al. in 1984. 11 The basic idea is to modify the energy functional employed in DFT by applying physically meaningful constraints. The constraint is enforced through an additional Lagrange multiplier, which on the level of the Kohn-Sham ͑KS͒ equations yields an additional term to the KS effective potential. In the present work we employ a local constraint via a projection technique that enables us to freeze the charge and spin of the chosen subsystems.
Having stated this general philosophy, two points deserve a critical discussion. First, we like to point out that in practice our nonadiabatic states may often be close to those of the diabaticity concept. Still, the formal definition is different, and in some cases, such as the example of the O 2 /Al͑111͒ system discussed below, there are notable differences. In contrast to recent other works in this area, 12,13 we therefore prefer to refer to our states simply as nonadiabatic states to underscore this difference to the established diabaticity concept. As a second point, we note that a solid, mathematical proof of the validity of constrained DFT does not exist. Still, spin-density functional theory, the Slater-Janak transition state, 11,14,15 the fixed-spin moment ͑FSM͒ approach, 16 and some other examples present frequently employed, important, and successful applications of essentially the same concept. Our approach represents a mild generalization of such applications, and our constraint is plausible and physically meaningful. The spins or charges of two interacting systems may be hindered to adjust or combine when there are selection rules, or when the interaction time is very short. Thus, the appropriate total-energy surfaces of the scattering event should ͑up to a certain distance͒ be that of conserved local spins or charges.
In a previous publication we have already used the present approach to resolve a long-standing problem in surface science, namely, the low-sticking probability of oxygen molecules at the Al͑111͒ surface. 17 However, the method is much more general, as we will illustrate below by also applying it to two other, nonperiodic model systems, namely, a scattering event of a sodium and a chlorine atom, and the interaction of a Cl 2 molecule with a magnesium cluster. In this respect, we also mention notable, independent work of Wu and Van Voorhis, 13,18 which is essentially equivalent to our approach 17 in that also in their work an additional potential is introduced to constrain electron numbers in welldefined parts of a system. As in our approach this potential is determined in a self-consistent way using formally identical methodology, and has been employed to study chargetransfer reactions.

II. LOCALLY CONSTRAINED DENSITY-FUNCTIONAL THEORY
Our locally constrained DFT ͑LC-DFT͒ approach starts by assigning electrons to defined subspaces of the total Hilbert space, e.g., to atoms or groups of atoms. This is done by employing a suitable projection scheme to distinguish the individual subsystems. In the following sections we will present this formalism for a system consisting of two subsystems called A and B, with straightforward generalization to more than two subsystems. Considering electrons and their spin, this leads to four electron numbers and N B ↓ , which uniquely define the PES V N A for the nonadiabatic quantum state with fixed spins and charges of the subsystems. Expressing the constraints in terms of auxiliary potentials, the electronic structure problem is solved self-consistently using DFT, i.e., the electronic structure is fully relaxed under the given constraints yielding

A. Definition of the subsystems
In order to assign electrons to the two subsystems, the Kohn-Sham single-particle wave functions are expanded into localized, atom-centered basis functions, e.g., Gaussians, Slater-type orbitals, or numerical orbitals. The implementation is therefore particularly convenient in codes employing such basis sets, whereas in codes based on other basis sets, such as plane waves, the implementation requires intermediate projection steps onto localized functions. [19][20][21] The Kohn-Sham orbital i with index i and spin index = ↑ , ↓ is thus written as a linear combination of all basis functions j , or in calculations using periodic boundary conditions as a linear combination of k-dependent Bloch basis functions j k = e ikr j , In the following our derivation will refer to the periodic case, while for finite systems the dependence on the k point index would be simply dropped.
In the case of atom-centered basis functions each basis function is uniquely assigned either to subsystem A or to subsystem B: All basis functions centered at atoms being part of subsystem A define the Hilbert space of subsystem A, and all basis functions centered at atoms being part of subsystem B define the Hilbert space of subsystem B. Every singleparticle wave function can then be projected onto the two Hilbert subspaces, which is done separately for each k point and each spin, and taking into account a nonorthogonality of the atomic orbitals by including the overlap matrix elements Here, p A,i k is the projection onto subsystem A, p B,i k the projection onto subsystem B, and n is the total number of basis functions, of which the first k =1, ... ,m are the basis functions of subsystem A. i,A k is defined as the A component of i k , i.e., all coefficients referring to basis functions of subsystem B are set to zero. Correspondingly, all coefficients referring to basis functions of subsystem A are set to zero in the functions i,B k with i = m +1, ... ,n, which define the complementary B component of i k . We note that the nor- This relation can be used to reduce the computational effort in that just one of the two double sums in Eq. ͑7͒ needs to be calculated and the other is obtained from the normalization condition. We note that in this scheme like in all projection schemes the actual values of p A,i k and p B,i k depend on the choice of the basis sets defining the Hilbert spaces of the subsystems.
Having split each single-particle state into an A part and a B part allows one to construct the partial densities of states ͑pDOSs͒ for subsystems A and B and for the two spin channels. Summing up the resulting four pDOSs over all occupied single-particle states i yields then the four electron numbers where f i k is the occupation number of the single-particle Kohn-Sham state i, typically chosen to be a Fermi function. With this, also the total spin S of each of the two subsystems is determined through the difference of the corresponding electron numbers

B. Constraining the electron numbers
Having established the various electron numbers, we proceed to introduce the constraint. Aiming to compute the nonadiabatic PES V N A ↓ , as well as four partial electron densities, which sum to the total electron density. If the Fermi energies were all degenerate at this stage, the chosen nonadiabatic charge and spin state of the system would correspond to the adiabatic ground state. However, typically the four Fermi energies are all different, reflecting the fact that there is no self-consistency between the total electron density and the effective Kohn-Sham potential.
In order to achieve this self-consistency, we choose to align the Fermi energies, still requiring that the electron numbers that can be filled into the four pDOSs up to the resulting common Fermi level remain unchanged. We therefore employ a method that shifts each of the four pDOSs independently. This is particularly important when hybridization between the two subsystems is present and a singleparticle state i contains nonzero expansion coefficients c ij k for basis functions belonging to subsystem A, as well as for basis functions belonging to subsystem B. In such cases, simply shifting the entire state i as is, e.g., done in the ⌬SCF method would not change the relative positions of the A and B Fermi energies. Instead, it is necessary to act separately on the basis functions in both subsystems. Only this gives the electronic structure enough flexibility to fully relax under the imposed constrained electron numbers, which in the end will yield a different expansion of state i in terms of A and B basis functions.
In practice, we first align the Fermi energies separately for each spin. Specifically, we request that ⑀ F,A = ⑀ F,B = ⑀ F , and appropriately shift the A and B Fermi energies to the Fermi energy ⑀ F . This is achieved by adding to the Kohn-Sham Hamiltonian auxiliary potentials that act differently on the A and B basis functions. These auxiliary potentials consist of a "strength" factor ⌫ and a projection operator P k onto the subspaces where the summation is done over the m basis functions spanning the Hilbert space of subsystem A and the n-m basis functions spanning the Hilbert space of subsystem B. In the chosen form P k symmetrically contains covariant and contravariant basis functions, 22,23 i k and i,k , respectively. The resulting matrix representations of P A k and P B k in the Hilbert space spanned by the Bloch basis functions is then Hermitian, which facilitates the implementation into existing DFT codes as further described below. As derived in the Appendix for subsystem A, the form of these matrices is as shown schematically in Fig. 1. The matrix element P A,ij k is equal to the overlap matrix element S ij k , if i and j both refer to basis functions assigned to subsystem A. If only i or j refer to a basis function assigned to subsystem A, the matrix element is Finally, if neither i nor j belong to subsystem A, P A,ij k is zero, reflecting that the pDOS of subsystem B is not affected by the auxiliary potential ⌫ A P A k . The auxiliary potential ⌫ B P B k has an analogous structure.
Since the purpose of the auxiliary potentials is to align the Fermi energies of subsystems A and B with ⑀ F , the obvious choice for the "strength" factors ⌫ A and ⌫ B is However, because of the resulting nonzero auxiliary potentials, the initial single-particle states are no longer solutions of the new effective Hamiltonian for each spin, comprised of both the Kohn-Sham Hamiltonian and the auxiliary potential. As a result, ⌫ A and ⌫ B must be determined selfconsistently. Diagonalization of the new effective Hamiltonians yields new eigenvectors to construct new partial densities of states, the Fermi levels of which define new strength factors through Eq. ͑11͒. This is repeated in a selfconsistency ͑SC͒ cycle, until the Fermi energies of subsystems A and B are aligned to an arbitrary precision ͑for each spin channel͒, The ensuing step of aligning the two different spin Fermi energies ⑀ F ↑ and ⑀ F ↓ is done in an analogous way, i.e., by adding another auxiliary potential of the form of Eq. ͑10͒. In this case, the matrix structure of the corresponding projection operator P k onto one spin subspace is simpler though. Since this subspace is spanned by all n basis functions, regardless of whether they are in subsystem A or B, the sum in Eq. ͑10͒ goes up to n, and the matrix representation of P k becomes simply the overlap matrix, cf. Fig. 1 with m = n. Adding an auxiliary potential ⌫ P k to the effective Hamiltonian resulting from the preceding alignment of the A and B Fermi energies, corresponds therefore to a mere shift of the eigenvalues, depending on the chosen strength factor ⌫ . Here we choose to shift the spin-up and spin-down pDOSs in opposite As before, this procedure has to be done in a self-consistent way, since adding the new auxiliary potential modifies the effective Hamiltonian. In fact, since the alignment of the A and B Fermi levels and the alignment of the spin-up and spin-down Fermi levels is not independent of each other, the two SC cycles must be nested. Discussing below how this can be implemented in a numerically efficient way into existing DFT codes, we note that once the double selfconsistency is achieved, we arrive at a final common Fermi level in the system and a new set of single-particle states i ͑with eigenvalues ⑀ i Ј and eigenvectors i Ј ͒ that is the selfconsistent solution to the effective Hamiltonian, containing the original Kohn-Sham Hamiltonian and the two auxiliary potentials. In each subspace spanned by the Bloch basis functions at one k point we therefore have , and ⌫ SCF as determined in the last cycle of the nested SC loops.

C. Numerical considerations
At this stage it is appropriate to recall what has been achieved so far. The imposed constraint of fixed electron numbers N A ↑ , N A ↓ , N B ↑ , and N B ↓ has been suitably transformed into external potentials. Adding these to the Kohn-Sham Hamiltonian led to the effective Hamiltonian of Eq. ͑13͒. As proven by the Hohenberg-Kohn theorem, the electron density resulting from occupying all single-particle states up to the common Fermi level corresponds therefore to the fully relaxed electronic structure under the given constraint. Calculating the total energy E N A connected to this electron density provides then ͑together with V nuc-nuc ͒ the nonadiabatic PES V N A , representing the chosen spin and electron numbers in the two subsystems.
As much as the total energy, also any other quantity that is a function of the electron density can be computed. However, care has to be taken, if the actual calculation involves terms that explicitly contain the single-particle eigenvalues. This becomes apparent when writing ⑀ i Ј k with Eq. ͑13͒ as which shows that ⑀ i Ј k contains actually two contributions, of which only the first, ⑀ i,KS Ј k , really reflects the modified electronic structure. The second term, ⌬⑀ i,pot Ј k , is spurious and results from the aligning potentials. When evaluating terms containing single-particle eigenvalues, one should therefore either directly use ⑀ i,KS Ј k or correct the spurious contribution through correspondingly computed terms containing ⌬⑀ i,pot Ј k .
This mainly applies to the Pulay force correction term, which depends explicitly on the single-particle eigenvalues. 25,26 Finally, we mention some numerical considerations. At first glance, it appears as if our LC-DFT formalism builds on the self-consistent solutions to the Kohn-Sham Hamiltonian and then requires two additional nested SC loops. However, for the latter self-consistency under the imposed constraint there is no need to start with self-consistent Kohn-Sham solutions. Additionally, the mere eigenvalue shift induced by the spin Fermi-level alignment step is nicely compatible with the structure of existing DFT codes. From a computational point of view, our algorithm can thus be implemented as one additional SC cycle at each iteration of the existing electronic SC cycle in a DFT code. Of course, all known approaches to improve the convergence of SC cycles such as the use of sophisticated mixing schemes can equally be applied to the additional cycle. Having implemented a Pulay mixing scheme, 24 our experience was that for the tested systems only the very first DFT iterations required a significant number of inner iterations, typically about 5-10, while close to the outer self-consistency also the inner self-consistency was often directly reached after the first iteration. For the example involving the Al͑111͒ surface, we even found the overall DFT convergence, i.e., the number of iterations in the outer SC cycle, frequently much improved compared to standard adiabatic calculations, because the oscillation of states around the Fermi level is reduced when controlling the electron numbers and thereby also the occupation numbers of the Kohn-Sham states close to the Fermi level. For the systems with localized electronic states, we found that using the spins or electron numbers in the subsystems as convergence criterion for the self-consistency was sometimes preferred to the equality of the Fermi energies. Overall, for the systems studied, the cost of the calculations using the constraint was thus often about the same and sometimes even lower than the cost of standard adiabatic calculations.

D. Comparison to ⌬SCF and fixed-spin moment
An important characteristic of the LC-DFT approach presented here is that it is parameter-free and only the set of the electron numbers in the four channels defining the nonadiabatic state of interest has to be specified. Under this constraint, the electronic structure is fully relaxed, i.e., the partial densities of states are not frozen and shifted statically. Instead, the single-particle states are flexible to vary the contribution of each basis function to each single-particle state freely. This can lead to a significant improvement compared to the two prominent and widely employed implementations of the constrained DFT concept, the ⌬SCF ͑Ref. 14͒ and the fixed-spin moment approach. 16 In the latter approach the system is only separated into a spin-up and a spin-down channel, which are filled independently. In LC-DFT we go a step further and allow one also to distribute the spin-up and spindown electrons in a well-defined way into the two subsystems, which permits a more general control over the spatial distribution of the electron and magnetization densities. The different results obtained with both methods are particularly obvious for the oxygen dissociation at the Al͑111͒ case described in Sec. III C below.
The improvement compared to the ⌬SCF method concerns extended systems. Both approaches have in common that they consider two subsystems and first analyze the pDOSs to determine the contributions of each subsystem to the individual single-particle states. However, in the ⌬SCF method the states with high A parts are then completely assigned to subsystem A, regardless of their B contributions, while the states with small A contributions are completely assigned to subsystem B. 28 A subsequent occupation of the states by 0 or 1 electrons is therefore only fully justified for the typically not interesting case of noninteracting subsystems, i.e., no hybridization. Otherwise it results in fractional effective occupation numbers of the individual subsystems, which necessarily introduces some uncertainty in the total energies obtained in the ⌬SCF method. The LC-DFT approach, on the other hand, allows for a physical rehybridization of the states under the imposed constraint. It thus conserves the electron numbers of both subsystems, while at the same time fully taking hybridization into account. This is the main difference between the present LC-DFT formalism and the ⌬SCF method, while in the limit of infinitely separated subsystems both approaches are equivalent.
It is finally also important to note that all methods discussed here, the FSM approach, the ⌬SCF method, and the LC-DFT formalism, intend to overcome limitations in the description of chemical processes by controlling the electron numbers. This obviously does not allow one to overcome approximations in the employed exchange-correlation ͑xc͒ functional. Local-density or gradient-corrected xc functionals are, e.g., known to cause inaccuracies in the energy splittings between different spin multiplets. [29][30][31] To this end, we note that in the limit of infinite separation between the subsystems ͑reactants͒, our approach reduces to the description of two isolated subsystems. The nonadiabatic states of defined spin or charge in our approach are then simply the corresponding excited states of the isolated subsystems, i.e., of the noninteracting reactants. It is worthwhile to point out that this is different to the diabaticity concept, which reduces in this limit of infinite separation to the adiabatic ͑excited͒ states of the interacting system. If there is a ͑unphysical͒ charge or spin transfer even at these distances ͓as in the example of O 2 /Al͑111͒ discussed below͔, it will be present in both the adiabatic and the diabatic description, but not in our nonadiabatic approach.

III. APPLICATIONS
As an important application of the LC-DFT method we now illustrate its use in the calculation of nonadiabatic PESs that are of interest in the investigation of dynamic processes. Specifically, we focus here on the scattering of atoms or molecules in crossed molecular beams or at solid surfaces. As a side effect this also shows how the method can be employed to suitably restrict unwanted electron transfer between weakly or noninteracting subsystems. In a DFT calculation, such an electron transfer will, e.g., occur whenever an occupied level of subsystem A is higher in energy than an unoccupied level of subsystem B. Alignment of the Fermi levels of the two systems will then lead to a fractional occupation of both states in the self-consistent solution, even if the distance between the two subsystems is macroscopic. A prominent example for this is the interaction of an oxygen molecule with the Al͑111͒ surface, where the unoccupied 2 *↓ orbitals of the molecule are lower than the Fermi level of the metal. 17,27 The resulting electron transfer lowers the total energy of the system and at macroscopic distances the latter does then not converge to the physically correct limit, given by the sum of the total energies of the isolated molecule and the isolated surface.
As an example for an extended periodic system, we correspondingly briefly discuss the interaction of O 2 with the Al͑111͒ surface. In addition, we also present calculated nonadiabatic PESs for two finite model systems, namely, for the scattering of Na and Cl, and the scattering of Cl 2 at a small Mg 4 cluster. All calculations have been carried out using the all-electron DFT code DMol 3 , which employs numerical atomiclike orbitals as basis functions. 32 Unless otherwise noted, we employ the DMol 3 standard basis set "all" consisting of atomic orbitals and polarization functions, 32 a realspace cutoff of 12 bohrs for the basis functions, and the PBE ͑Ref. 33͒ xc functional.

A. Na+ Cl
In the scattering of a sodium and a chlorine atom, two nonadiabatic states of interest are the ionic PES "Na + +Cl − " and the neutral PES "Na+ Cl." Since both neutral atoms are spin doublets in their ground states, there are two possible relative orientations of the spins in the latter case, yielding an overall singlet ͑antiparallel spins͒ or triplet ͑parallel spins͒ neutral state. Identifying each atom as one subsystem in our LC-DFT approach, Table I shows the constrained electron and spin numbers we used to represent each of these nonadiabatic states. The resulting PESs as a function of the interatomic separation are shown in Fig. 2, in which we additionally include the computed adiabatic ground state PES and the PES as resulting from a FSM calculation for an overall spintriplet state ͑N ↑ = 15, N ↓ =13͒. In both the latter PES and the neutral triplet PES there are therefore two unpaired electrons, but only the neutral triplet PES computed by LC-DFT has the additional control of locating one unpaired electron explicitly at each atom.
For bond lengths lower than 8 Å, the energetically most favorable nonadiabatic state is found to be the ionic PES, whereas for larger distances these are the degenerate singlet and triplet neutral PESs. The degeneracy of the latter two PESs is only lifted for distances smaller than 5 Å, at which the Pauli repulsion between the two unpaired spin-up electrons leads to a strong increase of the neutral triplet PES. Interestingly, the FSM curve is for all distances virtually degenerate to this neutral triplet PES, indicating that even without constraint one unpaired electron wants to stay at each atom. The minimum of the adiabatic PES, on the other hand, is somewhat lower than the minimum of the ionic PES, showing that the electron transfer in the adiabatic case differs slightly from the one electron imposed in the LC-DFT computation.
Another prominent feature of Fig. 2 is that for large interatomic separations the adiabatic PES is about 1 eV lower than the limit of neutral separated atoms. This is an illustration of the above described electron transfer problem in adiabatic calculations. Even for infinite interatomic distances, a small amount of charge is transferred from the sodium to the chlorine atom, since the lowest unoccupied 3p ↓ state ͑LUMO͒ of the latter, is lower in energy than the highest occupied 3s ↑ state ͑HOMO͒ of the sodium atom. In the selfconsistent calculation, electron density is consequently transferred, until the Fermi levels of the two atoms are aligned. At infinite separation between the two atoms, this charge transfer can be determined quantitatively by calculating the Na Kohn-Sham HOMO and Cl Kohn-Sham LUMO-level energies as a function of different occupations. The results obtained from calculations of the isolated charged atoms are displayed in Fig. 3 and show that HOMO and LUMO are only aligned after an electron transfer of 0.37e. This unphysical electron transfer at infinite separations is not possible in the LC-DFT approach by construction, explaining why in contrast to the adiabatic PESs the neutral triplet and singlet PESs approach the correct limit.  2. ͑Color online͒ Nonadiabatic PESs for the scattering of a Na and a Cl atom. Shown are the energies as a function of the interatomic distance Z. See Table I and the text for the constrained electron numbers defining the various nonadiabatic PESs. Additionally shown are the adiabatic ground state PES and the PES obtained with a FSM triplet calculation. The energy zero corresponds to the energy of the two isolated, neutral atoms.
Finally, we also employed this system to test the proper evaluation of the forces in the constrained LC-DFT calculations. Taking the example of the ionic PES, Fig. 4 shows the force on the Cl atom computed either analytically within the LC-DFT approach or numerically by differentiating the PES shown in Fig. 2. The agreement is excellent proving that forces can be computed accurately within the LC-DFT approach.

B. Mg 4 +Cl 2
As a simple example for subsystems consisting of groups of atoms we discuss the interaction of a Cl 2 molecule with a small metal cluster formed of a tetrahedron of four magnesium atoms. To illustrate the method we restrict ourselves here to computing the PESs as a function of the distance of a Cl 2 molecule approaching the Mg 3 plane of the cluster as explained in Fig. 5. For this we first relaxed the structure of the cluster and the Cl 2 molecule separately, and then held the resulting Mg-Mg bond lengths of 3.07 Å and the Cl-Cl bond length of 2.03 Å fixed in the subsequent calculations. Using the electron numbers compiled in Table I, we calculated the nonadiabatic PESs corresponding to the neutral, the ionicsinglet, and the ionic-triplet state. Together with the adiabatic PES, the resulting curves are shown in Fig. 5. At all distances, the nonadiabatic state corresponding to a neutral configuration exhibits the lowest energy, with the two ionic curves exhibiting significantly higher energies. Similar to the Na+ Cl case the latter two become degenerate at larger distances, when the Pauli repulsion affecting the ionic triplet curve becomes negligible. The closeness of the nonadiabatic neutral curve to the adiabatic result indicates only a comparably small electron transfer from the cluster to the molecule during the scattering process. In this case, there is therefore also only a small electron transfer problem and the adiabatic curve approaches the proper limit for large molecule-cluster separations. By differently occupying the Kohn-Sham HOMO and LUMO levels as done in Fig. 3, we indeed obtain only a very small electron transfer of 0.02e that is required to align the Fermi energies in this case.

C. O 2 dissociation at Al(111)
As a final example for an extended system, treated by periodic boundary conditions, we turn to the dissociation of an O 2 molecule at the Al͑111͒ surface. For this system the postulated dominant role of nonadiabatic effects 17,27,34 could not be verified until recently, since only empirical estimates of the underlying nonadiabatic PESs were available. [35][36][37] Using LC-DFT we now focus on the nonadiabatic neutral triplet PES, which is a suitable representation at large distances from the surface, where the gas-phase O 2 molecule will be in its spin-triplet ground state and the Al͑111͒ surface in a spinsinglet state. For the calculations we employed a ͑3 ϫ 3͒ Al͑111͒ slab consisting of seven aluminum layers, separated by a 30 Å vacuum. Oxygen is adsorbed at both sides of the slab to establish inversion symmetry, a real-space cutoff of 9 bohrs has been applied to the basis functions, and ten k points have been used to sample the irreducible wedge of the Brillouin zone. The electron numbers of the oxygen and the  Table I and the text for the constrained electron numbers defining the various nonadiabatic PESs. Additionally shown is the adiabatic ground-state PES, and the energy zero corresponds to the energy of the isolated neutral Cl 2 molecule and Mg 4 cluster.
aluminum subsystem used to define the neutral triplet PES are listed in Table I.
Discussing our results for the high-dimensional PES in detail elsewhere, 17,38 we illustrate the insights gained by the LC-DFT approach by concentrating on the two-dimensional dependence on the molecular bond length r and the centerof-mass distance of the molecule from the surface Z for a fixed molecular orientation and lateral position over the surface. Figure 6 shows corresponding "elbow plots" specifically for an O 2 molecule approaching the surface head-on and above an fcc site. In agreement with previous studies 27, 34 the adiabatic PES displayed in Fig. 6͑a͒ does not exhibit an energy barrier to dissociation, a finding that cannot be reconciled with the experimentally well-established low sticking coefficient for thermal molecules. 17,39 Suspecting a dominant role of nonadiabatic effects as the reason for this discrepancy, we turn to nonadiabatic representations, in which particularly the spin-triplet character of the gas-phase O 2 molecule is conserved. Figures 6͑b͒ and 6͑c͒ show corresponding PESs obtained with the FSM approach for an overall triplet state of the system and with our LC-DFT approach constraining the spin-triplet to the O 2 molecule, respectively. In contrast to the Na+ Cl neutral triplet PES discussed above, the FSM and LC-DFT approach now yield qualitatively different results. While no barrier is obtained in the prior method, the neutral triplet PES calculated with LC-DFT exhibits a clear energy barrier.
The reason for this difference is the different distribution of the magnetization density in the system. This is illustrated for a molecular configuration at the energy barrier in Fig. 7. In ͑a͒ the magnetization density computed for the free oxygen molecule ͓i.e., without the Al͑111͒ slab͔ in its spin-triplet ground state is shown, whereas ͑b͒ displays the result of an adiabatic calculation including the Al͑111͒ slab. In the latter case, neither the O 2 molecule, nor the metal atoms exhibit any spin-polarization, which is the most favored state for small molecule-surface separations. In the FSM calculation, the total spin of the system is forced to be a triplet, but as apparent from ͑c͒ the majority of this excess spin is not located at the O 2 molecule, but distributed over the entire metal slab. In contrast to this, in the LC-DFT result shown in ͑d͒ the triplet spin is localized at the oxygen molecule, reflecting the improved control over the spatial distribution of the magnetization density in the latter approach. The accumulation of spin-up density on the O 2 molecule repels the  6. ͑Color online͒ Two-dimensional cut ͑"elbow plot"͒ through the high-dimensional PES for the O 2 dissociation at Al͑111͒. The energy is shown as a function of the center-of-mass distance of the molecule from the surface Z and the oxygen-oxygen bond length r. The molecule approaches the surface head-on above an fcc site as explained in the insets. ͑a͒ Adiabatic calculation. ͑b͒ Triplet fixed-spin moment calculation. ͑c͒ Neutral triplet LC-DFT calculation. Only the latter PES exhibits an energy barrier. The energy difference between the contour lines is 0.2 eV and the small white circle denotes the molecular position discussed in Fig. 7. spin-up density of the metal slab towards the interior of the slab, while there is a strong accumulation of spin-down density close to the molecule. As a consequence, the metal slab is still in an overall singlet state in the LC-DFT calculation, which is obviously a better representation of the nonadiabatic state defined by an impinging triplet O 2 molecule compared to the FSM result. In addition, the LC-DFT approach overcomes the small charge-transfer problem present in this system, as well. Since the unoccupied 2 *↓ orbitals of the O 2 molecule are lower than the Fermi level of the metal, the adiabatic calculation yields a charge transfer of 0.01e to the O 2 molecule even at macroscopic distances from the surface.

IV. SUMMARY
In summary, we presented a locally constrained densityfunctional-theory ͑LC-DFT͒ approach, which allows one to confine electrons to subspaces of the Hilbert space, e.g., to selected atoms or groups of atoms. A major application of this technique is the computation of nonadiabatic potentialenergy surfaces, which we illustrated with examples treating the scattering of atoms and molecules at other atoms, clusters, or surfaces. Following the general formulation by Dederichs et al., 11 the electron confinement is realized by suitably introducing additional constraints to the electronic Hamiltonian. DFT is then used to obtain the fully relaxed electronic structure under the additional external potential imposed by the applied constraints. With the ⌬SCF and FSM methods as widely employed alternative implementations of this general concept, our LC-DFT method offers a more systematic approach to extended systems compared to ⌬SCF, and better control over the spatial distribution of the constraint electrons compared to FSM. This better spatial control allows one also to overcome the charge-transfer problem between widely separated subsystems that can occur in adiabatic DFT calculations.

ACKNOWLEDGMENTS
We thank Volker Blum for carefully reading the manuscript and useful suggestions. Stimulating discussions with Bengt I. Lundqvist and Eckart Hasselbrink are gratefully acknowledged.

APPENDIX: PROJECTION OPERATOR
In general, localized atom-centered basis functions such as atomic orbitals are not orthogonal, and a projection operator should be formulated in terms of covariant and contravariant basis functions. 22,23 The covariant Bloch basis functions ͉ i k ͘ and the contravariant Bloch basis functions ͉ jk ͘ are related to each other by the equations with S k being the overlap matrix. For covariant and contravariant basis functions the following orthonormality relation holds: where ␦ ij is the Kronecker symbol.
In principle there are two possible forms of the projection operator P A k into the subsystem A, which are equivalent.
and where the sum i =1, ... ,m runs over the m basis functions of subsystem A. However, expanded onto the Bloch basis functions, both these formulations for P A k yield non-Hermitian matrices. In order to facilitate the implementation into existing DFT codes, we therefore prefer to work with the following symmetrized form, which does lead to a Hermitian matrix: Inserting the expressions for the contravariant basis functions in Eq. ͑A1͒ enables us to express the projection operator entirely in terms of the known covariant basis functions ͑the atomic orbitals͒,

͑A5͒
Starting from this expression, we can then derive the matrix representation of this operator in the Hilbert space spanned by the Bloch basis functions