History of the DP code
(from my, VO, commentaries)

The DP code was born as my PhD thesis work around 1997 in Rome.

The thesis was dealing with beyond-LDA non-local approximations (NLDA), which use a different expansion with respect to GGA. The target was to improve on the calculated value of the dielectric constant in LDA or GGA. But I decided since the beginning that it was really the case to keep an eye also on the full dielectric function at ω ≠ 0. Indeed, at that time the first many-body Bethe-Salpeter equation optical spectra, in good agreement with the experiment, were appearing in Palaiseau. And, in my heart, the secret hope was that such non-local approximations could get at the same result, but using a simpler density-functional based theory. I have also to say that at that time I had no clear ideas about the differences among DFT, DFPT, DPFT and TDDFT; and also between ground state observables, like the dielectric constant, and excited state ones, like optical absorption. This confusion was removed only later by my supervisor Rodolfo Del Sole who, looking at some nice spectra I had produced, exclaimed: "But this is TDDFT!". And only at that time I realized that I had written a TDDFT code.

Before to start to implement such code, I begun to evaluate the possible options I had in front. One possibility was to numerically calculate the polarizability for q→0, that is, by using sets of Kohn-Sham wavefunctions calculated on different k-points grids, shifted by the small chosen q. Although my supervisors warmly advised such option as safer in the implementation and the results, I didn't like it at all. I indeed realized that such an option would have obliged me to work as a fool on many wavefunction files, with always the risk to confuse them, do errors, etc. My purpose was to have a code where you turn a key and it produces the result. Have the maximum possible automatization. Not to be slave of a code, of informatics. For this reason I decided to consider a more difficult possibility: to calculate analitically the limit q→0 for the polarizability. The more difficult term in such limit results from a commutator with the non-local part of the pseudopotential. But this solution, once implemented, allowed to work only on one wavefunctions file. And it was corresponding much more to my idea of the informatics at the service of physicists, and not the contrary. And this indeed has always inspired the philosophy behind this code.

Valerio Olevano

I started then to implement the main routines of the code. For the k and r-points symmetries, I was using a library of functions coming from the GW Godby's code, that I was also mantaining and evolving; for the FFT, a library by ICTP Nobile and Roberto, replaced very early by the Goedecker library; for the matrix inversions and multiplications, the TCMLIB, replaced later by LaPack. Of course, it was the debugging phase that, as usual, took the huge amount of time. The problem was that I had no code to compare with, or to look inside. The only comparison was the very last product: the LDA and RPA published dielectric constants of Silicon calculated with this or the other method. Finally, after a lot of work, all the bugs were identified and removed, and the dielectric constants produced were exactly those published.

The results with the non-local kernels however, were not such encouraging. Little improvements. After many trials using non-local kernels based on all the possible jellium local-field corrections (Hubbard, Geldart-Taylor, Vashishta-Singwi), all giving small changes with respect to RPA, it became clear that a local or a slightly non-local kernel can't have an effect. Indeed, in TDDFT the kernel appears always coupled to a χº which goes to 0 as O(q²) for q→0. It appeared necessary to have a ultra-nonlocal kernel showing a O(1/q²) behaviour for q→0, such as multiplied by χº the result is finite. The jellium kernel misses such a long-range term, and the same holds for all the jellium based approximations. This finding was also suggested in a paper by Gonze, Ghosez and Godby, who indicated a relation between this 1/q² divergence and the GW correction to the band-gap. Here however, on the optical spectra, excitonic effects are also playing a role.

I started then to implement a long-range kernel, with a 1/q² divergence weighted by an adjustable parameter, that I called α. And, as expected, this kernel did have an effect on the final spectra. Indeed it seemed to reproduce the effect seen when opening a gap, like in GW. But there was also a collateral effect on the spectrum shape, and I found no α able to improve on the shape the agreement with the experiment.

Although the results on optical spectra were disappointing, I was surprised by the good results on the Energy-Loss function (EELS), at the LDA level but even at the RPA. I started hence to calculate Energy-Loss spectra also at finite q, to compare with IXSS spectra. And to a much more exotic technique, CIXS. And again, the results were surprisingly good. The explanation of such result came only when I moved to Palaiseau, for a Post-Doc with Lucia Reining. Thank to the Bethe-Salpeter theory and the code she had developped, we could find that GW and excitonic effects cancel each other in the energy-loss function, contrary to optical absorption.

Lucia Reining

But in Palaiseau Lucia Reining had also the crucial idea to solve the kernel issue. She suggested to separate the two kernel tasks (and problems), that is the self-energy and the excitonic parts. She suggested to start from an already GW corrected polarizability, and to let to the long-range kernel only the task to reproduce the excitonic effects. Of course, in this case we had to assume negative values for the α weight. And this, with my surprise and skepticism, worked immediately!

My stay in Palaiseau allowed also to compare directly the results of the DP code and the Bethe-Salpeter EXC. And to verify eventual residual bugs in both. It became a question of principle between Lucia Reining and me what was the more correct code (or better to say, the less bugged). At the end we discovered a tiny difference in the results, when sampling at very fine frequency grids. But with my disappointement, it was the EXC code that was right. And Lucia Reining even indicated the bug in DP, and corrected it, such the two codes now gave exactly overimposable curves. Well, it was a reassuring confirmation that our previous published results, both mine on TDDFT and Lucia Reining's on Bethe-Salpeter, were correct, ufh.

The code evolved, toward a more and more user-friendly form, as I always targetted. Already in Rome, to solve the porting problems on the different machines in use there, an interface-to-libraries layer, driven by the makefile, was introduced, and the compilation became machine independent. I introduced the C shell to make it a really professional product, with a Unix command-like interface to the user. The C written parser also improved on this, with the introduction of a set of driving commands, later uniformized with those of ABINIT. This indeed allowed a lot of people to take the code, and use it with little effort for me to explain how it works; and for them to learn it.

The next important contribution to the code, however, came thank to Francesco Sottile and his thesis' work. He implemented in the code what later has been named the Nanoquanta kernel, which has finally solved the problem to have an ab initio TDDFT kernel able to reproduce excitonic effects on optical spectra, at a level of accuracy comparable to Bethe-Salpeter results. This kernel was a derivation done by Lucia Reining starting from a comparison between the Bethe-Salpeter equation and the TDDFT Casida's equation. But I must say that the same kernel was also previously derived by Rodolfo Del Sole along perturbation theory approximations holding into many-body theory. He showed me, before to leave to Palaiseau, this kernel and the promising results he got together with Gianni Adragna using a Tight Binding approach. But at that time the statistics was only on Silicon and GaAs and I was erroneously believing that the Tight-Binding adjustable parameters could do most of the job.

When we realized that the two kernels were in fact the same, we started to believe that an ab initio implementation was really worth to. This kernel has a really complicated, orbital dependent, shape. It costed Francesco also a lot of time in debugging. But at the end he succeeded! And the surprising nice result was that TDDFT and the DP code were even able to reproduce in solid Argon a series of three bound excitons within the photoemission band-gap, in agreement with the experiment and the Bethe-Salpeter result. When Francesco showed us this surprising result, we exclamed: Well so far TDDFT was considered to be the future. But now TDDFT is the present!

Francesco Sottile

Francesco Sottile also contributed with a further modernization of the code, introduction of an autoconfig to setup the compilation parameters, use of Perl scripts for the battery of tests, and so on. And he is also contributing in the optimization of this Nanoquanta kernel part.

Finally the code received also contributions from Fabien Bruneval, an attentive bug chaser in this as in the other codes I was mantaining. He also improved the efficiency of the code thank to the introduction in the pseudopotential part of the completely non-local representation instead of the semi-local. And from Andrea Cucca, who contributed to the parallelization of the code, still ongoing.