Accurate and efficient linear scaling DFT calculations with universal applicability

8.2 MB·English

by  Stephan Mohr

Accurate and efficient linear scaling DFT calculations with universal applicability
Stephan Mohr,1 Laura E. Ratcliff,1,2 Luigi Genovese,1,∗ Damien
Caliste,1 Paul Boulanger,1 Stefan Goedecker,3 and Thierry Deutsch1
1Univ. Grenoble Alpes, CEA, INAC-SP2M, F-38000 Grenoble, France
2Argonne Leadership Computing Facility, Argonne National Laboratory, Illinois 60439, USA
3Institut fu¨r Physik, Universita¨t Basel, Klingelbergstr. 82, 4056 Basel, Switzerland
(Dated: January 26, 2015)
Density Functional Theory calculations traditionally suffer from an inherent cubic scaling with
respect to the size of the system, making big calculations extremely expensive. This cubic scaling
can be avoided by the use of so-called linear scaling algorithms, which have been developed during
5
the last few decades. In this way it becomes possible to perform ab-initio calculations for several
1
tensofthousandsofatomsorevenmorewithinareasonabletimeframe. However,eventhoughthe
0
use of linear scaling algorithms is physically well justified, their implementation often introduces
2
some small errors. Consequently most implementations offering such a linear complexity either

n yield only a limited accuracy or, if one wants to go beyond this restriction, require a tedious fine
a tuning of many parameters. In our linear scaling approach within the BigDFT package, we were
J able to overcome this restriction. Using an ansatz based on localized support functions expressed

3 inanunderlyingDaubechieswaveletbasis–whichoffersidealpropertiesforaccuratelinearscaling
2 calculations–weobtainanamazinglyhighaccuracyandauniversalapplicabilitywhilestillkeeping
the possibility of simulating large systems with only a moderate demand of computing resources.

]
i
c
s I. INTRODUCTION and work with other quantities which are strictly local-
l- ized. One such quantity is the density matrix F(r,r(cid:48)),
r whichisrelatedtotheKohn-ShamorbitalsviaF(r,r(cid:48))=
t TheKohn-ShamformalismofDensityFunctionalThe-
m ory (DFT)1,2 has established itself as one of the most (cid:80)ifiψi(r)ψi(r(cid:48)) with fi being the occupation number
t. powerful electronic structure methods due to its good of orbital i. For insulators or metals at finite tempera-
a ture it can be shown that the elements of F(r,r(cid:48)) decay
balance between accuracy and speed and is thus popular
m exponentially with the distance between r and r(cid:48)10–16.
in various fields such as physics, chemistry, biology, and
- material sciences. Kohn-Sham DFT maps the problem Consequently, ifoneneglectselementswhicharebelowa
d given threshold, the number of non-zero elements scales
ofinteractingelectronsontoaproblemofnon-interacting
n
only linearly with respect to the size of the system, thus
o quasi-electrons: given a system containing N electrons,
paving the way towards an algorithm with this reduced
c the 3N-dimensional many-electron wave function is as-
[ sumed to be given by a single Slater determinant built complexity. During the past decades this fact has led to
from N orthonormal single particle orbitals ψ(r) which a number of algorithms which are capable of performing
1
linearscalingcalculations;anoverviewisgiveninRef.17.
v a priori extend over the entire system. The drawback
4 of this approach is that it has an inherent cubic scaling One particular approach is to write the density ma-
8 due to the orthonormality which is imposed on these so- trix in separable form by introducing a set of so-called
8 calledKohn-Shamorbitals. ToorthonormalizeasetofN support functions18,19; an idea that has already been
5 functions it is necessary to calculate the scalar products used in various codes such as ONETEP20, Conquest21,
0 among all of them, which has the complexity O(N2). In CP2K22 and SIESTA23. These support functions can
.
1 additionthecostofcalculatingonescalarproductispro- also be thought of as a localized basis to directly rep-
0 portionaltothesizeoftheunderlyingbasis,m . Since resent the Kohn-Sham orbitals. Even if there exists an
5 basis
both N and m are proportional to the total size of idealsetoflocalizedfunctions–namelythemaximallylo-
1 basis
the system – typically indicated by the number of atoms calizedWannierfunctions24 –thisisofcoursenotknown
:
v n –oneendsupwithacomplexityO(n3 ). Thisbehav- beforehand.
at at
i
X ior is common to all programs using a systematic basis Another important aspect is the choice of the under-
set, be it plane waves3–5, finite elements6 or wavelets7. lying basis which is used to give a numerical representa-
r
a Even when not using a systematic basis set, e.g. in the tion of the support functions. Ideally, as we would like
case of Gaussians8 or atomic orbitals9, the complexity toefficientlydescribelocalizedfunctions,suchabasisset
remains O(n3 ) due to the matrix diagonalization which should at the same time feature orthogonality and com-
at
is required in a straightforward implementation. Due to pact support. This can indeed be offered by Daubechies
this limitation this standard approach is only suitable wavelets25,makingthemanidealbasissetforlinearscal-
for systems containing a few hundred or at most a few ing calculations. Moreover, wavelets are able to yield a
thousand atoms. reasonable precision with only moderate values for the
Inordertoobtainanalgorithmwhichcanscalelinearly numericalgridspacing,thuskeepingthenumberofbasis
with respect to the size of the system, one thus has to functions relatively small for a systematic approach. In
abandontheconceptoftheextendedKohn-Shamorbitals addition the multiresolution properties of wavelets allow
2
for an adaptive resolution and thus to have a finer mesh Here H(r) is the Kohn-Sham Hamiltonian, defined by
only in those regions – close to the atoms – where it is
1
actually required. H(r)=− ∇2+V (ρ(r),r)+V (r), (2)
As for a given system the localized Wannier functions 2 KS PSP
are unknown, there are two possibilities: either use a
where V (ρ(r),r) = V (r)+(cid:82) ρ(r(cid:48)) dr(cid:48) +V (r) is
larger number of support functions and hope that the KS ext |r−r(cid:48)| XC
mostimportantfeaturesoftheWannierfunctionscanbe the Kohn-Sham potential containing the external poten-
captured in this way, or to use a smaller number of sup- tial, the Hartree potential and the exchange-correlation
port functions and optimize them in situ (namely within potential, and VPSP is the pseudopotential representing
apredefinedlocalizationregion)togetascloseaspossible the ions in the system. In BigDFT the pseudopotentials
totheWannierfunctions. Wehaverecentlyimplemented are norm-conserving GTH-HGH27 pseudopotentials and
a linear scaling code based on Daubechies wavelets, dis- their Krack variants28, possibly enhanced by a nonlinear
tributed in the BigDFT package26, where we decided to core correction29.
followthesecondapproach,inthiswayobtainingsupport The linear scaling version of BigDFT is based on a
functions which are optimally adapted to their chemical separable ansatz for the density matrix. Using so-called
environmentandthusyieldaveryhighaccuracy. Conse- support functions φα(r) and a density kernel K the den-
quently we are in this way able to drastically reduce the sity matrix can be written as
sizeofthesubspaceinwhichthedensitymatrixisrepre-
(cid:88)
sented. Together with the aforementioned properties of F(r,r(cid:48))= φα(r)Kαβφβ(r(cid:48)). (3)
Daubechies wavelets this leads to an optimally reduced α,β
number of degrees of freedom and enables us to perform
Thechargedensity, whichenterstheKohn-ShamHamil-
calculations on thousands of atoms while only requiring
tonian,isgivenbythediagonalpartofthedensitymatrix
moderate amounts of memory and CPU time.
and thus reads
In this paper we will present the progress we have re-
cently made for this linear scaling code. Thanks to the (cid:88)
ρ(r)= φ (r)K φ (r). (4)
α αβ β
previously presented ingredients our code offers the abil-
α,β
ity to perform DFT calculations which at the same time
yield an astonishingly high accuracy, are universally ap- Using the orthonormality relations between the support
plicable, exhibit a linear complexity, and require only a functions and their duals, (cid:82) φ (r)φ˜ (r)dr = δ , it fol-
α β αβ
verysmallamountofcomputationalresources. Inpartic- lows that the elements of the density kernel are given
ularitofferstheadvantagethatonedoesnotneedtodoa by
tedious fine tuning of various parameters and in particu-
(cid:90)(cid:90)
lar of the basis, but can straightforwardly get reasonable K = φ˜ (r)F(r,r(cid:48))φ˜ (r(cid:48))drdr(cid:48). (5)
αβ α β
results using a set of default values.
We will show that even with a relatively simple set of
WritinganalogouslytheelementsoftheHamiltonianma-
input parameters one can calculate energies with an ab-
trix as
solute accuracy of the order of 10meV/atom for a large
variety of isolated systems. From a computational point (cid:90)
H = φ (r)H(r)φ (r(cid:48))dr, (6)
of view, as the approach scales linearly with respect to αβ α β
systemsize,itispossibletousetheconceptof“CPUmin-
utes per atom” to evaluate the computational resources it follows that the band structure energy is given by
needed for a full DFT calculation. We will first present
E =tr(KH). (7)
the basic principles of our approach and then demon- BS
stratealloftheaforementionedpropertiesbyvariousex-
The algorithm thus consists of two parts: finding a good
amples.
setofsupportfunctionsφ ,whichgiverisetotheHamil-
α
tonianmatrixH,anddeterminingthedensitykernelKin
thedual(andneverexplicitlycalculated)basisφ˜ (r). In
α
II. THEORY AND ALGORITHM
order to reach linear scaling, the support functions must
be localized and the density kernel sparse. This can be
A. Overview of the algorithm achieved by introducing cutoff radii around the center of
eachsupportfunctionandsettingtozeroallcomponents
Using the formalism based on the density matrix which lie outside of the region; for the density kernel an
F(r,r(cid:48)), the band structure energy, which is the central elementK issettozeroifthecentersofφ andφ are
αβ α β
quantity to be calculated, is given within the framework fartherawaythanthesumofthekernelcutoffsofthere-
of Kohn-Sham DFT by the expression gions α and β – these kernel cutoffs actually correspond
to the cutoffs of the dual support functions (cf. Eq. 5)
(cid:90) (cid:12)
EBS = H(r(cid:48))F(r,r(cid:48))(cid:12)(cid:12)(cid:12) dr(cid:48). (1) acuntdoffarsetihnemgesneelvraels,asbeiettlahregedristchuasnsiothneisnupSpeco.rtIIfuBn.ctTiohne
r=r(cid:48)
3
B. Optimization of the support functions and
density kernel
The support functions are optimized by minimizing a
target function. Ideally this target function should lead
tostronglylocalizedsupportfunctionswhichatthesame
time yield a very high accuracy. For the latter property
the correct quantity to be minimized is the band struc-
ture energy of Eq. (7). For the first property the choice
is not unique; one possibility would be
(cid:88)
Ω= (cid:104)φ |H |φ (cid:105), (8)
α α α
α
where H is the Kohn-Sham Hamiltonian including a
α
confinement, i.e. H =H+c (r−R )4, with R being
α α α α
thecenterofthelocalizationregion. Inordertocombine
the two properties we define our target function as
(cid:88) (cid:88)
Ω= K (cid:104)φ |H |φ (cid:105)+ K (cid:104)φ |H|φ (cid:105), (9)
αα α α α αβ α β
α β(cid:54)=α
where the prefactor for the confinement is smoothly re-
duced during the run, as explained in more detail in
Ref. 26. In this way we have a strong confinement in
the beginning, leading to a decent localization, and still
obtainahighprecisionaswecorrectlyminimizetheband
structureenergyintheend. Sincethedecreaseofthecon-
finement is done automatically – taking into account the
FIG. 1. Flowchart of the linear scaling algorithm. It consists properties of the system –, this procedure is universally
ofadoubleloopstructure: intheinnerfirstloopthesupport applicable without the need for any fine tuning.
functions are optimized with a fixed density kernel, in the Incontrasttootherlinearscalingcodesemployingthe
secondinnerloopthedensitykernelisoptimizedusingafixed same support function approach20–23 we decided to im-
setofsupportfunctions. Thesetwoloopsaretheniteratedin pose an additional constraint and to keep the support
an outer loop until overall convergence is reached. functions approximately orthogonal. In other terms, we
optimize φ such that the overlap matrix
α
S =(cid:104)φ |φ (cid:105) (10)
αβ α β
is close to the identity matrix. As the kernel cutoff is
related to the extension of the dual support functions
support functions are again expanded in terms of an un- (see Eq.(5)), the sparsity of K is governed by the cut-
derlying wavelet basis, as is explained in more details in off implicitly used for these dual functions. By choosing
Ref. 26. thesupportfunctionstobeorthogonalthedualandnon-
dualentitiesare identicaland thecutoffsfor thesupport
functionsandthedensitykernelmatrixelementscancon-
A flowchart of the algorithm is shown in Fig. 1. It sequently be chosen along the same lines, leading to a
consistsofanesteddoubleloopstructure,withoneouter degree of sparsity for the kernel which is comparable to
loopandtwoinnerloops. Inthefirstinnerloop,thesup- that of the overlap matrix itself.
portfunctionsareoptimizedusingafixedkernelandthus In addition, quasi-orthogonal support functions lead
also using a fixed potential and Hamiltonian. Therefore to an overlap matrix which is close to the identity and
it is not advisable to perform too many iterations in this whose inverse can thus be cheaply calculated using poly-
loop,butrathertoexiteveniftheconvergencethreshold nomial expansions. This property is in particular impor-
has not been reached yet. In the second inner loop, the tant for the Fermi Operator Expansion30,31 (FOE) pro-
density kernel is optimized using a fixed set of support cedure that we use to determine the density kernel. The
functions. In contrast to the first inner loop, this second efficiency of this approach relies heavily on the sparsity
one is done self-consistently, i.e. the Hamiltonian is up- of the matrix S−1/2HS−1/2. Even if H is very sparse
datedineachiteration. Thesetwoloopsaretheniterated thanks to the cutoff radii of the support functions, this
in the outer loop until overall convergence is reached. propertyismuchlesspronouncedforS−1/2 inthecaseof
4
non-orthogonal support functions. Thus, if one wants to
keep a high degree of sparsity also for the matrix prod-
uct, it is indispensable to work with a set of orthogonal
support functions.
Inspiteofall thesebenefits, itmustbestatedthator-
thogonality and locality are in general two contradicting
properties; consequently a strict enforcement of the or-
thogonality might lead to convergence problems. There-
fore we only perform an explicit orthogonalization in the
verybeginning; inthefollowingtheorthogonalityisonly
approximatelyconservedbytheuseofaLagrangemulti-
plier in the support function gradient. This is enough to
keep the overlap matrix diagonally dominant – the off-
diagonal elements being typically of the order of 0.1 –
and thus to maintain the aforementioned benefits.
For the determination of the density kernel the code
offers several possibilities, each one having its particular
strengthsandthusareasofapplication. Inthispaperwe
always used the FOE method, as it is the only method
thatallowstoperformcalculationsforverylargesystems.
For more details we again refer to Ref. 26.
III. ASSESSMENT OF THE ACCURACY
FIG. 2. Illustration of the six systems used for the compari-
son of the linear and cubic energies and forces: Vitamin B ,
12
In order to assess the precision of our code, it has to Chlorophyll, Si H , C , (H O) , DNA.
87 76 60 2 100
be compared with a reference. We use as such the cubic
version of the code7, which has been shown to yield very
accurate results29. Thus, by demonstrating that the lin- seem to be fairly universal, meaning that any other sys-
ear version of BigDFT is able to reproduce the results of tem can be readily calculated without the need for any
thecubic scalingversion, wecan be suretohavereached adjustments.
a very high level of precision. All calculations were done TheresultsaresummarizedinTab.I.Ascanbeseen,it
usingfreeboundaryconditionsandthePBE32 functional ispossibletogetwiththeseinputparametersafairlycon-
unless otherwise stated. stant energy offset between the linear and cubic versions
of about 10meV/atom. The difference of the force norm
is typically of the order of 10−3hartree/bohr. These are
A. Comparison of energies and forces morethanreasonablevaluesinviewofthefactthatthey
were obtained without any fine tuning. Last but not
least, it is worth noting that support functions are op-
One of the most important features of our code is the
timized such that Pulay forces33 are absent in our ap-
ability to yield a high level of precision without the need
proach26, and the evaluation of the forces is therefore
for doing a tedious fine tuning of the parameters. As an
straightforward.
illustration of this property, we show the energies and
forces calculated with both the linear and cubic scaling
version for six systems with rather different electronic
structures; they are depicted in Fig. 2. The chosen pa- B. Calculations of energy differences
rameterswereverysimilarforallsystems,withthemain
difference being the various cutoff radii used for the sup- Absolute energies are always somewhat arbitrary as
portfunctionregionsofthedifferentatomkinds. Tokeep they are only defined up to an additive constant. There-
it simple, the choice was just based on the row in which fore energy differences are more meaningful as this am-
the element appears in the periodic table: 5.5bohr for biguity vanishes. For the linear scaling version there is
elementsofthefirstrow, 6.3bohrforthoseofthesecond in addition the benefit that the offset caused by the fi-
row, and 7.0bohr for those of the third row. The kernel nite cutoff radii cancels to a large degree and it is thus
cutoff was set to 9.0bohr; an automatic adjustment for possible to get very close to the exact result. The value
technical reasons was done by the code if the value was of 10meV/atom mentioned in the previous section is the
smaller than the cutoff radius of the support function error in the absolute energy and can thus be considered
plus8timesthegridspacing,seeRef.26fordetails. The as an upper bound; in fact it is however possible to ob-
number of support functions per atom was 1 for H, 9 for tainamuchhigheraccuracythanonemightthinkatfirst
Si, and 4 for all the other elements. These parameters sight.
5
cubic linear difference
energy force norm energy force norm energy force norm
hartree hartree/bohr hartree hartree/bohr meV/atom hartree/bohr
Vitamin B −926.78 2.13·10−3 −926.71 1.93·10−2 11.43 1.72·10−2
12
Chlorophyll −476.70 3.05·10−3 −476.64 1.24·10−2 12.33 9.38·10−3
Si H −386.79 7.80·10−4 −386.69 1.01·10−2 17.20 9.27·10−3
87 76
C −341.06 2.69·10−4 −341.02 7.36·10−3 17.23 7.09·10−3
60
(H O) −1722.99 5.23·10−1 −1722.87 5.26·10−1 10.89 2.80·10−3
2 100
DNA −4483.12 5.60·10−1 −4482.84 5.63·10−1 10.69 2.40·10−3
TABLE I. Comparison of the energies and forces for different systems, calculated using the linear and cubic versions. Linear
energies have an offset of about 10meV/atom with respect to the cubic version, and the deviations of the linear forces from
their cubic counterparts is of the order of 10−3hartree/bohr. The first four systems were close to their geometrical ground
state, whereas the last two were unrelaxed. All systems were calculated with similar parameters.
pure impurity difference
eV eV eV
cubic -53179.7148 -53251.5162 71.8014
linear -53168.4248 -53240.3371 71.9123
difference 11.2900 11.1791 0.1109
relative difference - - 0.15%
TABLE II. Energy differences between the pure and the im-
pure Si wire, for both the linear and the cubic version. The
energyoffsetbetweenthelinearandthecubicversioncancels
toalargedegreeandthusyieldsaveryaccurateresultforthe
energy difference.
FIG.3. Thetwosiliconwires(consistingof660atoms)which
Asanexample,wecalculatedtheenergydifferencebe-
are used to benchmark the calculation of energy differences.
tweenahydrogen-passivatedsiliconwireinitspurestate
Whereastheoneontheleftsideisapurewire,thereisadefect
and another one which contains a defect. As a defect we
intheoneontherightsidewithonesiliconatomhavingbeen
choseasubstitutionalatom; oneofthesiliconatomswas replacedbyaphosphorusatom. Theredcirclehighlightsthis
replaced by phosphorus. These systems have been used substitutional atom.
asacasestudytodeterminethebindingenergyofimpu-
ritiesinsemiconductorsusingchargedDFTcalculations,
seeRef.34. ThetwoconfigurationsaredepictedinFig.3.
C. Geometry optimizations
Apart from the fact that silicon is a rather delicate sys-
tem, this defect is even more challenging as it requires
the comparison of two systems with a different number A good test to check at the same time the accuracy
of electrons. The calculations were performed with the of the energies and the forces is to perform a short ge-
analogous parameters as those used for the benchmarks ometry optimization for both the linear and the cubic
in Sec. IIIA, again demonstrating their universal char- version, starting from the same non-equilibrated struc-
acter. From the results in Tab. II it becomes clear that ture. If the forces calculated by the linear version are
energy differences can indeed be calculated with a very accurate enough, this should lead to identical trajecto-
high precision. Whereas the offset between the linear ries and thus to a parallel evolution of the energies and
andcubicversionintheabsoluteenergyisabout11.2eV the forces. As an example we took an alkane consisting
– which is about 17meV/atom, in agreement with the of 302 atoms. The results are shown in Fig. 4. As can
results of Sec. IIIA –, the discrepancy between the two be seen, the offset in the energy between the linear and
version becomes as little as 0.11eV for the energy differ- the cubic version remains more or less constant through-
ence between the pure and the impure wire. Comparing out the entire geometry optimization, and is again of
this mismatch with the correct value of 71.8eV one gets the order of 10meV/atom. In addition, the forces are
a relative error of only 0.15%. pretty much identical, with variations of the order of
6
-689.10
V) 11.4 E -E
e linear cubic
gy (hartree) ---666888999...432000 energy (m 1111..02 0 40 ite 8ra0tion s120 160 lcinuebaicr
er
n
e -689.50
-689.60 FIG. 5. The 12 energetically lowest configurations for the
0 20 40 60 80 100 120 140 160 fullerenes B C .
12 48
iterations
100
hartree/bohr) 10-1 force norm(Ha/bohr)1100--43 0 20 |4F0lin 6ea0r- F80c u1b0ic0| 120 140 160 vceFuneigrbes.riicgo5env.teiiScrssauinlcolohyntbwlaoobewrlopeesnrtte-ocsseartnrerbtupoctrnhotudefruueecnlsle3ee5rrtgehoneefteicsBcoa1or2rrreeCdce4atr8ir,vneedgsrueyfplotirdscettolheifdcetah1tin2ee
orm ( 10-2 iterations system; previous studies36,37 have even led to different
e n conclusions. In order to successfully describe the ener-
forc 10-3 lcinuebaicr gtoetcicatocrhdeenrienrggyofditffheerseenccoensfi–giu.er.atsiiognnsalist–isoifntdhiespoerndsearbolef
0 20 40 60 80 100 120 140 160
1meV/atom, which is possible with the cubic version of
iterations
BigDFT used in Ref. 35. The energy spectra for these
FIG.4. Evaluationoftheenergy(upperpanel)andtheforces 12 configurations are shown in Fig. 6, for three different
(lowerpanel)asafunctionofthenumberofiterationsinthe setups: the cubic version using the PBE functional, the
geometry optimization. For the energy, the offset remains linear version using as well the PBE functional, and the
rather constant throughout the entire run, leading to two cubicversionusingtheLDA38functional. Ascanbeseen,
parallel curves; as can be seen from the inset, the offset only there are – even if the main features of the system are
varieswithinlessthan1mev/atom. Thecurvesfortheforces
wellcaptured–notabledifferencesbetweenthecubicand
are more or less superposed, with a difference of the order of
the linear version: whereas the ground state is the same,
10−4hartree/bohr. Thelargeplateauinthemiddlewherethe
there has been some reordering of the excited states; in
force remains constant is due to the fact that an adjustment
particular the energetic gap between the ground state
of the bond lengths can only start at the two endpoints of
thealkaneandthenslowlypropagatelikeawavetowardsthe and the first local minimum is wrong. This shows that
center of the molecule. for this system the linear scaling approach using finite
cutoff radii is not appropriate. However the third panel
demonstrates that these boron-carbon fullerenes are in-
only 10−4hartree/bohr. Together this shows that the deed an extremely delicate system. As can be seen, the
two trajectories performed by the linear and the cubic energeticorderingisalsoconsiderablymodifiedifthecu-
versions are identical, demonstrating the high precision bic version is used with another functional: using LDA
of the forces calculated by the linear version. Again it instead of PBE, the energy levels of the ground state
must be stressed that these results were obtained with andthefirstexcitedoneremainidentical, butthehigher
thesamestandardsetofparameterswhichhadbeenused levelsarecompletelyjumbled. Duetothisveryhighsen-
in the previous examples. sitivity it is thus not surprising that the linear version
is not able to satisfactorily reproduce the results of the
cubicversion, asthenoiseintroducedbythefinitecutoff
D. Energetic ordering radii is higher than the signal we look for.
The statement that our approach is universally appli-
cabledoesnotmeanthatitcanblindlybeappliedtoany IV. SCALING WITH SYSTEM SIZE
system. There are situations where the “noise” intro-
ducedbythefinitecutoffradiiislargerthanthe“signal” Thegoalofthisworkwastoobtainacodewhichatthe
oneislookingfor;inthiscasetheresultsofourapproach sametimeyieldsaveryhighprecisionandscaleslinearly.
– as of any method relying on such a truncation – must The first property has already been shown in detail in
be interpreted with care. By universally applicable we Sec. III; now it remains to demonstrate the second prop-
rather mean that we can straightforwardly – i.e. without erty. To this end we performed single point calculations
the need of doing a tedious fine tuning of the input pa- forrandomlygeneratedwaterdropletsofvarioussize,for
rameters and the basis set – apply our method to those boththelinearandthecubicversion. Therunswereper-
cases where the “signal to noise ratio” is large enough formed in parallel, using in total 6400 cores. The results
even when finite cutoff radii are used. for the total runtime and the memory peak are shown
As an example for a system where the linear scaling in Fig. 7. As can be seen both quantities clearly exhibit
7
5000
total runtime linear
0.7 0.7 0.7 4500 linear extrapolation, refetroetnacl eru 3n0ti0m0e a ctoumbics
4000
3500
0.6 0.6 0.6 ds) 3000
n
o
ec 2500
s
me ( 2000
0.5 0.5 0.5 ti
1500
1000
energy(eV)0.4 energy(eV)0.4 energy(eV)0.4 50 00 2000 4000 60n0u0mbe r8 o0f0 a0tom 1s0000 12000 14000
0.3 0.3 0.3 30
memory peak linear
memory peak cubic
linear extrapolation, reference 3000 atoms
25
0.2 0.2 0.2
B) 20
G
k (
a
0.1 0.1 0.1 pe 15
y
or
m
me 10
0 0 0
5
cubic, PBE linear, PBE cubic, LDA
0
2000 4000 6000 8000 10000 12000 14000
FIG. 6. Energetic ordering of the lowest configurations of
number of atoms
B C . From left to right: cubic version with PBE, linear
12 48
version with PBE, cubic version with LDA. All energy levels
FIG.7. Totalruntimeforonesinglepointcalculation(upper
havebeenshiftedsuchthatthegroundstateisalwaysatzero.
panel)andmemorypeakperMPItask(lowerpanel)forwater
The colors correspond always to the same level, demonstrat-
droplets of various size, ranging from 600 to 15000 atoms.
ing that the energetic ordering is not always the same.
The linear scaling version indeed exhibits linear scaling, as
indicated by the linear extrapolations. The small deviations
aremainlycausedbyaslightlydifferentnumberofiterations
a strict linear scaling. The cubic version, on the other required to reach convergence.
hand, shows a much steeper increase and does not per-
mit calculations beyond about 2000 atoms. Due to the
28
spherical geometry the degree of sparsity of the matrices DNA
is rather low and it is therefore more difficult for the lin- 26 water
earscalingversiontoexploitthisproperty. Nevertheless, 24
even for the smallest droplet the linear scaling version is m 22
o
faster than the cubic one. By extrapolating the curves er at 20
fortheruntimethecrossoverpointforthissystemcanbe s p 18
estimated to be at about 200 atoms, which is in a range ute 16
n
still easily accessible by the cubic version of BigDFT for mi
U 14
production runs (see e.g. Ref. 39). P
C 12
Moreover it becomes clear that the computational re-
10
sources consumed by the linear version are only moder-
8
ate. For a complete single point calculation for 10000
6
atoms only a few thousand CPU hours are required, 5000 10000 15000 20000 25000
which is a rather small amount considering the capac- number of atoms
ities of current supercomputers. This demonstrates that
we were able to drastically reduce the number of degrees FIG. 8. CPU minutes per atom needed for a fully self-
consistentcalculationforDNAfragmentsandwaterdroplets
of freedom while keeping a very high level of accuracy.
ofdifferentsizes. Asexpected,thisquantityremainsconstant
Foratrulylinearscalingregimeitisalsointerestingto
above a certain size, with the prefactor mainly influenced by
look at the ”CPU minutes per atom” which are needed
the geometry of the system to be calculated.
to perform a fully self-consistent calculation. This quan-
8
tity is shown in Fig.8 for DNA fragments and the water 22 110
CPU minutes per atom (left axis)
droplets used for the aforementioned benchmark. As is 20 total walltime (right axis) 100
to be expected for a linear scaling code, this quantity 18 90
remains constant above a certain critical size where the
16 80
linear scaling parts of the code start to dominate. More-
oiwnvafletureeronndcereoscpatlhneetsepatrsheialfyancstfeooerr,htrhoeewsuDltthNineAggiefnroamhgimegtherenyrtsov.faltuheessfyosrttehme PU minutes 111024 567000 minutes
C
8
40
6
30
V. PARALLELIZATION
4
20
2
A. Parallel efficiency 10
4000 8000 12000 16000 20000 24000
total number of cores
Even if the linear scaling version of BigDFT requires
30
only moderate resources, it is indispensable to have an memory peak MPI master
efficient parallelization scheme in order to keep the run- ideal
times short and thus to have the possibility of perform- 25
ingadvancedcalculationssuchasgeometryoptimizations
within a reasonable time frame. To this end we paral- B) 20
G
lelized our code with a mixed distributed/shared mem- k (
a
ory parallelization scheme using MPI and OpenMP. It pe 15
y
is worth noting that the shared memory parallelization or
m
is not just an additional speedup. Reducing the num- me 10
ber of MPI tasks and in turn increasing the number of
OpenMP threads helps to improve the MPI load balanc- 5
ing and to reduce communication overhead; moreover it
can substantially reduce the memory peak per compute 0
node and thus help in situations where this resource is 4000 8000 12000 16000 20000 24000
total number of cores
critical.
Reachinganefficientparallelizationforalinearscaling
FIG. 9. CPU minutes per atom and total walltime (upper
code is not easy. First of all we have a small number of
panel) and memory peak (lower panel) as a function of the
degrees of freedom and thus little workload than can be
totalnumberofcoresforaDNAfragmentconsistingof14300
shared among the cores. Moreover there are also load atoms. Therunswereallperformedusing8OpenMPthreads
balancing problems which can arise due to the different for each MPI task, thus only the latter number was varied.
sizes and surroundings of the localization regions. Nev-
ertheless we were able to reach a degree of parallelism
which is more than sufficient to efficiently perform large
efficient load balancing in such a situation, resulting in
calculations. As an illustration we show the scaling for a
many cores being idle most of the time. Consequently
fully self-consistent calculation of a DNA fragment con-
this degradation in the parallel speedup will be shifted
sistingof14300atoms. ThenumberofMPItasksranged
towards larger values if one increases the system size.
from 160 to 3200, each one being split into 8 OpenMP
threads; thus in total we have a range of 1280 to 25600 Neverthelessitisworthnotingthatevenforthelargest
cores. The results for the runtime and the memory peak number of cores there is still some speedup, as can be
are shown in Fig. 9. For the runtime we show the CPU seenfromthetotalwalltimedecreasingsteadily. Evenin
minutes per atom as well as the total walltime. this range, the CPU time per atom remains of the same
As can be seen from the CPU minutes per atoms re- order as the one needed at the crossover point with the
mainingalmostconstantthespeedupisverygooduptoa cubic code. As calculations with the cubic BigDFT code
few thousand cores. On the other hand this quantity in- arealreadyaccessibleinthisrange(aspointedoutinthe
creases for very large number of cores, revealing that the previous section), production runs of very large systems
speedup becomes only moderate in this range. However might become feasible by linearly scaling the computing
this does not mean that the code is poorly parallelized. resources needed at the crossover point. The memory
The reason is rather that, despite of the considerable parallelization, shown in the lower panel, does not suffer
size of the system, the number of degrees of freedom is from any degradation, and we come close to a perfect
so small that using too many cores simply results in a scaling. In summary these results show that our code
verysmallnumberofoperationstobeexecutedbyasin- has an excellent level of parallelism as long as one keeps
gle core and thus to a poor ratio of computation and a good balance between the size of the system and the
communication. Moreoveritismoredifficulttoreachan computational resources that are utilized.
9
B. Parallelization of the sparse matrices
Some strategies to obtain an efficient parallelization
have already been outlined in Ref. 26. Here we present
a new concept which we developed to reach an efficient
handling of the sparse matrices.
The problem of the transposed approach26 for calcu- FIG.10. Illustrationoftheconceptofthematrixtaskgroups
latingtheoverlapmatrix(andsimilarlytheHamiltonian foratoysystem. Thereare8MPItaskswhichareregrouped
matrix) is that it requires a global reduction operation in3taskgroups. Eachtaskgroupisinitiallydefinedsuchasto
at the end. This can pose a bottleneck both from the form a matrix subset, and subsequently extended to include
viewpoint of runtime and memory. Moreover this global all the processes which share data of the initial subset. Care
is taken in not including each MPI task in more than two
reduction is wasteful due to the locality which is inher-
different taskgroups.
ent to the linear scaling approach. Each MPI task only
needs to know a small portion of the entire matrices and
a global reduction is thus not necessary. To circumvent overlapcalculation gather/compress
this issue, the MPI tasks are regrouped in taskgroups, seconds seconds
sharing a portion of the entire matrix – obviously this withouttaskgroups 88.9 266.9
portion is chosen such that it contains those parts of the withtaskgroups 25.2 16.4
global matrix which are actually needed by each indi-
TABLE III. Speedups offered by the matrix taskgroups for
vidual task. In addition the taskgroups are defined such
two specific operations.
that each MPI task belongs to at most two taskgroups.
This has two advantages: first of all each MPI task only
needs to hold a copy of a part of the global matrix, thus
VI. CONCLUSION
reducing the memory requirements, and secondly the re-
duction only needs to be performed within a taskgroup.
Thanks to the use of non-blocking collective MPI rou- Using a set of strictly localized and quasi-orthogonal
tines a task can participate in two reduction operations support functions which are adapted in-situ to their
at the same time without the risk of serializing the code. chemical environment we were able to develop a code
which can perform accurate DFT calculations with a
strictlinearscalingwithrespecttothesizeofthesystem.
The concept is visualized in Fig. 10 for a toy system Thanks to the compact support property of the underly-
consisting of 8 MPI tasks. The matrix subparts which ingDaubechieswavelets,thisreductionofthecomplexity
areneededbyeachMPItaskareindicatedbyrectangles. with respect to the traditional cubic approach does not
The MPI tasks are initially regrouped in 3 taskgroups: come at the cost of a loss of precision. Indeed, we are
taskgroupIisresponsibleforthepartneededbythetasks able to get an excellent level of accuracy with a set of
1-2, taskgroup II for those needed by the tasks 3-6, and standard input parameters and can thus get rid of the
taskgroup III for those needed by the tasks 7-8. The tediousfinetuningwhichisoftenneededforotherO(N)
taskgroups are then enlarged in order to guarantee that approaches. We believe this would considerably simplify
the reduction is done correctly. Therefore taskgroup I the usage of this code by end-users, as this fine-tuning
includes the tasks 1-4, taskgroups II includes the tasks usually restricts the usage O(N) approaches to a com-
2-8, and taskgroup III includes the tasks 5-8. munity of specialists.
Obviously the toy system is too small to get a real Moreover we were able to drastically cut down the
benefit from the taskgroups. However for large systems degrees of freedom needed for an accurate calculation,
there can be a substantial gain. As an example we show thus reducing the prefactor of the scaling law and con-
in Tab. III the timings obtained with and without the sequently considerably diminishing the amount of com-
matrix taskgroups for two specific operations. As can be putational resources required. Even for systems con-
seen,thecalculationoftheoverlapandHamiltonianma- taining 10000 atoms a complete single point calculation
trix (both computation and communication) can be ac- can be done using only of the order of one thousand
celerated by more than a factor of 3, demonstrating that CPU hours. Furthermore the linear scaling behavior of
actuallymostofthetimewasspentinthereductionand the program motivates the notion of “CPU minutes per
notinthecomputationitself. Theotheroperation,which atom”, making it easy to estimate the computational
is required at the end of the FOE procedure to build up resources needed for the simulation of a given system
thedensitykerneloutofthepartialresultscalculatedby of any size. Together with the efficient parallelization
each task and is almost entirely communication based, scheme this low consumption of computational resources
can even benefit from a speedup of 16. also paves the way towards more expensivetasks suchas
10
geometry optimizations. ing under Contract No. ANR-2010-COSI-005-01 NEW-
Acknowledgments. S. M. acknowledges a postdoc- CASTLE. Computing time has been provided by the
toral fellowship from the CEA Nanosciences horizontal ”Curie” national GENCI-IDRIS supercomputing center
aid programme (project BigPOL). T.D and L. G. are in- under contract No. i2014096905. Computing time of the
debtedtotheFrenchNationalResearchAgencyforfund- CSCS Swiss project s499 is also acknowledged.
∗ [email protected] PhysRevLett.86.5341
1 P.HohenbergandW.Kohn,Phys.Rev.136,B864(1964) 17 S.Goedecker,Rev.Mod.Phys.71,1085(Jul1999),http:
2 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965) //link.aps.org/doi/10.1103/RevModPhys.71.1085
3 X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, 18 W. Hierse and E. B. Stechel, Phys. Rev. B 50,
F.Bottin,P.Boulanger,F.Bruneval,D.Caliste,R.Cara- 17811 (Dec 1994), http://link.aps.org/doi/10.1103/
cas, M. Coˆt´e, T. Deutsch, L. Genovese, P. Ghosez, M. Gi- PhysRevB.50.17811
antomassi, S. Goedecker, D. Hamann, P. Hermet, F. Jol- 19 E. Herna´ndez and M. J. Gillan, Phys. Rev. B 51,
let, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, 10157 (Apr 1995), http://link.aps.org/doi/10.1103/
M.Oliveira,G.Onida,Y.Pouillon,T.Rangel,G.-M.Rig- PhysRevB.51.10157
nanese,D.Sangalli,R.Shaltaf,M.Torrent,M.Verstraete, 20 C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C.
G.Zerah,andJ.Zwanziger,Comput.Phys.Commun.180, Payne, J. Chem. Phys. 122, 084119 (2005)
2582 (2009) 21 D.R.BowlerandT.Miyazaki,J.Phys.: Condens.Matter
4 P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, 22, 074207 (2010)
C.Cavazzoni,D.Ceresoli,G.L.Chiarotti,M.Cococcioni, 22 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello,
I.Dabo,A.D.Corso,S.deGironcoli,S.Fabris,G.Fratesi, T. Chassaing, and J. Hutter, Comput. Phys. Commun.
R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, 167, 103 (2005)
M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, 23 J.M.Soler,E.Artacho,J.D.Gale,A.Garc´ıa,J.Junquera,
R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, P. Ordej´on, and D. Sa´nchez-Portal, J. Phys.: Condens.
C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, Matter 14, 2745 (2002)
A.Smogunov,P.Umari,andR.M.Wentzcovitch,J.Phys.: 24 N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847
Condens. Matter 21, 395502 (2009) (1997)
5 M.D.Segall,P.J.D.Lindan,M.J.Probert,C.J.Pickard, 25 I. Daubechies et al., Ten lectures on wavelets, Vol. 61
P.J.Hasnip,S.J.Clark,andM.C.Payne,J.Phys.: Con- (SIAM, 1992)
dens. Matter 14, 2717 (2002) 26 S. Mohr, L. E. Ratcliff, P. Boulanger, L. Genovese,
6 J. E. Pask and P. A. Sterne, Modell. Simul. Mater. Sci. D. Caliste, T. Deutsch, and S. Goedecker, The
Eng. 13, R71 (2005) Journal of Chemical Physics 140, 204110 (2014),
7 L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, http://scitation.aip.org/content/aip/journal/jcp/
S. A. Ghasemi, A. Willand, D. Caliste, O. Zilberberg, 140/20/10.1063/1.4871876
M.Rayson,A.Bergman,andR.Schneider,TheJournalof 27 C.Hartwigsen,S.Goedecker,andJ.Hutter,Phys.Rev.B
Chemical Physics 129, 014109 (2008), http://link.aip. 58, 3641 (1998)
org/link/?JCP/129/014109/1 28 M. Krack, Theor. Chem. Acc. 114, 145 (2005)
8 M. Valiev, E. Bylaska, N. Govind, K. Kowalski, 29 A. Willand, Y. O. Kvashnin, L. Genovese, A. Va´zquez-
T.Straatsma,H.V.Dam,D.Wang,J.Nieplocha,E.Apra, Mayagoitia, A. K. Deb, A. Sadeghi, T. Deutsch, and
T. Windus, and W. de Jong, Comput. Phys. Commun. S. Goedecker, J. Chem. Phys. 138, 104109 (2013)
181, 1477 (2010) 30 S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122
9 V.Blum,R.Gehrke,F.Hanke,P.Havu,V.Havu,X.Ren, (1994)
K.Reuter,andM.Scheffler,Comput.Phys.Commun.180, 31 S. Goedecker and M. Teter, Phys. Rev. B 51, 9455 (1995)
2175 (2009) 32 J.P.Perdew,K.Burke,andM.Ernzerhof,Phys.Rev.Lett.
10 J.D.Cloizeaux,Phys.Rev.135,A685(Aug1964),http: 77, 3865 (1996), ISSN 0031-9007
//link.aps.org/doi/10.1103/PhysRev.135.A685 33 P. Pulay, Molecular Physics 17, 197 (1969),
11 J.D.Cloizeaux,Phys.Rev.135,A698(Aug1964),http: http://www.tandfonline.com/doi/pdf/10.1080/00268976900100941,
//link.aps.org/doi/10.1103/PhysRev.135.A698 http://www.tandfonline.com/doi/abs/10.1080/
12 W.Kohn,Phys.Rev.115,809(Aug1959),http://link. 00268976900100941
aps.org/doi/10.1103/PhysRev.115.809 34 Y. M. Niquet, L. Genovese, C. Delerue, and T. Deutsch,
13 R. Baer and M. Head-Gordon, Phys. Rev. Lett. 79, Phys. Rev. B 81, 161301 (Apr 2010), http://link.aps.
3962 (Nov 1997), http://link.aps.org/doi/10.1103/ org/doi/10.1103/PhysRevB.81.161301
PhysRevLett.79.3962 35 S. Mohr, P. Pochet, M. Amsler, B. Schaefer, A. Sadeghi,
14 S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett. 82, L. Genovese, and S. Goedecker, Phys. Rev. B 89,
2127 (Mar 1999), http://link.aps.org/doi/10.1103/ 041404 (Jan 2014), http://link.aps.org/doi/10.1103/
PhysRevLett.82.2127 PhysRevB.89.041404
15 S. Goedecker, Phys. Rev. B 58, 3501 (Aug 1998), http: 36 I. Garg, H. Sharma, K. Dharamvir, and V. Jindal, Jour-
//link.aps.org/doi/10.1103/PhysRevB.58.3501 nal of Computational and Theoretical Nanoscience 8,
16 L. He and D. Vanderbilt, Phys. Rev. Lett. 86, 642 (2011), http://www.ingentaconnect.com/content/
5341 (Jun 2001), http://link.aps.org/doi/10.1103/ asp/jctn/2011/00000008/00000004/art00018