Fragment Approach to Constrained Density Functional Theory Calculations using Daubechies Wavelets
Laura E. Ratcliff, Luigi Genovese, Stephan Mohr, Thierry Deutsch
FFragment Approach to Constrained Density Functional Theory Calculations usingDaubechies Wavelets
Laura E. Ratcliff,
1, 2, ∗ Luigi Genovese, Stephan Mohr, and Thierry Deutsch Argonne Leadership Computing Facility, Argonne National Laboratory, Illinois 60439, USA Univ. Grenoble Alpes, CEA, INAC-SP2M, L Sim, F-38000, Grenoble, France (Dated: March 1, 2018)In a recent paper we presented a linear scaling Kohn-Sham density functional theory (DFT) codebased on Daubechies wavelets, where a minimal set of localized support functions is optimized insitu and therefore adapted to the chemical properties of the molecular system. Thanks to thesystematically controllable accuracy of the underlying basis set, this approach is able to providean optimal contracted basis for a given system: accuracies for ground state energies and atomicforces are of the same quality as an uncontracted, cubic scaling approach. This basis set offers, byconstruction, a natural subset where the density matrix of the system can be projected. In thispaper we demonstrate the flexibility of this minimal basis formalism in providing a basis set thatcan be reused as-is , i.e. without reoptimization, for charge-constrained DFT calculations within a fragment approach. Support functions, represented in the underlying wavelet grid, of the templatefragments are roto-translated with high numerical precision to the required positions and used asprojectors for the charge weight function. We demonstrate the interest of this approach to expresshighly precise and efficient calculations for preparing diabatic states and for the computational setupof systems in complex environments.
I. INTRODUCTION
Density functional theory (DFT) is arguably themost popular approach to electronic structure calcula-tions for a wide range of systems. However, it suffersfrom various well-known limitations, like for example theself-interaction problem which can result in electrondelocalization errors, and the fact that it is in principlea ground state theory only. For these reasons, the DFTformalism has been extended in the form of constrainedDFT (CDFT) to include an additional constraint onthe density, so that the lowest energy state satisfying agiven condition can instead be found. When a reasonableguess for such a condition is at hand, it can therefore beused both to find a particular excited state of the systemand to localize the electronic density in such a way asto prevent spurious delocalization, and thus provides away of overcoming the above problems. Of course, time-dependent DFT (TDDFT) can be used to find multipleexcited states, however when one is interested in a partic-ular excited state CDFT can be advantageous, especiallygiven that the additional costs associated with adding aconstraint are relatively low. Furthermore, TDDFT alsosuffers from self interaction problems, and can give inac-curate results for certain types of excited states, includingcharge transfer excitations.Constrained DFT has been implemented in a num-ber of codes, using both localized basis sets and planewaves , and has been successfully applied in a vari-ety of contexts, including charge constrained moleculardynamics , the calculation of the correct energy align-ment of metal/molecule interfaces and the calculationof electronic coupling matrix elements . For a generaloverview of CDFT see Ref. 13.As with all DFT calculations, the choice of basis sethas a large impact on both the accuracy and computa- tional cost of CDFT. One way of accessing large systemsis to reformulate the standard cubic scaling approach toDFT in terms of localized orbitals, or ‘support functions’,which we will discuss further in the following section. Wewish to perform CDFT calculations on large systems us-ing such an approach, whilst maintaining the high accu-racy associated with systematic basis sets. As such, werequire a basis set which is at the same time localizedand systematic. For this reason, we have chosen to use aDaubechies wavelet basis set , as it is a systematic basisset exhibiting the desired properties of compact supportin both real and Fourier space and can be chosen to beorthogonal. Wavelet basis sets have an inherent flexibil-ity, in that they allow for multiresolution grids, whichis particularly useful for inhomogenous systems. Com-bined with the ability to explicitly treat charged systemsin open boundary conditions, wavelets provide an idealbasis set for accurate CDFT calculations of large systems.In this paper, we show that the combination of a sup-port function approach with a wavelet basis set allowsfor the definition of a flexible fragment based approach toCDFT, which can further reduce the computational cost,particularly for very large systems. In this approach,a set of support functions are optimized for an isolated(small) molecule, or ‘fragment’, and reused as a fixed ba-sis in a larger system containing many of these moleculese.g. a solvent. It is then straightforward to associate theconstrained charge with a given fragment using a L¨owdinlike definition of the CDFT weight function. However, inthe larger system each molecule may well have a differentorientation and so the support functions, which are de-scribed in terms of the fixed wavelet grid, cannot simplybe duplicated for each molecule.Therefore, we have developed a scheme to reformat thesupport functions for arbitrary roto-translations using in-terpolating scaling functions. This interpolation, thanks a r X i v : . [ c ond - m a t . m t r l - s c i ] M a r -1.5-1.0-0.50.00.51.01.5 -6 -4 -2 0 2 4 6 8x φ (x) ψ (x) FIG. 1. Least asymmetric Daubechies wavelet family of order2 m = 16; both the scaling function φ ( x ) and wavelet ψ ( x )differ from zero only within the interval [1 − m, m ]. to the properties of the underlying basis set, results inonly a negligible loss of accuracy and so the support func-tions can be directly reused, reducing the computationalcost by an order of magnitude compared to optimizingthe support functions from scratch for the full system.In the following sections we will first summarize ourapproach to large scale DFT calculations using localizedsupport functions represented in a wavelet basis set , asimplemented in the BigDFT electronic structure code .We will then outline our implementation of CDFT, fol-lowing which we will explain our fragment approach, in-cluding a description of the reformatting scheme, vali-dating our method with calculations on prototypical sys-tems. Finally, we will present an application of CDFT forthe fullerene C in two different environments, throughwhich we will demonstrate the flexibility and potential ofa fragment based approach. II. METHODOLOGYA. Linear scaling DFT with wavelets
We and others have recently presented a newly de-veloped method for DFT calculations on large systems,which combines the use of a minimal localized basisof ‘support functions’ with the use of an underlyingwavelet basis set . This method has been implementedin BigDFT, which uses the orthogonal least asymmet-ric Daubechies family of order 16, which are depicted inFig. 1. The Kohn-Sham (KS) orbitals are expressed interms of the support functions via a set of coefficients c αi : | Ψ i (cid:105) = (cid:88) α c αi | φ α (cid:105) , (1)where the support functions are represented directly inthe wavelet basis set localized on a 3 dimensional grid,so that they can be thought of as adaptively contractedwavelets. Rather than working directly with the KS or-bitals, we instead work in terms of the density matrix, ρ ( r , r (cid:48) ), which is itself defined in terms of the supportfunctions and the density kernel, K αβ : ρ ( r , r (cid:48) ) = (cid:88) α,β φ α ( r ) K αβ φ β ( r (cid:48) ) . (2) The density matrix has been shown to decay exponen-tially with distance for systems with a gap thanks tothe so-called nearsightedness principle , and thus aformulation in terms of the density matrix allows us totake advantage of this to achieve linear scaling with thenumber of atoms in the system, thereby avoiding the cu-bic scaling of standard approaches to DFT. From thisthe charge density is calculated directly from the sup-port functions and density kernel. Similarly, the bandstructure energy and the charge of the system can becalculated from the density kernel via: E BS = Tr [ KH ] , N = Tr [ KS ] (3)where H indicates the Hamiltonian matrix in the basisof the support functions, and S is the support functionoverlap matrix.The support function formalism allows one to map thedegrees of freedom of KS orbitals into a localized descrip-tion, that can be directly put in relation with atomic po-sitions. In practice, the support functions are truncatedwithin spherical localization regions with a user-definedradius, and some additional truncation must be appliedto the density kernel which is then exploited via sparsematrix algebra to achieve a fully linear scaling algorithm.Therefore, to some extent, a given support function φ α can be associated to the atom a where its localizationregion is centered. In order to achieve accurate results,both the support functions and density kernel are opti-mized during the calculation, that is the energy is min-imized with respect to both quantities. Providing thelocalization regions are sufficiently large, this results in aminimal localized basis set with an accuracy equivalentto the underlying basis set.The general scheme is common to other basis optimiza-tion for density-matrix minimization based linear scal-ing DFT codes, e.g. ONETEP and Conquest , withthe addition of a few novel features. These include theapplication of a confining potential to the KS Hamilto-nian, which ensures the support functions remain local-ized. In order to apply the confining potential consis-tently, we also enforce an approximate orthogonality con-straint on the support functions, in contrast with otherapproaches which use fully non-orthogonal support func-tions . Furthermore, the properties of the descriptionin the wavelet basis are such that the algorithm guar-antees that the Pulay contribution to the atomic forcescan safely be neglected, so the forces can be calculatedaccurately and cheaply.The method can be divided into two key components:the optimization of the support functions (either withor without a confining potential), and the optimizationof the density kernel. This latter point can be achievedvia a choice of schemes, as detailed previously . Theseinclude a direct minimization approach, where the coeffi-cients c αi are first updated using DIIS or steepest descentsto minimize the band structure energy and then used toconstruct the density kernel, and the Fermi Operator Ex-pansion (FOE) method, where the density kernel is ex-pressed as a function of the Hamiltonian matrix whichcan be evaluated numerically using a Chebyshev polyno-mial expansion. Close attention has also been paid tothe parallelization of the code, such that massively par-allel machines can be exploited to perform large scalecalculations (see also Ref. ). For charged calculations,we have found the use of the direct minimization methodto be the most suitable, due to a reduction in the oc-currence of charge sloshing during convergence, and theflexibility afforded by working with the wavefunction co-efficients rather than directly with the density kernel asin the FOE method. B. Atomic charge analysis
The mapping between electronic and localized degreesof freedom which is provided by the support function for-malism allows one to perform an accurate atomic chargeanalysis, meaning that each atom is assigned a partialnet charge, such that the electrostatic properties of thesystem are conserved. Obviously this conservation isonly possible within a certain limit, as one is mappinga continuous quantity (the electronic charge) to a dis-crete quantity (the atomic point charges). If, however,the error introduced by this mapping onto point chargesis small enough, the system under investigation can bereasonably approximated by a simple setup of chargedpoint particles, which paves the way for future applica-tions such as coupling different levels of accuracy withinthe same calculation.Given the overlap matrix S and the density kernel K ,the partial charge located on atom a can be defined bythe so-called L¨owdin charge: q a = ( a ) (cid:88) α (cid:48) (cid:16) S / KS / (cid:17) α (cid:48) α (cid:48) , (4)where the sum runs over all support functions α (cid:48) whichare located on atom a . Obviously (cid:80) a q a = tr( KS ) = N ,i.e. the total charge (the monopole) is conserved. In orderto check whether higher multipoles can also be conserved,we compared the dipole moment calculated using this ap-proach with that calculated using the continuous chargedensity. In addition a comparison with the dipole mo-ment calculated with the cubic scaling version of BigDFTwas done as a reference. The values for a strongly polar-ized molecule (H O) and a non-polarized one (C ) aregiven in Tab. I.As a second test of the reliability of our method wedirectly compared the atomic point charges with thosecalculated by performing a Bader charge analysis of thecharge density calculated using the cubic scaling ap-proach. As can be seen from Tab. II, the differencesbetween the exact results are smaller with our approach.In particular, for the C fullerene, reasons of symmetryimpose that the charge should be equally distributed andno atom should carry a net charge. As can be seen, the L¨owdin procedure comes closer to this result than theBader analysis. C. Constrained DFT
The general idea of constrained DFT is to force acharge to remain localized in a given region of the simu-lation space. This is achieved via the addition of a La-grange multiplier term to the Kohn-Sham energy func-tional which enforces a given constraint on the resultingelectronic density, so that rather than being the ground-state density of the system, the density instead corre-sponds to a particular excited state. This Lagrange mul-tiplier can also be thought of as an additional appliedpotential, otherwise referred to as the constraining po-tential. The constraint can also take a number of otherforms, but for the purposes of this work we are inter-ested only in constraining the charge. The new functionaltherefore becomes: W [ ρ, V c ] = E KS [ ρ ] + V c (cid:18)(cid:90) w c ( r ) ρ ( r ) d r − N c (cid:19) , (5)where E KS is the Kohn-Sham energy functional, V c is theaforementioned Lagrange multiplier, N c is the requiredcharge within the specified region and w c ( r ) is a weightfunction which defines this region. The weight functionand N c are defined in advance, however the value of V c which correctly enforces the constraint must be foundduring the calculation. Wu and Van Voorhis demon-strated that the lowest energy state for which the con-straint is correctly applied is in fact a maximum withrespect to V c and so it becomes possible to efficientlydetermine the correct V c . It is also straightforward toadd multiple constraints to the system, and indeed oneis frequently interested in constraining the charge differ-ence between two regions. This feature is important forthe simulation of charge-transfer excitations within theCDFT formalism.Rewriting the new functional (Eq. 5) in density matrixform as : W [ ρ, V c ] = E KS [ ρ ] + V c (Tr [ Kw c ] − N c ) , (6)the charge constraint is easily added to the existing algo-rithm in BigDFT. This construction requires the weightmatrix, w αβ , which is defined via the weight function as w αβ = (cid:90) φ α ( r ) w c ( r ) φ β ( r ) d r . (7)It then remains to define the weight function, for which anumber of different schemes exist. The support functionapproach used here lends itself to a L¨owdin like defini-tion, which is analogous to that used above to determineatomic charges. Using this approach we directly con-struct the weight matrix via w αβ = (cid:16) S PS (cid:17) αβ , (8) cubic linearexact dipole exact dipole point charge approx.H O (0.463, -0.506, -0.186) (0.466, -0.510, -0.187) (0.606, -0.668, -0.247)norm 0.711 0.716 0.935 d ex · d (-0.0004, -0.0004, -0.0004) (-0.025, -0.025, -0.025) (-0.055, -0.055, -0.055)TABLE I. Dipole moments calculated using the exact charge density for the cubic and linear scaling approaches, respectively,and using the partial atomic point charges. All values are given in atomic units.cubic – Bader linear – L¨owdinH O (-1.27, 0.61, 0.67) (-0.83, 0.42, 0.41)C . ± .
045 0 . ± . O we in-dicate the values on all three atoms, for C we give the meanof the absolute values together with the standard deviation.All values are given in atomic units. where S is the overlap matrix between support functionsand P is a projection matrix, defined as 1 for all supportfunctions belonging to the region where a constraint isbeing applied, and 0 elsewhere. Alternatively, if one isconstraining a charge difference between two regions, itshould be set to 1 on one of the regions, − V c for a given charge constraint value, N c . There are twopossible approaches to the optimization. In the first ap-proach, one can find the optimum value of V c at each stepof the self-consistent density optimization, i.e. the groundstate density is updated in an outer loop, with V c updatedin an inner loop. Alternatively, the second approach con-sists of fully minimizing the functional W [ ρ, V c ] of Eq. 6with respect to the density for a fixed value of V c , up-dating V c and finding the new minimum density, thenrepeating to convergence, i.e. the maximization with re-spect to V c is performed in an outer loop with the groundstate density found self-consistently in an inner loop. Wechose the latter, as it was observed to be more stable.We use Newton’s method to update V c , with the secondderivative calculated using a finite difference approach. D. Fragment approach
The combination of the novel features described aboveand the use of a wavelet basis set make this approachideal for the application of CDFT to large systems. Inparticular the ability to reuse the support functions canresult in significant savings for e.g. geometry optimiza-tions and calculations on charged systems, as previously support functionoptimization support functionreformatting
FIG. 2. The fragment approach as illustrated for a cluster ofwater molecules: the support functions are initially optimizedfor an isolated water molecule and then duplicated for a col-lection of water molecules, avoiding the need for optimizationin the larger system. demonstrated . Furthermore, this idea of support func-tion reuse can be extended to a fragment based approach,which is similar to the so-called fragment orbital methodwhich has been used to calculate electronic coupling ma-trix elements .The central idea is to take a group of atoms, or morespecifically an isolated molecule, and fully optimize thesupport functions. These support functions are then usedas a fixed basis for a system containing several molecules,as illustrated in Fig. 2 for a simple example. We refer tothe initial molecule for which the support functions wereoptimized as the ‘template’ molecule.As the support functions are kept fixed in the frag-ment approach, w αβ need only be calculated once at thestart of the calculation, after which it remains fixed. Fur-thermore, due to the quasi-orthogonality of the supportfunctions, when the fragment approximation is justified S can in general be calculated using a Taylor approxi-mation, and so the calculation of the weight matrix addsvery little overhead to the calculation.In this work we focus on systems where the respectivefragments are well defined, and thus the support func-tions generated from the isolated fragments can be usedfor the full system with a minimal impact on the ac-curacy. In cases where electrons are being added to afragment, it is important to ensure that the lowest unoc-cupied molecular orbital (LUMO) and, if necessary, thenext few states in energy are sufficiently well representedby the support function basis. As discussed elsewhere ,this can be achieved using the direct minimization for-malism to optimize a few additional states during the iso-lated calculation without adding a charge to the system.For systems where the fragments are less well defined theimplementation could in principle be extended to furtheroptimize the support functions for the combined system,either in a neutral state or while the charge constraintis being enforced. In such cases, the L¨owdin approach isalso expected to be less accurate, and so it would be de-sirable to use an alternative form for the weight function.Finally, it should be mentioned that there are somesubtleties related to the initial guess for the charge den-sity. This depends on the initial density kernel, whichis constructed from the fragment KS orbitals. For neu-tral calculations it is straightforward to use the fragmentorbitals and occupancies directly from the isolated calcu-lations, however for charged calculations some additionalinput is required. One approach would be to occupy thefragment orbitals in order of their energies, however thiscan lead to charge distributions which are significantlydifferent from the required constraint. This can resultin slow convergence, or even worse, problems with localminima. A better approach should therefore take intoaccount the effect of the constraining potential on thefragment orbital energies. This can be done by assigningoccupation numbers so that any excess/deficit in chargeis localized on the same fragment as the constraint, sothat the initial density already satisfies the charge con-straint. Alternatively, the risk of encountering a localminimum can be reduced by adding a degree of noiseto the fragment orbitals, or by completely randomizingthe initial guess, subject to the correct overall charge.However, such an approach is in general much slower toconverge, and thus the latter strategy is generally notrecommended. E. Reformatting scheme for roto-translations
As the support functions are defined in terms of anunderlying grid of wavelets, in order to implement afragment approach it is necessary to have a scheme forreformatting the support functions due to a change inatomic positions. We are frequently interested in situa-tions where a molecule has been rotated and translated,for example when calculating electronic coupling matrix elements in a dimer for varying angles between the twomonomers. Therefore, we have developed and imple-mented a scheme for reformatting the support functions,given the axis and angle of rotation between some initialand final positions for a given fragment mass center.A fragment of the system is defined by the user viaa list of atomic positions, which should of course be inbijection with the atom list defining the template frag-ment. Therefore the first problem is to identify the com-bination of translation and rotation which sends the tem-plate fragment to the position of the system’s fragment.As a first step, two reference systems are chosen suchthat the fragment center of mass is in the same position.This operation is equivalent to finding the translation be-tween the template and the system. We then have twolists of atomic positions, { R T,Sa } , where T and S labelthe template and system fragment respectively, and thesubscript a , indicating the atom, ranges from 1 to N , thenumber of atoms in the fragment.If the system fragment is a rigid displacement of thetemplate (i.e. its internal coordinates are unchanged),the rotation matrix R we seek is such that R Sa = (cid:80) Nb =1 R ab R Tb . In general, we should assume there is aslight modification of the internal coordinates, as the ge-ometry of the fragment might be affected by the interac-tion with the environment. In this case, the matrix R issuch as to minimize the cost function J ( R ) = 12 N (cid:88) a =1 || R Sa − N (cid:88) b =1 R ab R Tb || . (9)The determination of the matrix R in such a mannerconstitutes a version of the well-known Kabsch algo-rithm (also known as Wahba’s problem) , which canbe solved by the Singular-Value Decomposition of the3-by-3 matrix B ij = (cid:80) Na =1 ( R Sa ) i ( R Ta ) j . After havingfound two matrices U and V and a diagonal matrix S such that B = USV t , the optimal rotation is R = UDV t D = diag (1 , , det( U ) det( V )) . (10)The value of J ( R ) defined in (9) might then be usedto quantify the validity of the rigid transformation ap-proach. In the case where its value is below a giventhreshold (fixed to 10 − in our case), we may proceedwith the reformatting of the template basis functions,which will be denoted by (cid:12)(cid:12) φ Tα (cid:11) in what follows.As described in Refs. , from the expression of (cid:12)(cid:12) φ Tα (cid:11) in a Daubechies wavelets basis set, the so-called“magic-filter” transformation can be used to define a realspace representation of the basis functions, given in termsof one-dimensional interpolating scaling functions (ISF) φ Tα ( x, y, z ) = (cid:88) i,j,k c ijk ϕ i ( x ) ϕ j ( y ) ϕ k ( z ) , (11)where ϕ i ( x ) ≡ ϕ ( x/h − i ) is one element of the ISF ba-sis set, which is constituted of uniform translations ofthe mother function ϕ ( t ) over the points of a uniformgrid of spacing h , covering the entire simulation domain.These points are labeled by indices ijk . The points r ijk = ( hi, hj, hk ) therefore lie within the box contain-ing the support of φ Tα ( r ).This real-space expression is optimal in the sense thatit preserves the same moments of the original represen-tation given in Daubechies wavelets. The interpolatingproperty of the ISF basis set is such that c ijk = φ Tα ( r ijk ).Let us suppose we have a one-dimensional function ex-pressed in ISF, namely f ( x ) = (cid:80) i f i ϕ i ( x ). We knowthat f i = f ( hi ). If we want to translate the function f by a displacement ∆, and express this function in theISF basis, we have f ( x + ∆) = (cid:80) i f (cid:48) i ϕ i ( x ), with f (cid:48) i = f ( hi + ∆) = (cid:88) j f i − j t ∆ j , (12)where the filter t ∆ j = ϕ ( j + ∆ /h ) implements the (uni-form) translation. This filter has a limited extension (thesame as the function ϕ ( x )) and of course t hkj = δ j, − k .Imagine now we have a different ISF basis set { ϕ I (˜ x ) } defined on a uniform grid spacing of separation ˜ h anda reference frame ˜ x I = ˜ hI , which is related to x by amore complicated transformation ˜ x ( x ) of the coordinatespace. If this transformation can be inverted, by x (˜ x ),then a new function ˜ f (˜ x ) ≡ f ( x (˜ x )) can be defined in thisframe. For each grid point I , it is then always possible tofind ¯ i in the old frame such as to minimize the absolutevalue of ∆ I ≡ x (˜ x I ) − h ¯ i .Using the above relations we might approximate˜ f (˜ x ) (cid:39) (cid:80) I ˜ f I ϕ I (˜ x ), where˜ f I = ˜ f (˜ x I ) = f ( h ¯ i + ∆ I ) = (cid:88) j f ¯ i − j t ∆ I j . (13)If the transformation ˜ x is a continuous function of x which varies slowly enough, this is in general a rathergood approximation (see Fig. 3).This framework can be easily generalized to a roto-translation in three dimensions. Indeed, we would like toestimate the function φ Sα (˜ r ) ≡ φ Tα ( r (˜ r )) (cid:39) (cid:88) I,J,K ˜ c IJK ϕ I (˜ x ) ϕ J (˜ y ) ϕ K (˜ z ) , (14)where the coordinates ˜ r = (˜ x, ˜ y, ˜ z ) are defined as˜ r = R · r (15)where R is calculated by Eq. (10). In addition a rigidshift vector s = ( s x , s y , s z ) is defined as the differencebetween the coordinates of the center of mass of the twofragments. If the rotation is the identity matrix, thetemplate reference frame is then ˜ r = r + s . As in theone dimensional case presented above, the interpolationdepends on the inverse mapping r (˜ r ). We detail in thefollowing a procedure to identify such a function.The coefficients ˜ c IJK of φ Sα (˜ r ) can be found in threesteps. We first start by considering the transformation law for ˜ x . This transformation can be thought of as afunction of the template coordinates r :˜ x ( x, y, z ) = R x + R y + R z . (16)In the same spirit as Eq. (13), we may invert Eq.(16) withrespect to one template coordinate t = x, y, z into ˜ x inthe system’s reference frame. The choice of the variable t depends on the entries of the rotation matrix, and it isin general given by the coordinate which is multiplied bythe coefficient of the highest absolute value in Eq. (16).This choice guarantees that the t variable is the one forwhich ˜ x − t is slowly varying. Let us imagine t = x forthis example. We can define the function φ (1) α (˜ x, y, z ) = φ Tα ( x (˜ x, y, z ) − s x , y, z )= (cid:88) I,j,k ˜ c I,j,k ϕ I (˜ x ) ϕ j ( y ) ϕ k ( z ) , (17)by proceeding for all j, k , as described in Eq.(13), to de-fine the coefficients ˜ c I,j,k . The second step is related tothe expression of ˜ y . Depending on the choice of the vari-able t in the first step, we have to consider one of thesethree relations: R ˜ y = R ˜ x + R y − R z , (18) R ˜ y = R ˜ x − R x + R z , (19) R ˜ y = R ˜ x + R x − R y , (20)which hold when in the first step t = x, y, z respectively.These relations can be derived using the orthogonality ofthe rotation matrix R . This function can now be invertedwith respect to one of the old variables. Again, this choicewill depend on the values of the coefficients multiplyingeach variable.In our example, we have to consider the relation (18)as we have chosen t = x in the first step. We choose toinvert the relation with respect to z , having z = z (˜ x, ˜ y, y ).In this case we will have, as a second step φ (2) α (˜ x, ˜ y, y ) = φ (1) α (˜ x, y, z (˜ x, ˜ y, y ) − s z )= (cid:88) I,J,j ˜ c I,j,J ϕ I (˜ x ) ϕ J (˜ y ) ϕ j ( y ) . (21)In the third step, the remaining variable, (which is y forthe illustrated example), can be directly obtained fromthe inverse relation r = R − · ˜ r = R t · ˜ r , (22)which is easier to express as R is an orthogonal matrix.In our case, the final result is therefore φ Sα (˜ x, ˜ y, ˜ z ) = φ (2) α (˜ x, y (˜ x, ˜ y, ˜ z ) − s y , ˜ z ) . (23)We recall that the definition of φ (1 , α depends on the or-der of the operation. Here we have chosen to interpolatefirst with respect to x , then z and y . The best choice oforder depends only on the entries of the matrix R .
1. Accuracy
In order to assess the accuracy of the reformattingscheme, we have applied it to a water molecule under-going a series of rotations. Support functions were gen-erated for a template water molecule, using a dense gridwith a spacing of 0 .
132 ˚A, and were then reused for watermolecules in a variety of different orientations using a lessdense grid with a spacing of 0 .
185 ˚A. As a point of com-parison, calculations were also performed for each orien-tation fully optimizing the support functions with a gridspacing of 0 .
185 ˚A. This allows us to quantify both theerror introduced by the support function reformattingand the errors due to representing the wavefunctions ona fixed grid, i.e. the so-called ‘eggbox effect’. The eggboxeffect of the standard cubic scaling approach is also pre-sented. The computational setup has been chosen suchthat the difference in ground state energies between thecubic and support function approaches is of the order of1 meV/atom.The results are shown in Fig. 3, where we can see thatthe eggbox effect is of the order of 0.1 meV/atom. Asboth the cubic and linear scaling approaches use the sameunderlying grid, the variation is similar in each case. Theerror due to the interpolation also remains small – lessthan a few meV/atom. Importantly, the overall error forthe reformatted calculations remains of the same order ofmagnitude as that due to the selected localization radiiof the support functions.
III. RESULTS
Below we present results for three different systems,where the first two can be validated against CDFT im-plementations in other codes. In each case we use thelocal density approximation (LDA) exchange-correlationfunctional and HGH pseudopotentials within isolatedboundary conditions. For carbon, nitrogen and oxygenwe use four support functions per atom, for hydrogen weuse one, and for zinc we use nine. A. N Wu and Van Voorhis have previously studied N andso this system provides a useful test case. We used agrid spacing of 0 .
185 ˚A with support function radii of7 . E − E r e f ( m e V / a t o m ) θ ( ° ), u(0,0,1) (0,1,0) (0,1,1) (1,0,0) (1,0,1) (1,1,0) (1,1,1)cubic eggboxlinear eggboxtemplate FIG. 3. Plot showing the energy variation for a watermolecule rotated through different angles ( θ ) and axes of ro-tation ( u ). Results are shown for the standard cubic scalingapproach (‘cubic eggbox’), fully optimized support functions(‘linear eggbox’) and a fixed support function basis generatedfor a template molecule (‘template’). The cubic reference isthe energy at the initial orientation calculated using the cu-bic approach, for the linear and template approaches it is thesame quantity calculated in the fully optimized support func-tion basis. There is a roughly constant error of 1 meV/atomin the support function basis compared to the cubic scalingapproach. Selected orientations are shown along the bottom. -0.4-0.3-0.2-0.10.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 V c N c Δ E ( e V ) FIG. 4. Change in energy with respect to an unconstrainedcalculation and applied potential value for differing chargeseparations in N . used B3LYP we do not expect exact agreement. Fur-thermore, we should recall that the fragment approachpresented here is aimed at systems where the donor andacceptor are well separated, whereas the use of supportfunctions optimized for an isolated nitrogen atom is nec-essarily an approximation in this case. Nonetheless, wehave successfully reproduced the correct trends for boththe energy and the Lagrange multiplier. FIG. 5. The ZnBC-BC model complex.
B. ZnBC
As a more elaborate test case we take thezincbacteriochlorin-bacteriochlorin (ZnBC-BC) complex,which has also been studied previously in some detail,both with CDFT and other approaches, e.g. Refs. 36and 37. This system is ideally suited to our approach asthe donor and acceptor are clearly separated. Further-more, TDDFT has been shown to give incorrect energiesfor the ZnBC + -BC − and ZnBC − -BC + charge transfer(CT) excited states and so the advantages of CDFTare clear. It has previously been demonstrated that thedifferences between a (1,4)-phenylene-linked ZnBC-BCcomplex and a model complex where the link is elimi-nated are small ; for simplicity we therefore choose touse the latter, where the distance between the two previ-ously linked carbon atoms is 5.84˚A, as depicted in Fig. 5.Taking the coordinates from Ref. 37, we relaxed the iso-lated ZnBC and BC molecules separately, then built themodel complex without further relaxation. We used agrid spacing of 0 .
185 ˚A and localization radii of 5 .
82 ˚A. Toassess the accuracy of the fragment support functions wecompare the neutral energies for the model complex withthose obtained using cubic scaling BigDFT. The resultsare shown in Tab. III, where we can see that the errorfor both the model complex and the isolated molecules isless than 1 meV/atom.The energies for the two CT excited states relativeto the unconstrained DFT ground state are 3.71 eV forZnBC + -BC − and 3.98 eV for ZnBC − -BC + , which isconsistent with previous results . We can also gainsome insight into the nature of these CT states by plot-ting the difference in the electronic density between theneutral and constrained calculations, as in Fig. 6. Notonly is the charge transfer characteristic clear, the plotfor ZnBC + -BC − also shows remarkably good agreementwith previous calculations that used the significantlymore expensive Bethe-Salpeter approach , which con-firms that CDFT can be used to obtain physically rel-evant CT excitons, and provide a reliable estimation ofthe corresponding excitation energies.We have also plotted the relationship between the con-straining potential, V c , the total energy relative to theunconstrained calculation, ∆ E , and the charge differencebetween the two molecules, N c . This is shown in Fig. 7, cubic frag. diff.(eV) (meV)ZnBC − . − .
604 22 . − . − .
980 17 . − . − .
575 54 . − -BC + - − .
592 -ZnBC + -BC − - − .
860 -TABLE III. Energies for isolated ZnBC and BC, the neu-tral model ZnBC-BC complex and the two lowest energy CTstates, as calculated using standard BigDFT (‘cubic’) and thefragment approach (‘frag.’). Where applicable the differencebetween the two approaches is also indicated (‘diff.’). (a)ZnBC − -BC + (b)ZnBC + -BC − FIG. 6. Density differences between the neutral and chargedcalculations for the two charge transfer states. Red (blue) in-dicates an increase (decrease) in the electronic charge densitywith respect to the neutral. where N c = 1 corresponds to ZnBC + -BC − and N c = − − -BC + ; our results agree well withprevious calculations , despite the use of a differentexchange-correlation functional. This test highlights therobustness of the method – in order for the correct valueof V c to be found within a minimal number of iterationsof the constraint loop, there should be a smooth rela-tionship between a given V c and the resulting N c . If forcertain values of V c the convergence is insufficient, suchthat the final charge deviates from the correct value, thiswill negatively impact the search for the correct V c . Weobserved that in general such a smooth curve is straight-forward to obtain, given a reasonable initial guess for thedensity kernel and therefore charge density. As discussedin Section II C, this can be achieved by defining the ini- -0.10.00.1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 V c N c Δ E ( e V ) FIG. 7. Applied potential value and change in energy com-pared to an unconstrained calculation for differing charge sep-arations in the ZnBC-BC model complex. tial occupancies in a manner which is consistent with thedesired charge difference.
C. C In order to accurately calculate material properties itis important to account for environmental effects, e.g. byincluding a solvent or neighboring molecules in a molec-ular material. However, this can considerably increasethe cost of a simulation, as in the case of large systemsin solution where the solvent must fill a correspondinglylarge volume. Various strategies have been developed forreducing the cost, for example by using implicit solvationmethods , however it is frequently desirable to treatexplicitly the environmental degrees of freedom. Thanksto the fragment approach, the treatment of solvents andother surrounding molecules can readily be achieved inBigDFT with relatively low cost, as we will demonstratethrough the example of the fullerene C in two differentenvironments: when in an aqueous solution and whensurrounded by other C molecules. For each system weconstrain a charge of ± molecule inorder to determine the environmental impact on the ion-ization potential (IP) and electron affinity (EA).For traditional DFT calculations with semilocal func-tionals like LDA it is well known that the above quanti-ties are badly estimated by the frontier orbitals, i.e. theHOMO (LUMO) for the IP (EA). Therefore, in order toextract physically meaningful information, one must ei-ther use more expensive beyond-DFT approaches, or in-stead calculate the IP and EA using the so-called ∆SCFmethod. This is made possible when explicitly chargedcalculations are available, i.e. when charged and neutralcalculations have energies that can be measured withrespect to a common reference. The treatment of theelectrostatic potential which is included in the BigDFTcode makes such a comparison possible . This latter ap- proach results in values which match experiment muchbetter than traditional semilocal functionals, indeed ourresults for the IP and EA of the isolated molecule agreevery well with the experimental values of 7.6 eV and2.7 eV respectively. In this case, we wish to apply the∆SCF approach to a molecule in an environment, which,as we will show, is easily achieved using CDFT.On the other hand, with unconstrained DFT calcula-tions the use of the ∆SCF approach is much more deli-cate when studying environmental effects with LDA: asthe charge tends to be overly delocalized, charged cal-culations do not simply represent a perturbation fromthe isolated values, as discussed in more detail below. Inother words, the calculated energy differences do not cor-respond to the IP and EA of C in an environment, butto a completely different quantity. If one wishes to calcu-late this quantity it is therefore essential to use CDFT.
1. Computational details
There have been a number of previous studies of C inwater, both experimental and theoretical , howeverthey have mainly focused on neutral fullerenes. Previousresearch has indicated the existence of a first hydrationshell surrounding C containing between 60 and 65 wa-ter molecules , we have therefore chosen to restrictourselves to systems containing 66 water molecules. Wepresent results for three example structures, which aredepicted in Fig. 8(a). They were generated by insertingthe C into water droplets where the water moleculeswere deposited with random orientations at random po-sitions subject to the room-temperature density of wa-ter. The structures were then relaxed until the RMSforces were below 10 meV/˚A. For the environment offullerenes, we limit the cost of the simulations by includ-ing six nearest neighbor fullerenes only, so that the sys-tem is arranged as a three dimensional cross, as depictedin Fig. 8(b). Each of the fullerenes was considered in itsgas-phase structure.The fragment calculations were performed with a gridspacing of 0 .
185 ˚A, while the template calculations wereperformed using a denser grid of 0 .
132 ˚A to ensure ac-curate reformatting; we used support function radii of4 .
23 ˚A. These values have been chosen such as to en-sure the applicability of the L¨owdin approach for theweight matrix on the central C whilst preserving ab-solute accuracy of the unconstrained calculations to theorder of 3 meV/atom, see data in Sec. III C 2. In orderto ensure the support functions are sufficiently accuratefor the negatively charged calculations for C , the tem-plate calculation was performed optimizing three addi-tional states (to account for degeneracies).We continue to use the LDA functional as the differ-ence between the LDA and other treatments like PBE forthe IP and EA of fullerenes has previously been shownto be negligible, and reasonable agreement with exper-iment has also been observed . We also neglected the0 (a)The three configurations A, B and C (left to right) of C in H O. (b)C with its six nearest neighbors.The blue lines are drawn between themolecular centers along the axes. FIG. 8. The different environments for C . modelling of dispersive terms on the C – C interac-tions due to their negligible impact on frontier orbitaleigenvalues and on total energy differences in chargedcalculations.
2. Testing the fragment approach
The results for the C structure with a center to cen-ter distance of 10 ˚A are shown in Tab. IV with the cor-responding values for the isolated molecule. We also in-clude the cubic scaling results as they allow us to assessthe accuracy of the fragment approach for this system.As anticipated, for the isolated molecule the fragment er-ror is of the order of 0.2 eV, which is about 3 meV/atom.In order to confirm that this accuracy is preserved, wealso compared unconstrained cubic and fragment calcu-lations for the seven C structure. To find an uncon-strained solution we built the initial guess from the frag-ment densities in such a manner that the charge wasequally distributed among fragments; the final solutionremained close to this charge distribution. We found val-ues of -3.681 eV and 6.707 eV, i.e. the difference with thecubic results is of the same magnitude as for the isolatedmolecule (see Tab. IV).We have also investigated the effect of varying theseparation between the molecules by repeating the con-strained fragment and (unconstrained) cubic calculationswith center to center separations ranging from 10 ˚A to20 ˚A, which corresponds to a shortest distance betweenmolecules of 3.1 ˚A to 13.1 ˚A. The results are plotted inFig. 9. As expected, for large separations the results tendtowards the isolated values. For the unconstrained calcu-lations there is an abrupt change between two differentstates, whereas with CDFT there is not only a smoothtrend, but also an exponential relationship with distance,proving that the fragment approach is sufficiently preciseto capture such trends. isolated in H O in C (10 ˚A) Q cubic frag. A B C cdft cubic − − . − . − . − . − . − . − . .
648 7 .
783 7 .
262 8 .
033 7 .
837 7 .
526 6 . E Q − E forC when isolated and in the two environments. Two valuesare given for the isolated C : that of the fragment approach,which in this case merely refers to a fixed support functionbasis as only one fragment is present, and the cubic scalingreference. For the results in water, constrained fragment re-sults are given. For the nearest neighbor results (‘in C ’),results are presented for both the constrained fragment and(unconstrained) cubic approaches. The unconstrained resultsexhibit stronger deviation from the isolated values, showingthat the environment is not correctly modeled as it is notacting as a perturbation of the system. Units are in eV. Thus far, we have only considered calculations witha shifting of the template support functions, however wealso wish to demonstrate the effectiveness for rotated sup-port functions. To this end, for a distance of 10 ˚A thesix outer fullerenes were collectively rotated along the z -axis by angles of 15 ◦ , 45 ◦ and 90 ◦ , with the orientationof the central molecule remaining unchanged. This wasfound to have a negligible impact on the IP and EA, witha difference in the values for the various orientations ofaround 0 .
01 eV for the constrained fragment calculationscompared with 0 .
05 eV for the unconstrained cubic calcu-lations. Such values are too small to be significant com-pared to the errors associated with the basis. In order todetermine whether the energies are truly unaffected bythe orientation, it would be necessary to account for thedispersion effects which are not captured by the LDA.However, the fact that no spurious errors are introducedby the rotating of the support functions serves to furtherconfirm the accuracy of the reformatting scheme.1 Δ E ( e V ) R (Å) cubic, Q=−1frag. cdft, Q=−1cubic, Q=+1frag. cdft, Q=+1
FIG. 9. Variation in electron affinity and ionization potentialfor increasing separations, where the distance is measured be-tween the centers of neighboring molecules. The energy plot-ted is relative to the isolated value in the respective basis, i.e.∆ E = ( E Q isol. − E ) − ( E Q full − E ), where Q indicates thecharge state, ‘full’ refers to the 7 molecule system and ‘isol.’refers to the isolated molecule.
3. Comparison of environments
We now return to the results in Tab. IV in order tocompare the effect of the two environments. As expected,for both environments the constrained values remain rel-atively close to the isolated results, certainly much closerthan the unconstrained results for the fullerene environ-ment. In a sense, when the constraint is enforced, thepresence of the surrounding molecules could be thoughtof as a perturbation on the isolated state, although thestrength of the perturbation is clearly much stronger forthe water. To further explore this, we have also plotteddifferences in the converged electronic densities betweenthe neutral and charged calculations in Fig. 10. The ef-fects of the charge constraint are clearly visible, with theexcess charge distributed across the molecules for the cu-bic calculation and much more clearly localized for theconstrained calculations. Furthermore, the charge dif-ference on the central molecule for the constrained cal-culations clearly retains the same character as the re-spective isolated density, with the excess or deficit ofcharge also resulting in an induced dipole on the neigh-boring molecules. As expected, the impact of the wateris stronger than the neighboring fullerenes, where thecloser proximity and stronger dipole moment of the wa-ter molecules results in a stronger deviation from the iso-lated density difference. Similar behavior has also beenobserved for the water structures which are not depicted.As can be seen from the variation in IP and EA forstructures A, B and C (Tab. IV), the effect of the water isnot only stronger than that of the neighboring fullerenes,but the dependence on the structure is also quite sig-nificant. Furthermore, we have also observed that theresulting energies are strongly affected by the choice of weight function, which is not the case for the fullereneenvironment. There is some freedom in the procedurefor optimizing the template support functions, e.g. thelocalization radii and the number of additional states in-cluded; we tested a few of the different options. For thefullerene environment the variation in calculated IP andEA due to the choice of template parameters were smalland systematic, whereas for the aqueous environment thevariation was much stronger. Indeed, the fragment ap-proach provides an ideal setup to explore the impact ofdifferent weight functions – such a strong dependencewould be harder to detect when considering only two orthree different choices. In future the choice of weightfunction could also be decoupled from the fragment basisto allow for a more thorough exploration of its influencein constraining the charge.Of course, in order to correctly assess the impact ofthe environment, one should go beyond the model struc-tures used here, both in terms of the size and procedureused to generate them. Furthermore, in the case of thewater, proper sampling should be performed over a num-ber of different configurations, which should be generatedat the correct temperature e.g. using molecular dynam-ics or Monte Carlo simulations. However, asidefrom the generation of input structures and any even-tual relaxations of the atomic coordinates, the fragmentcalculations are quick and easy to perform, requiring lit-tle additional setup aside from the template calculations.Furthermore, we have demonstrated both the accuracyand flexibility of the fragment approach for such systems.As such, given appropriate atomic coordinates this workcould easily be extended in future to a large number ofconfigurations for both environments, or indeed appliedto other fullerenes or solvents.
IV. CONCLUSION
We have presented a method for constrained DFTcalculations on large systems, using a fragment basedscheme. This has been implemented in the BigDFT elec-tronic structure code within a recently developed frame-work, which uses a basis of localized support functionsrepresented in an underlying wavelet grid to achieve lin-ear scaling behavior with respect to system size whileretaining the systematic accuracy of the underlying grid.The division of a given system into fragments (ideallydistinct molecules), each with its own associated sup-port functions leads to a natural approach to CDFT,where the charge is constrained to a given fragment viaa L¨owdin like definition of the weight function. ThisL¨owdin approach can also be used to straightforwardlycalculate atomic charges, as we have demonstrated.Furthermore, by using a reformatting scheme whichenables the reuse of support functions for identical frag-ments, irrespective of their position or orientation in thesystem, we are able to further reduce the cost of simula-tions by an order of magnitude, as the support functions2 (a)isolated fragment, Q = − Q = − Q = − Q = − Q = +1 (f)cubic, Q = +1 (g)constrained fragment, Q = +1 (h)constrained fragment, Q = +1 FIG. 10. Density differences between the neutral and charged calculations for the fullerene when isolated, when surrounded byother fullerenes (with a center to center distance of 10 ˚A) and when in water (structure A). The densities are plotted on thecentral plane with a logscale, with red (blue) indicating an increase (decrease) in the electronic charge density with respect tothe neutral. can be separately optimized for each ‘template’ fragmentand used as a fixed basis in the system of interest. Theproperties of the wavelet basis set ease the implementa-tion of reformatting a numerical field given in real space,so that we were able to implement this scheme in a man-ner which is both efficient and accurate. The flexibilityof this method, together with the ability of the BigDFTcode to treat systems with many atoms (see e.g. Ref. ),makes it ideally suited for both neutral and charged cal-culations on very large systems at modest computationalcost.We have presented results from two previously stud-ied systems in order to validate our approach, as well asan example application. For this latter point, we haveperformed calculations on C in two different environ-ments, namely a model nearest neighbor system contain-ing seven fullerenes and in an aqueous solution. The ef-fects of the constraint are clearly visible in the electronicdensities, which we have compared to the unconstrainedand isolated results. We have also shown that the pres-ence of water has both a stronger impact on the resultsand a stronger dependence on the choice of weight func-tion than the presence of neighboring fullerenes.The reformatting approach described here has anotherkey benefit aside from reducing the cost of such calcu- lations: as the basis set for each fragment will remainequivalent following the reformatting, the computationalsetup provided by our approach is ideal where Hamilto-nian matrix elements of the whole system have to beconsidered. For example, for electronic coupling ma-trix elements (‘transfer integrals’) between two identicalmonomers, the basis set for each monomer will remainequivalent following the reformatting, so that there is noambiguity in the sign of the coupling matrix elements.In contrast, for support functions optimized from scratchthere is no guarantee that the phase for both the supportfunctions and wavefunction coefficients will be identicalbetween the two molecules and thus the sign of the trans-fer integral cannot be determined. Indeed, we will bepublishing results in the near future for such an appli-cation based on the framework presented in the currentwork .In the future we also hope to extend this work to per-mit calculations on realistic nanoscale devices, using amulti-scale approach. In the first instance this wouldinvolve changing the definition of a fragment to an in-dividual atom, which naturally leads to a DFT basedtight-binding like method. This formalism also allowsthe correct definition of an ‘embedded’ approach, wheredifferent regions of a simulation cell are treated at differ-3ent levels of precision, e.g. for a point defect in a bulksemiconductor, with higher accuracy close to the defect.Work is ongoing in this direction. ACKNOWLEDGEMENTS
We acknowledge funding from the European projectMMM@HPC (RI-261594), the CEA-NANOSCIENCE BigPOL project and the ANR projects SAMSON (ANR-AA08-COSI-015) and NEWCASTLE (ANR-2010-COSI-005-01). This research used resources of the ArgonneLeadership Computing Facility at Argonne National Lab-oratory, which is supported by the Office of Science of theU.S. Department of Energy under contract DE-AC02-06CH11357. CPU time was also provided by IDRIS(project i2014096905). ∗ lratcliff@anl.gov P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964) W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965) J. P. Perdew and A. Zunger, Phys. Rev. B , 5048 (1981) Y. Zhang and W. Yang, J. Chem. Phys. , 2604 (1998) P. H. Dederichs, S. Bl¨ugel, R. Zeller, and H. Akai, Phys.Rev. Lett. , 2512 (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997(1984) A. M. P. Sena, T. Miyazaki, and D. R. Bowler, J. Chem.Theory Comput. , 884 (2011) H. Oberhofer and J. Blumberger, J. Chem. Phys. ,064101 (2009) J. ˇRez´aˇc, B. L´evy, I. Demachy, and A. de la Lande, J.Chem. Theory Comput. , 418 (2012) A. M. Souza, I. Rungger, C. D. Pemmaraju, U. Schwingen-schloegl, and S. Sanvito, Phys. Rev. B , 165112 (2013) Q. Wu and T. Van Voorhis, J. Chem. Phys. , 164105(2006) H. Oberhofer and J. Blumberger, J. Chem. Phys. ,244105 (2010) B. Kaduk, T. Kowalczyk, and T. Van Voorhis, Chem. Rev. , 321 (2012) I. Daubechies,
Ten Lectures on Wavelets (SIAM, 1992) S. Mohr, L. E. Ratcliff, P. Boulanger, L. Genovese, D. Cal-iste, T. Deutsch, and S. Goedecker, J. Chem. Phys. ,204110 (2014) L. Genovese, A. Neelov, S. Goedecker, T. Deutsch,S. A. Ghasemi, A. Willand, D. Caliste, O. Zilberberg,M. Rayson, A. Bergman, and R. Schneider, J. Chem. Phys. , 014109 (2008) W. Kohn, Phys. Rev. Lett. , 3168 (1996) E. Prodan and W. Kohn, PNAS , 11635 (2005) S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett. , 2127(1999) W. Kohn, Int. J. Quantum Chem. , 229 (1995) L. He and D. Vanderbilt, Phys. Rev. Lett. , 5341 (2001) R. Baer and M. Head-Gordon, Phys. Rev. Lett. , 3962(1997) C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C.Payne, J. Chem. Phys. , 084119 (2005) D. R. Bowler and T. Miyazaki, J. Phys.: Condens. Matter , 074207 (2010) S. Mohr, L. E. Ratcliff, L. Genovese, D. Cal-iste, P. Boulanger, S. Goedecker, and T. Deutsch,arXiv:1501.05884(submitted) Q. Wu and T. Van Voorhis, Phys. Rev. A , 024502 (2005) K. Senthilkumar, F. C. Grozema, F. M. Bickelhaupt, andL. D. A. Siebbeles, J. Chem. Phys. , 9809 (2003) K. Senthilkumar, F. C. Grozema, C. F. Guerra, F. M. Bick-elhaupt, F. D. Lewis, Y. A. Berlin, M. A. Ratner, andL. D. A. Siebbeles, J. Am. Chem. Soc. , 14894 (2005) W. Kabsch, Acta Crystallogr. A , 827 (Sep 1978) G. Wahba, SIAM Rev. , 409 (1965) F. L. Markley, J. Astronaut. Sci. , 245 (1988) A. Neelov and S. Goedecker, J. Comp. Phys. , 312(2006) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. , 566(1980) C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev. B , 3641 (1998) Q. Wu and T. Van Voorhis, J. Chem. Theory Comput. ,765 (2006) I. Duchemin, T. Deutsch, and X. Blase, Phys. Rev. Lett. , 167801 (2012) A. Dreuw and M. Head-Gordon, J. Am. Chem. Soc. ,4007 (2004) J. Tomasi and M. Persico, Chem. Rev. , 2027 (1994) D. A. Scherlis, J.-L. Fattebert, F. Gygi, M. Cococcioni,and N. Marzari, J. Chem. Phys. , 074103 (2006) J. Fosso-Tande and R. J. Harrison, Chem. Phys. Lett. , 179 (2013) J. Dziedzic, H. H. Helal, C.-K. Skylaris, A. A. Mostofi, andM. C. Payne, Europhys. Lett. , 43001 (2011) A. Cerioni, L. Genovese, A. Mirone, and V. A. Sole, J.Chem. Phys. , 134108 (2012) D. Lichtenberger, K. Nebesny, C. Ray, D. Huffman, andL. Lamb, Chem. Phys. Lett. , 203 (1991) J. Devries, H. Steger, B. Kamke, C. Menzel, B. Weisser,W. Kamke, and I. Hertel, Chem. Phys. Lett. , 159(1992) J. A. Zimmerman, J. R. Eyler, S. B. H. Bach, and S. W.McElvany, J. Chem. Phys. , 3556 (1991) X.-B. Wang, C.-F. Ding, and L.-S. Wang, J. Chem. Phys. , 8217 (1999) R. Rivelino and F. de Brito Mota, Nano Lett. , 1526(2007) R. Rivelino, A. M. Maniero, F. V. Prudente, and L. S.Costa, Carbon , 2925 (2006) P. Scharff, K. Risch, L. Carta-Abelmann, I. Dmytruk,M. Bilyi, O. Golub, A. Khavryuchenko, E. Buzaneva,V. Aksenov, M. Avdeev, Y. Prylutskyy, and S. Durov, Car-bon , 1203 (2004) L. Li, D. Bedrov, and G. D. Smith, J. Chem. Phys. ,204504 (2005) Y. I. Prylutskyy, A. S. Buchelnikov, D. P. Voronin,V. V. Kostjukov, U. Ritter, J. A. Parkinson, and M. P.Evstigneev, Phys. Chem. Chem. Phys. , 9351 (2013) S. Banerjee, J. Chem. Phys. , 044318 (2013) N. Choudhury, J. Chem. Phys. , 034502 (2006) M. L. Tiago, P. R. C. Kent, R. Q. Hood, and F. A. Re-boredo, J. Chem. Phys. , 084311 (2008)55