Table Of Content

An exact Coulomb cutoff technique for supercell calculations Carlo A. Rozzi,1,2 Daniele Varsano,3,2 Andrea Marini,4,2 Eberhard K. U. Gross,1,2 and Angel Rubio1,3,2 1Institut fu¨r theoretische Physik, Freie Universität Berlin, Arnimallee 14, D-14195 Berlin, Germany 2European Theoretical Spectroscopy Facility (ETSF) 3Departamento de F´ısica de Materiales, Facultad de Ciencias Qu´ımicas, UPV/EHU, Centro Mixto CSIC-UPV/EHU and Donostia International Physics Center, E-20018 San Sebastián, Spain 4Dipartimento di Fisica, Università “Tor Vergata”, Roma, Italy 6 (Dated: 4th February 2008) 0 0 We present a new reciprocal space analytical method to cutoff the long range interactions in 2 supercell calculations for systems that are infinite and periodic in 1 or 2 dimensions, extending n previous works for finite systems. The proposed cutoffs are functions in Fourier space, that are a used as a multiplicative factor to screen the bare Coulomb interaction. The functions are analytic J everywhere but in a sub-domain of the Fourier space that depends on the periodic dimensionality. 2 We show that the divergences that lead to the non-analytical behaviour can be exactly cancelled when both theionic and theHartree potential are properly screened. This techniqueis exact, fast, ] andveryeasytoimplementinalreadyexistingsupercellcodes. Toillustratetheperformanceofthe r newscheme,weapplyittothecaseoftheCoulombinteractioninsystemswithreducedperiodicity e h (as one-dimensional chains and layers). For those test cases we address the impact of the cutoff in t differentrelevantquantitiesforgroundandexcitedstateproperties,namely: theconvergenceofthe o ground state properties, the static polarisability of the system, the quasiparticle corrections in the . t GW scheme and in the binding energy of the excitonic states in the Bethe-Salpeter equation. The a results are very promising. m - PACSnumbers: 02.70.-c,31.15.Ew,71.15.-m,71.15.Qe d n o I. INTRODUCTION ofsystems with the same kind ofreduced periodicity are c theclassesofthepolymers,andofthesolidswithdefects. [ Throughout this paper we call nD-periodic a 3D ob- 1 Plane waves expansions with periodic boundary con- ject, thatcanbe consideredinfinite andperiodic inn di- v ditions have been proven to be a very effective way mensions,beingfiniteintheremaining3 ndimensions. 1 to exploit the translational symmetry of infinite crystal − In order to simulate this kind of systems, a commonly 3 solids,inordertocalculatethepropertiesofthebulk,by adopted approach is the supercell approximation.1 0 performing the simulations in one of its primitive cells 1 In the supercell approximation the physical system is only1. The use of plane waves is motivated by several 0 treated as a fully 3D-periodic one, but a new unit cell facts. First, the translational symmetry of the poten- 6 (the supercell) is built in such a way that some extra 0 tials involved in the calculations is naturally and easily empty space separates the periodic replica along the di- / accounted for in reciprocal space, through the Fourier t rection(s) in which the system is to be considered as fi- a expansion. Second, veryefficient and fast algorithms ex- nite. Thismethodmakespossibletoretainalltheadvan- m ist (like FFTW2) that allow us to calculate the Fourier tages of plane waves expansions and periodic boundary - transforms very efficiently. Third the expansionin plane conditions. Yet the use ofa supercellto simulate objects d waves is exact, since they form a complete set, and it is n that are not infinite and periodic in all the directions, only limited in practice by one parameter, namely, the o leads to artifacts, even if a very large portion of vacuum maximum value of the momentum, that determines the c is interposed between the replica of the system in the : size of the chosen set. Fourth, in many cases, the use v non-periodic dimensions. of Born-von Karmán periodic boundary conditions sup- i Infact,thestraightforwardapplicationofthesupercell X pliesaconceptuallyeasy(thoughartificial)waytogetrid method generates in any case fake images of the original r of the dependence of the properties of a specific sample system, that can mutually interact in several ways, af- a on its surface and shape, allowing us to concentrate on fecting the results of the simulation. It is well known the bulk properties of the system in the thermodynamic that the response function of an overall neutral solid of limit.3 molecules is not equal, in general, to the response of the However, mainly in the last decade, increasing inter- isolated molecule, and converges very slowly to it, when est has been developed in systems at the nano-scale,like the amount of vacuum in the supercell is progressively tubes, wires, quantum-dots, biomolecules, etc., whose increased.7,8 physicaldimensionalityis, for allpracticalpurposes,less For instance, the presence of higher order multipoles than three.4 These systems are still 3-dimensional (3D), can make undesired images interact via the long range buttheirquantumpropertiesarethoseofaconfinedsys- part of the Coulomb potential. In the dynamic regime, teminoneormoredirections,andthoseofaperiodicex- multipolesarealwaysgeneratedbytheoscillationsofthe tendedsystemintheremainingdirections. Otherclasses charge density even in systems whose unit cell does not 2 carryany multipole in its groundstate. This is the case, In a finite system the potential is usually required to forexample,whenweinvestigatetheresponseofasystem be zero at infinity. In a periodic system this condition in presence of an external oscillating electric field. ismeaningless,sincethe systemitself extendsto infinity. Thingsgoworsewhentheunitcellcarriesanetcharge, NeverthelessthegeneralsolutionofEq. (1)inbothcases since the total charge of the infinite system represented is known in the form of the convolution by the supercell is actually infinite, while the charges at n(r′) the surfaces of a finite, though very large system always V(r)= d3r′, (2) generate a finite polarisationfield. This situation is usu- r r′ sZpZacZe | − | ally normalised in the calculation by the introduction of a suitable compensating positive backgroundcharge. that it is referred from now on as the Hartree potential. Another common situation in which the electrostat- It might seem that the most immediate way to build ics is known to modify the ground state properties of thesolutionpotentialforagivenchargedistributionisto the system occurs when a layeredsystem is studied, and computetheintegralinrealspace,butproblemsimmedi- an infinite array of planes is considered instead of a sin- atelyariseforinfinitesystems. Infact,thedensitycanbe gle slab, being in fact equivalent to an effective chain of reducedtoaninfinitesumoverdeltachargedistributions capacitors.9 qδ(r r′), These issuesbecome particularlyevidentinallthe ap- − proaches that imply the calculation of non-local opera- q V(r)= , (3) tors or response functions, because, in these cases, two r L n supercells may effectively interact even if their charge Xn | − | densitiesdonotoverlapatall. Thisisthecase,forexam- and the integral in Eq. (2) becomes an infinite sum as ple, of the many-body perturbation theory calculations well, but this sum is in general only conditionally, and (MBPT), and, in particular, of the self-energy calcula- notabsolutelyconvergent.10 ThesumofEq.(3)apoten- tions at the GW level.6,7 tialthatis determineduptoa constantfora neutralcell However we are usually still interested in the disper- with zero dipole moment, while the corresponding sum sionrelationsoftheelementaryexcitationsofthesystem for the electric field is absolutely convergent. A neutral along its periodic directions, and those are ideally dealt cellwithanonnulldipolemoment,ontheopposite,gives with using a plane waves approach. Therefore, the ideal a divergent potential, and an electric field that is deter- path to keep the advantagesof the supercellformulation minedup to anunknownconstantelectric field(the sum inplanewaves,andtogainadescriptionofsystemswith for the electric field is conditionally convergent in this reduced periodicity free of spurious effects is to develop case). a technique to cut the Coulomb interaction off out of a Even if, in principle, the surface terms have to be al- desired region. This problem is not new and has been ways taken into account, in practice they are only rel- addressednowforaverylongtime andindifferentfields evant when we want to calculate energy differences be- (condensed matter, classicalfields, astrophysics,biology, tweenstateswithdifferenttotalcharge. Thesetermscan etc). Severaldifferent approacheshave been proposedin be neglected in the case of a neutral cell whose lowest the past to solve it, however a complete review of them nonzero multipole is quadrupole11. As in the present isbeyondthescopeofthispaper. Theaimofthepresent work we are interested in macroscopic properties of the work is to focus on the widely used supercell schemes to periodicsystem,thosesurfaceeffectsareneverconsidered show how the image interactioninfluences both the elec- inthediscussionthatfollows. Howeverthissample-shape tronicgroundstatepropertiesandthedynamicalscreen- effects play an important role for the analysis of differ- ingintheexcitedstateof0D-,1D-,2D-periodicsystems, ent spectroscopies as, for example, infrared and nuclear and to propose an exact method to avoid the undesired magnetic resonance. interaction of the replicas in the non-periodic directions. A major source of computational problems is the fact The paper is organised as follows: in Sec. II the ba- thatthesuminEq. (3)isveryslowlyconvergingwhenit sics of the plane wave method for solids are reviewed, in issummedinrealspace,andthisfacthashistoricallymo- Sec.III the new methodis outlined, inSec.IVthe treat- tivatedtheneedforreciprocalspacemethodstocalculate ment of the singularities is explained, in Sec. V some it. It was Ewald who first discovered that, by means of applications of the proposed technique are discussed. anintegraltransform,the sumcanbe splitintwoterms, and that if one is summed in real, while the other in reciprocal space, both of them are rapidly converging.12 II. THE 3D-PERIODIC CASE The point of splitting is determined by an arbitrary parameter. The main problem of electrostatics we are facing here Letusnowfocusonmethods ofcalculatingthe sumin canbereducedtothe problemoffinding solutionstothe Eq. (3) purely based on the reciprocal space. Poissonequationforagivenchargedistributionn(r),and If we consider a periodic distribution of charges with given boundary conditions density n(r) such that n(r) = n(r + L ), with L = n n 2V(r)= 4πn(r). (1) n L ,n L ,n L , and n ,n ,n Z, it turns out x x y y z z x y z ∇ − { } { } ∈ 3 that the reciprocal space expression for a potential like known methods, tipically used in molecular dynamics simulations, are the multipole-correction method21, and V(r)= n(r′)v(r r′ )d3r′, (4) the particle-mesh method22, whose review is beyond the | − | scope ofthe presentwork,andwe refer the readerto the ZZZ space original works for details. Differently from what happens for the Ewald sum, in a 3D-periodic system, can be written as the method that we propose to evaluate the sum in V(G )=n(G )v(G ), (5) Eq. (3) entirely relies on the Fourier space and amounts n n n to screening the unit cell from the undesired effect of where we have used the convolution theorem to trans- (some of) its periodic images. The basic expression is form the real space convolution of the density and the Eq. (5), whose accuracy is only limited by the max- Coulomb potential into the product of their reciprocal imum value GN of the reciprocal space vectors in the space counterparts. Here G = n G ,n G ,n G sum. Since there is no splitting between real and recip- n x x y y z z are the multiples of the primitive {reciprocal space vec}- rocal space, no convergence parameters are required. tors 2π, and v(G ) is the Fourier transform of the long Our goal is to transform the 3D-periodic Fourier rep- Ln n resentation of the Hartree potential of Eq. (5) into the range interaction v(r), evaluated at the point G . For n modified one the Coulomb potential it is V˜(G )=n˜(G )v˜(G ) (8) n n n 4π v(G )= . (6) n G2 such that all the interactions among the undesired peri- n odicreplicaofthesystemdisappear. Thepresentmethod Fourier transforming expression (5) back into real is a generalisation of the method proposed by Jarvis et space we have, for a unit cell of volume Ω, al.23 for the case of a finite system. In order to build this representation, we want to: 1) V(r)= 4Ωπ n(GG2n)exp(iGn·r). (7) dsyesfitneem,aosuctreoefnwinhgicrhegthioenreDisanrooCunoduloeamcbhicnhtearragcetiionn;th2e) nX6=0 n calculate the Fourier transform of the desired effective interactionv˜(r) that equals the Coulombpotential in , At the singular point nx = ny = nz = 0 the potential D and is 0 outside V isundefined,but,sincethevalueatG=0corresponds D totheaveragevalueofV inrealspace,itcanbechosento be anynumber,correspondingtothe arbitrarinessinthe 1 if r V˜(r)= r ∈D . (9) choice of the static gauge (a constant) for the potential. (0 if r / Observe that the same expression can be adopted in the ∈D case of a charged unit cell, but this time, the arbitrary Finallywemust3)modifythedensityn(r)insuchaway choice of v(G) in G = 0 corresponds to the use of a that the effective density is still 3D-periodic, so that the uniform backgroundneutralising charge. convolution theorem can be still applied, but densities belonging to undesired images are not close enough to interact through v˜(r). III. SYSTEMS WITH REDUCED PERIODICITY The choice of the region for step 1) is suggested by D symmetryconsiderations,anditisasphere(orradiusR) It has been shown13 that the slab capacitance effect for finite systems, an infinite cylinder (of radius R) for mentionedinthe introductionactuallyis a problemthat 1D-periodic systems, and an infinite slab (of thickness cannot be solved by just adding more vacuum to the su- 2R) for 2D-periodic systems. percell. This has initially led to the development of cor- Step 2) means that we have to calculate the modified rections to Ewald’s original method14, and then to rig- Fourier integral orous extensions in 2D and 1D.15,16 The basic idea is to V˜(G)= v˜(r)e−iG·rd3r= v(r)e−iG·rd3r. restrictthe suminreciprocalspacetothe reciprocalvec- tors that actually correspond to the periodic directions sZpZacZe ZDZZ of the system. These approach are in general of order (10) O(N2)17,18, but they have been recently refined to order Still we have to avoid that two neighbouring images in- O(NlnN).19,20 Another class of techniques, developed teract by taking them far away enough from each other. so far for finite systems, is based on the expansion of Then step 3) means that we have to build a suitable su- the interaction into a series of multipoles (fast multipole percell, and re-define the density in it. method).24,25,26 Withthistechniqueitispossibletoeval- Let us examine first step 2), i.e. the cutoff Coulomb uateeffectiveboundaryconditionsforthePoisson’sequa- interaction in reciprocal space. We know the expression tionsatthecell’sboundary,sothattheuseofasupercell of the potential when it is cutoff in a sphere.23 It is is not required at all, making it computationally very 4π efficient for finite24,25 and extended systems26. Other v˜0D(G)= 1 cos(GR) . (11) G2 − (cid:2) (cid:3) 4 ThelimitR convergestothebareCoulombtermin does not exist, since for G = 0, the cutoff has a finite z thesenseofa→di∞stribution,while,sincelim v˜0D(G)= value, while it diverges in the limit G 0 G→0 k → 2πR2, there is no particular difficulty in the ori- 1 gin. This scheme has been successfully used in many v˜2D(G ,G ) for G >0, G 0+. (15) applications10,23,25,27,28. k z ∼ G2k k z → The 1D-periodic case applies to systems with infinite So far we haven’t committed to a precise value of the extent in the x direction, and finite in the y and z di- cutoff length R. This value has to be chosen, for each rections. The effective Coulomb interaction is then de- dimensionality, in such a way that it avoids the interac- fined in real space to be 0 out of a cylinder of radius tion of any two neighbour images of the unit cell in the R having its axis parallel to the x direction. By per- non-periodic dimension. forming the Fourier transformation we get the following In orderto fix the values of R we must choosethe size expressionfor the cutoff coulombpotentialin cylindrical ofthesupercell. Thisleadsustothestep3)ofourproce- coordinates:31 dure. Werecallthatevenoncethelongrangeinteraction iscutoffoutofsomeregionaroundeachcomponentofthe 4π v˜1D(G ,G )= 1+G RJ (G R)K (G R) system, this is not sufficient yet to avoid the interaction x ⊥ G2 ⊥ 1 ⊥ 0 | x| (cid:20) among undesired images. The charge density has to be modified, or, equivalently, the supercell has to be built G RJ (G R)K (G R) , (12) −| x| 0 ⊥ 1 | x| in such a way that two neighbouring densities along ev- (cid:21) ery non periodic direction do not interact via the cutoff whereJ andK aretheordinaryandmodifiedcylindrical interaction. Bessel functions, and G = G2 +G2. It is easy to see how this could happen in the simple ⊥ y z case of a 2D square cell of length L: if both r and r′ It is easy to realise that, sqince the K functions damp belongstothecell,thenr,r′ L,and r r′ √2L(see the oscillations of the J functions very quickly, for all ≤ | − |≤ the schematic drawing in Fig. 1). If a supercell is built practical purposes this cutoff function only acts on the that is smaller than (1+√2)L, there could be residual very first smaller values of G, while the unscreened 4π G2 interaction, and the cutoff would no longer lead to the behaviour is almost unchanged for the larger values. exact removal of the undesired interactions. Unfortunately, while the J (ξ) functions have a con- n Let us call A the unit cell of the system we are work- stantvalueforξ =0,andthewholecutoffiswelldefined 0 ing on, and = A ,i = , , the set of all the for G⊥ = 0, the K0(ξ) function diverges logarithmically cells in the sAystem{. iIf the−sy∞ste·m··is∞nD}-periodic this set for ξ 0. Since, on the other hand, K (ξ) ξ−1 for small→ξ, 1 ≈ only includes the periodic imagesof A0 inthe n periodic directions. Let us call the set of all the non-physical B images of the system, i.e. those in the non-periodic di- v˜1D(G ,G ) log(G R) for G >0, G 0+. x ⊥ ∼− x ⊥ x → rections. Then = R3. Obviously, if the system is (13) 3D-periodic =AR∪3,Band ∅. This means that the limit limG→0+v(G) does not exist IngeneralAwewantto alBlow≡the interactionofthe elec- for this cutoff function, and the whole G = 0 plane x trons in A with the electrons in all the cells A , is ill-defined. We will come back to the treatment of the 0 i ∈ A but not with those B . To obtain this we define the singularitiesinthenextsection. Wenoticethatthisloga- i ∈B supercell C A such that, i rithmicdivergenceisthecommondependenceonewould 0 ⊇ 0 ∀ get for the electrostatic potential of a uniformly charged if r A , and r′ A then r r′ C 1Dwire29. It is expected that bringingchargeneutrality ∈ 0 ∈ i | − |∈ 0 (16) in place would cancel this divergence (see below). (if r∈A0, and r′ ∈Bi then |r−r′|∈/ C0 The 2D-periodic case, with finite extent in the z di- (see Fig. 1 for a simplified 2D sketch). The new density rection, is calculated in a similar manner. The effective n˜(r) is such that Coulomb interaction is defined in real space to be 0 out of a slab of thickness 2R symmetric with respect to the if r A then n˜(r)=n(r), xy plane. In Cartesian coordinates we get ∈ 0 (17) (if r C0, and r / A0 then n˜(r)=0. ∈ ∈ v˜2D(Gk,Gz)= G4π2 1+e−GkR|GGz|sin(|Gz|R) TtiohnessdizeepeLnCdsofonthtehesuppeerri-ocdeilcl idnimthenesnioonna-pliteyrioodfitchdeisryecs-- (cid:20) k tem. In order to completely avoid any interaction, even e−GkRcos(Gz R) , (14) in the case the density of the system is not zero at the − | | (cid:21) cell border, it has to be where G = G2 +G2. L =(1+√3)L for finite k x y C In the limqit R the unscreened potential 4π is LC =(1+√2)L for 1D-periodic (18) recovered. Similar→ly t∞o the case of 1D, the limit GG2 0 LC =2L for 2D-periodic →  5 (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) twheekcnoonwstianntasdcvoannsciesttehnattlythfoersuthmemtwuostpboetefinntiitaels,.wSeicnacne (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) iprnarcnolvgueiddepe,dasrotthtoaoftsbwpoeetahfikn,pdaoltlaetnmhteieatilnhsfioondnitttioehseseinspataomratehtefeosooeutcitnotngh.setalnotnsg, (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) expInlowithedattfoololobwtasiwnethsheoewxahcotwcachnacerglleatnioeuntsrawlihteyncaonpebre- L ating with the cutoff expression of Sec. III in Fourier (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) space. Thetotalpotentialofthe systemisbuiltinthefollow- ingway: weseparateoutfirstshortandlongrangecontri- (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) bapuoGtteiaonuntsisasitlaogntehncheeariarogtneeidcdpbeyontsteithnyitsiandl+eb(nyrs)iat=yddiZisnVgexap(n(rd−)sa=u2brZ2tr)ea.rcf(Ttairnh)ge. + r The ionic potential is then written as (1+ (2))L erf(ar) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) V(r)=∆V(r) Z , (19) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)p(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) where a is chosen so that ∆−V(r) isrlocalised within a (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) psprhesesrieonofofratdhieusiornaic, psmotaelnletriatlhiannretchieprcoeclallsiszpea.ceTihse ex- (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) V(G)=4π+∞rsinG(Gr)∆V(r)dr− exp(cid:16)G−24Ga22(cid:17) Z 0 Figure 1: Schematic description for the supercell construc-  (20) tion in a 2D system. The upper sketch corresponds to the which,forG=0givesafinite contributionfromthefirst 2D-periodic case (i.e. a 2D crystal). The middle sketch cor- term,anda divergentcontributionfromthe secondterm responds to a 0D-periodic system, and the bottom one to a 1D-periodic. In the 0D-periodic case the systems in different cells do not interact, while in the 1D-periodic the chains do +∞ notinteract,butwithineachchaintheinteractionofallofits V(G=0)=4π r2∆V(r)dr . (21) elements is permitted.. −∞ Z 0 The first is the contribution of the localised charge, and Actually, since the required super-cell is quite large, a iseasilycomputed,sincetheintegrandiszeroforr >r . a compromisebetweenspeedandaccuracycanbeachieved ThesecondtermiscancelledbythecorrespondingG=0 in the computation, using parallelepiped super-cell with termintheelectronicHartreepotential,duetothecharge L = 2L for all the cases. This approximation rests C neutralityofthesystem. Thistriviallysolvestheproblem on the fact that the charge density is usually contained of the divergences in 3D-periodic systems. in a region smaller than the cell in the non-periodic di- Nowletusconsidera1D-periodicsystem. TheHartree rections, so that the spurious interactions are, in fact, term alone in real space is given by avoided,evenwithasmallercell. Therefore,onthebasis ofthisapproximation,wecanchoosethevalueofthecut- off length R always as half the smallest primitive vector V(x,y,z)= n(Gx,y′,z′) in the non-periodic dimension. XGx ZΩZ v(Gx,y y′,z z′)eiGxx′dy′dz′. (22) × − − IV. CANCELLATION OF THE SINGULARITIES Invoking the charge neutrality along the chain axis, we have that the difference between electron and ionic den- Themainpointinthe procedureofeliminating thedi- sities satisfies vergencesinallthecasesofinterestistoobservethatour final goal is usually not to obtain the expression of the [n (G =0,y,z) n (G =0,y,z)]dydz =0. Hartree potential alone, because all the physical quan- ion x − el x ZZ tities depend on the total potential, i.e. on the sum of (23) the electronic and the ionic potential. When this sum is Unfortunately, the cutoff function in Eq. (12) is diver- considered we can exploit the fact that each potential is gent for G = 0. So the effective potential results in an x defined up to an arbitrary additive constant, and choose undetermined 0 form. However, we can work out an ·∞ 6 analyticalexpressionforitbydefiningfirstafinite cylin- G 0D-periodic v˜0D(G)= drical cutoff, but then bringing the size of the cylinder >0 4π[1−cos(GR)] to infinity. This way, as a first step, we get a new cutoff G2 interactionin a finite cylinder of radius R, and length h, 0 2πR2 assuming that h is much larger than the cell size in the periodic direction. In this case the modified finite cutoff Gx G⊥ 1D-periodic v˜1D(Gx,G⊥)= potential includes a term G4π2 1+G⊥RJ1(G⊥R)K0(GxR) >0 any v˜1D(G ,r,h) log h+√h2+r2 , (24) (cid:2)−GxRJ0(G⊥R)K1(GxR) x ∝ r ! 0 >0 −4π 0RrJ0(G⊥r)log(r)dr (cid:3) which, in turn gives, for the particular plane Gx =0, 0 0 −πRR2(2log(R)−1) R Gk Gz 2D-periodic v˜2D(Gk,Gz)= v˜1D(Gx =0,G⊥)≈−4π rJ0(G⊥r)log(r)dr >0 any G4π2 1+e−GkR GGkz sin(GzR)−cos(GzR) Z h (cid:16) (cid:17)i 0 J1(G⊥R) 0 >0 G4π2z [1−cos(GzR)−GzRsin(GzR)] +4πRlog(2h) . (25) G⊥ 0 0 −2πR2 The effective potential is now split into two terms, but Table I: Reference Table summarising the results of the cut- only the second one depends on h. The second step is offworkforcharge-neutralsystems: finitesystems(0D),one- achieved by going to the limit h + , to obtain the dimensionalsystems(1D)andtwo-dimensionalsystems(2D). → ∞ exact infinite cutoff. By calculating this limit, we no- The complete reciprocal space expression of the Hartree po- tice that only the second term in the right hand side tential is provided. For the 1D case, R stands for the radius of Eq. (25) diverges. This term is the one that can be of the cylindrical cutoff whereas in the 0D case is the radius droppeddue tochargeneutrality(infactithasthe same of thespherical cutoff. In 2D stands for half thethickness of theslab cutoff (see text for details). formforthe ionicandelectronicchargedensities). Thus, for the cancellation to be effective in a practical implementation, we have to treat on the same way both the ratio Lx between the in-plane lattice vectors. ionic and Hartree Coulomb contributions. Of course the Ly first term in the right hand side Eq. (25) has always to 4π be taken into account, affecting both the long and the v˜2D(G =0,G ) [1 cos(G R) G Rsin(G R)] short range part of the cutoff potentials. k z ≈ G2z − z − z z Following this procedure, we are able to get a consid- (α+√1+α2)(1+√1+α2) sin(G R) z +8hlog . erable computational advantage, compared, e.g., to the α G ! z methodoriginallyproposedbySpataruet al.,30 sinceour (27) cutoffisjustananalyticalfunctionofthereciprocalspace coordinates, and the evaluation of an integral for every The G=0 value is value of G ,G is not needed. The cutoff proposed in x ⊥ Ref.30isactuallyaparticularcaseofourcutoff,obtained v˜2D(G =0,G =0)= 2πR2 (28) k z − by using the finite cylinder for all the components of the Gvectors: inthis casethe quadratureinEq.(25) has to To summarise, the divergences can be cancelled also beevaluatedforeachG ,G ,andG ,andaconvergence in 1D-periodic and 2D-periodic systems provided that x y z study in h is mandatory (see discussion in Sec. VB, and 1) we apply the cutoff function to both the ionic and Fig. 5). the electronic potentials, 2) we separate out the infinite contributionasshownabove,and3)weproperlyaccount In the 1D-periodic case, the G = 0 value is now well for the short range contributions as stated in Table I. defined, and it turns out to be lim v˜(G ,G ) G⊥→0 x ⊥ The analyticalresults of the present workare condensed in Tab. I: all the possible values for the cutoff functions v˜1D(G =0,G =0)= πR2(2log(R) 1). (26) x ⊥ − − are listed there as a quick reference for the reader. The analogous result for the 2D-periodic cutoff is obtained by imposing finite cutoff sizes h = αh = h V. RESULTS x y (much larger than the cell size), in the periodic direc- tionsx andy,anddroppingtheh-dependentpartbefore The scheme illustrated above has been implemented passing to the limit h + . The constant α is the both in the real space time-dependent DFT code → ∞ 7 OCTOPUS27, and in the plane wave many-body- 60 perturbation-theory (MBPT) code SELF32. The tests Ionic Hartree have been performed on the prototypical cases of in- Total 40 finite chains of atoms along the x axis. The compar- isons are performed between the 3D-periodic calculation ) V 20 (physicallycorrespondingtoacrystalofchains),andthe e ( 1D-periodic case (corresponding to the isolated chain) al both in the usual supercell approach,and within our ex- nti 0 e act screening method. The discussion for the 2D cases ot P-20 follows the same path as for the 1D case, while results for the finite systems have already been reported in the literature.23,25 We addressed different properties to see -40 the impactofthe cutoffateachlevelofcalculation,from the ground state to excited state and quasiparticle dy- -60-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 namics. y axis (a.u.) A. Ground state calculations 80 Allthecalculationshavebeendonewiththereal-space 60 implementation of DFT in the OCTOPUS27 code. We have used non-local norm-conserving pseudopotentials34 40 ) to describe the electron-ion interaction and the local- V e density approximation (LDA)35 to describe exchange- al ( 20 correlation effects. The particular choice of exchange- nti correlationorionic-pseudopotentialdoesnotmatterhere ote 0 as we want to assert the impact of the Coulomb cutoff P and this is independent of those quantities. Moreover, -20 Ionic Hartree we have used a grid of 0.38 a.u. for Si and Na. Total In this case the footprint of the interaction of neigh- -40 bouring chains in the y and z direction is the dispersion -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 of the bands in the corresponding direction of the Bril- y axis (a.u.) louin zone. However it is known that, if the supercell is largeenough,thebandsalongtheΓ X directionareun- − Figure 2: Calculated total and ionic and Hartree potentials changed. This is in apparent contradiction with the fact for a 3D-periodic (top) and 1D-periodic (bottom) Si chain that the radial ionic potential for a wire (that asymp- . totically goes like ln(r) as a function of the distance r from the axis of the wire) is completely different from the crystal potential. potential now behaves like it is expected for a potential The answerto this contradictionis clearif we perform of a chain, i.e. diverges logarithmically, and is clearly a cutoff calculation. In fact the overall effect on the oc- different from the latter case. Nevertheless the sum of cupied states turns out to be cancelled by the Hartree the ionic and Hartree potential is basically the same as potential, i.e. by the electron screening of the ionic po- for the 3D-periodic system. tential, but two different scenarios are visible as soon as In the static case the two band structure are then ex- the proper cutoff is used. pected, and are found to be the same, confirming that, In Fig. 2 (top) it is shown the ionic potential, the as far as static calculations are performed, the super- Hartree potential, and their sum for a Si atom in a par- cell approximation is good, provided that the supercell allelepiped supercell with side lengths of 2.5, 11, and 11 is large enough (see Fig. 3). In static calculations, then, a.u. respectively in the x, y and z directions. No cut- the use of our cutoff only has the effect of allowing us to off is used here. The ionic potential is roughly behaving eventually use a smaller supercell, what provides clear like 1 in the areanot too closeto the nucleus (where the computational savings. In the case of the Si-chain a r pseudopotential takes over). The total potential, on the full 3D calculation would need of a cell size of 38 a.u. other hand, falls off rapidly to an almost constant value whereasthecutoffcalculationwouldgivethesameresult at around 4 a.u. from the nuclear position, by effect of with a cell size of 19 a.u. Of course, when more delocal- the electron screening. ized states are considered,like higher energy unoccupied Fig. 2 (bottom) shows the results when the cutoff is states,largerdifferencesareobservedwithrespecttothe applied (the radius of the cylinder is R = 5.5 a.u. such supercell calculation. that there is zero interaction between cells). The ionic In Fig. 4 a Na chain with lattice constant 7.5 a.u. is 8 5 B. Static polarisability After the successful analysis of the ground state prop- 0 erties with the cutoff scheme, we have applied the modified Coulomb potential to calculate the static polaris- ) V ability of an infinite chain in the Random Phase Ap- e (F -5 proximation (RPA). As a test case we have considered E - a chain made of hydrogen atoms, two atoms per cell E at a distance of 2 a.u. The lattice parameter was 4.5 -10 a.u. Forthissystemwehavealsocalculatedexcitedstate propertiesinmany-bodyperturbationtheory,inparticu- lar the quasiparticle gap in Hedin’s GW approximation6 and the optical absorptionspectra in the Bethe-Salpeter -15 0 0.1 0.2 0.3 0.4 0.5 framework7,36 (see subsections below). All these calcu- q /G x x lations have been performed in the code SELF.32 The polarisability for the monomer, i.e. a finite system, in Figure 3: Silinear chain in a supercell size of 4.9x19x19 a.u. theRPAapproximationincludinglocalfieldeffectsisde- fined as 1 Ω α= lim χ (q) . (29) −q→0q2 00 4π 5 where χ (q) is the interacting polarisation function GG′ that is solution of the Dyson like equation 4 V) 3 χGG′(q)=χ0GG′(q)+ χ0GG′′(q)v(q+G′′)χG′′G′(q). e (EF 2 XG′′ (30) E- and χ0 is the non interacting polarisation function ob- 1 tained by the Adler-Wiser expression.37 v(q+G) are the Fourier components of the Coulomb interaction. Note 0 that the expression of α in Eq. 29 is also valid for calculations in finite systems, in the supercell approximation, -1 andthe dependence fromthe wave-vectorq is due to the 0 0.1 0.2 0.3 0.4 0.5 representation in reciprocalspace. q /G x x In the top panel of Fig. 5 we compare the values of the calculated polarisability α for different supercell Figure 4: Effect of the cutoff in a Na linear chain in a su- sizes. α is calculated both using the bare Coulomb percell size of 7.5x19x19 a.u. The bands obtained with an v(q+G) = 4π and the modified cutoff potential of ordinarysupercellcalculationwithnocutoff(dashedline)are |q+G|2 compared to the bands obtained applying the 1D cylindrical Eq.(12) (the radius of the cutoff is always set to half cutoff (solid line). As it is explained in the text, only the the inter-chain distance). The lattice constant along the unoccupied levels are affected by thecutoff. chainaxisiskeptfixed. Usingthecutoffthestaticpolar- isability already converges to the asymptotic value with an inter-chaindistance of 25 a.u., while without the cut- offtheconvergenceismuchslower,andtheexactvalueis approximated to the same accuracy for much larger cell consideredinacellof7.5x19x19a.u.,andtheeffectofthe sizes(beyondthecalculationsshowninthetopofFig.5). cutoff on the occupied and unoccupied stated is shown. We must stress that the treatment of the divergences As expected, the occupied states are not affected by the in this case is different with respect to the case of the use of the cutoff, since the density of the system within Hartreeand ionic potential cancellationfor ground-state the cutoff radius is unchanged, and the corresponding calculations (i.e. charge neutrality). In fact, while in bandisthe sameasitis foundforanordinary3Dsuper- the calculation for the Hartree and ionic potential the cell calculation with the same cell size. However there is divergingtermsaresimply droppedbyvirtueoftheneu- a clear effect on the bands corresponding to unoccupied tralising positive background, here the h-dependence in states, and the effect is larger the higher is the energy of Eq.(25)canberemovedonlyfortheheadcomponentby the states. In fact the high energy states, and the states virtueofthevanishinglimitlim χ0 (q)=0,whilefor q→0 00 inthecontinuumaremoredelocalized,andthereforethe the other G = 0 components we have to resort to the x effect of the boundary conditions is more sensible. expression of the finite cylindrical cutoff as in Eq.(25). 9 A finite version of the 1D cutoff has been recently ap- Interchain distance (a.u.) plied to nanotube calculations.30,38 This cutoff was ob- 200 15 20 25 30 35 40 tained by numerically truncating the Coulomb interac- bare Coulomb cut off Coulomb tion along the axis of the nanotube, in addition to the radial truncation. Therefore the effective interaction is limited to a finite cylinder, whose size can be up to a 150 hundred times the unit cell size, depending on the den- ) 3 sity of the k-point sampling along the axis.5 The cutoff u. a. axial length has to be larger than the expected bound ( α exciton length. 100 InthebottompartofFig.5wecomparetheresultsob- tainedwithouranalyticalcutoff(Eq.(12))withitsfinite special case as proposedin Ref. 30. We observe that the valueofthestaticpolarisabilitycalculatedwiththefinite 50 cutoff oscillates aroundan asymptotic value, for increas- 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 ing axial cutoff lengths. The asymptotic value exactly supercell volume (a.u.3) coincideswith the valuethatis obtainedwithourcutoff. Westressthatwealsoresorttothefiniteformofthecut- offonly forthe divergingofcomponents ofthe potential, 155 thuswenote thatthereis aclearnumericaladvantagein 90 180h (a.u.) 270 360 170 using our expression, since the cutoff is analytical for all vtuarlueehsaesxtcoebpetnautmGexri=ca0ll,yaenvdaltuhaetecdofrorerstphoensedipnoginqtusaodnrlay-. 3α (a.u.)115600 In the inset of the bottom part of Fig. 5 it is also shown 150 the convergences of the polarisability obtained with our 3u.) 14010 20k-points 30 40 cutoff with respect to the k-points sampling. The sam- a. ( plingisunidimensionalalongtheaxialdirection. Observe α thatthecalculationusingourcutoffisalreadyconverged for a sampling of 20 k-points. In the upper axis it is also 145 indicated the corresponding maximum allowed value of the finite cutoff length in the axial direction that has been used to calculate the G =0 components. x Finitesizeeffectsturnouttoberelevantalsoformany- 0 100 200 300 400 h (a.u.) bodyperturbationtheorycalculations. Forthesametest system (linear H -chain), in the next two subsections, 2 we consider the performance of our cutoff potential for Figure 5: Top: Polarisability per unit cell of an H2 chain in RPAapproximation as a function of thesupercell volume. the calculation of the quasiparticle energies in GW6 ap- The solid line joins the values obtained with the cutoff po- proximation and in the absorptionspectra in the Bethe- tential, while the dashed lines joins the values obtained with Salpeter framework.7,36 thebareCoulombpotential. Thecutoffradiusis8.0a.u. The inter-chaindistanceisindicatedinthetopaxis. Bottom: Po- larisability of the H2 chain calculated with the finite cutoff 1. Quasiparticles in the GW approximation potential of Ref.30. In abscissa different values of the cutoff lengthalongthechainaxis. Thedashedstraightlineindicates the value obtained with the cutoff of Eq.(12).In the inset we In the GW approximation, the non-local energy- showtheconvergenceofthepolarisability withrespecttothe dependent electronic self-energy Σ plays a role sim- k-pointssamplingalongthechainaxisobtainedwiththecut- ilar to that of the exchange-correlation potential of offofEq.(12). Intheupperaxisitisindicated themaximum DFT. Σ is approximated by the convolution of the allowed lengthhfor eachk-pointsampling usedin thecalcu- one electron Green’s function and the dynamically lation of theGx =0 componentsby Eq.(25). screened Coulomb interaction W. We first calculate the ground state electronic properties using the DFT code ABINIT.33 These calculation are performed in LDA35, and pseudopotentials34 approximation. An energy cutoff of 30 hartree has been used to get converged results. The LDA eigenvalues and eigenfunctions are then been calculated at the first order of perturbation the- used to construct the RPA screened Coulomb interac- ory in Σ V .40 Dividing the self-energy in an ex- xc tion W, and the GW self-energy. The inverse dielec- changeΣ −andacorrelationΣ parts( φDFT ΣφDFT = ptroiclemapatprrioxxǫim−G,1aGt′iohna3s9baenedntchaelcquulaatseipdaurtsiicnlgeethneerpgliaessmhaovne- hinφgDj FreTp|Σrexxs|eφnDitaFtTioin+fhoφrDjtFhTe|Σsecl|cfφ-DienFeTrgi)yh,winjegae|tplta|hneief-owllaoivwes- 10 basis set: Interchain distance (a.u.) 18 25 32 38 7 6 d3q bare Coulomb hnk|Σx(r1,r2)|n′k′i=− (2π)3 v(q+G)× V) 6.5 5 cut off Coulomb ×ρnn1(Xqn1,GBZz)ρ⋆n′n1(qX,GG)fn1k1 (31) point (e 6 40 5 Rc [a.u.]10 15 X 5.5 and e h at t 5 1 d3q p hnk|Σc(r1,r2,ω)|n′k′i= 2 (2π)3( v(q+G′) P Ga 4.5 Xn1BZz GXG′ Q 4 dω′ ρ (q,G)ρ⋆ (q,G′) ǫ−1 (q,ω′) × nn1 n′n1 2π GG′ 3.50 2000 4000 6000 Z 3 supercell volume (a.u. ) f 1 f n1(k−q) n1(k−q) + − ×hω−ω′−ǫLn1D(kA−q)−iδ ω−ω′−ǫLn1D(kA−q)+iδi) Figure 6: Convergence of the GW quasiparticle gap for the (32) H2chainasafunctionofthecellsize,usingthebareCoulomb potential(dashedline)andthecutoffpotential(solid line).In where ρ (q+G) = nkei(q+G)·r1 n k and the in- the inset the behaviour of the GW quasiparticle gap as a tegral innnth1e frequencyhdom| ain in Eq|.(132)1ihas been an- function of the value of the cutoff radius for a supercell with inter-chaindistanceof32a.u. isshown. Theplateauobtained alytically solved considering the dielectric matrix in the plasmon pole mode: ( ǫ−1 (ω) = δ +Ω /(ω2 around a radius of 8 a.u. (i.e. one fourth of the supercell G,G′ G,G′ G,G′ − size) corresponds to the situation in which the radial images ω˜G2,G′)). of chains no longer mutually interact, and the calculation is In order to eliminate the spurious interactionbetween converged. Increasingtheradiusaboveapproximately12a.u. differentsupercells,leavingthebareCoulombinteraction the interaction is back and produces oscillations in the value unchanged along the chain direction, we just introduce of the gap. the expression of Eq.(12) in the construction of Σ and x Σ , and also in the calculation of ǫ−1 . As we did for c GG′ thecalculationofthestaticpolarisability,thedivergences Fig.6. Notice that still at 38 a.u. inter-chain distance appearing in the components (G = 0) cannot be fully the GW gap is underestimated by about 0.5 eV. A sim- x removedand for such components we resortto the finite ilar trend (but with smaller variations) has been found version of the cutoff potential Eq.(25). by Onida et al.28, for a finite system (Sodium Tetramer) In Fig. 6 the convergence of the quasiparticle gap at using the cutoff potential of Eq.(11). Clearly there is a the X point is calculated for different supercell sizes in strong dimensionality dependence of the self-energy cor- the GW approximation. A cutoff radius of 8.0 a.u. has rection. Thenon-monotonicbehaviourversusdimension- beenused. Whenthecutoffpotentialisused,60kpoints ality of the self-energy correction has also been pointed inthe axisdirectionhas beennecessaryto getconverged out in Ref. 41 where the gap-correction was shown to results. In the inset of Fig. 6 we show the behaviour of have a strong component of the surface polarisation. thequasiparticlegapinfunctionofthecutoffradius. We observethatforR >6a.u. aplateauisreached,and,for c R > 12 a.u., a small oscillation appears due to interac- c 2. Exciton binding energy: Bethe-Salpeter equation tionbetweenthetailsofthechargedensityofthesystem withitsimageintheneighbourcell. Differentlyfromthe Startingfromthequasiparticleenergieswehavecalcu- DFT-LDA,calculationforneutralsystems,wherethesu- lated the optical absorption spectra including electron- percell approximation turns out to be good, as we have hole interactions solving the Bethe-Salpeter equation36. discussed above, we can see that the convergence of the Thebasissettodescribetheexcitonstateiscomposedby GW quasiparticle correction turns out to be extremely productstates of the occupiedand unoccupied LDA sin- slow with respect to the size of the supercell and huge gle particle states and the coupled electron-hole excited supercells are needed in order to get converged results. states S = A a† a 0 ,where 0 istheground This is due to the fact that in the GW calculation the | i cvk cvk ck vk| i | i addition of an electron (or a hole) to the system induces state of the system. Acvk is the probability amplitude of finding an exPcited electron in the state (ck) and a hole chargeoscillationintheperiodicimagestoo. Itisimpor- tant to note that the slow convergence is caused by the in (vk), and it satisfies the equation correlationpartoftheself-energy(Eq.(32)),whiletheex- changepartisrapidlyconvergentwithrespecttothecell (ǫQP ǫQP)A + K A =E A ck − vk vck vck,v′c′k′ v′c′k′ S vck size. The use of the cutoff Coulomb potential really im- vck,v′c′k′ X proves drastically the convergence as it is evident from (33)