Dear Federico,

As the name suggests, the fundamental quantity in DFT is the density, which all observable quantities are functionals of.

The Kohn-Sham approach is simply a clever way to generate the density n(ṟ), based on the assumption that there is a system of non-interacting particles which has everywhere the same value of n as the interacting one of interest. For non-interacting Fermions in the ground state, the N-particle wave-function is indeed a single Slater determinant composed of the lowest N single-particle orbitals. These orbitals however have "a priori" no physical meaning. They are only part of the tool needed to construct the density as sum over their modulus squared and, since the particles don't interact, there is no exchange energy for that system.

I enclose a few pages from a book I wrote with Veljko Zlatic on the Modern Theory of Thermoelectricity, which, I hope will be helpful to you.

Best regards, René Monnier

Thank you very much. I have read it and consulted some more and I think my questions are partly answered. I'm leaving this for someone else who might wonder:

-The solution wavefunction for a K-S system is a Slater determinant, yes. But that doesn't matter very much since what we want is the density, as the sum of the density of all individual wavefunctions.

-The K-S (or should I say Roothaan) equations are solved simultaneously as a system and that yields the LCAO coefficients and the orbital energy eigenvalues simultaneously.

Still unsure why the molecular wavefunction does take into account electron exchange by being a Slater orbital but at the same time it is derived using the Vxc potential which already includes it.

At the same time, the simple sum of each MO's densities clearly isn't a slater determinant so the exchange effect comes only from vxc.

I think, there is still some confusion here!

Your first point is correct. You want to construct the density of the true (interacting) system from the orbital densities of a non-interacting system

in an external potential to be determined.

Your second point is incorrect. The Rothaan equations solve the Hartree-Fock problem in a given basis. This is NOT what is done in Kohn-Sham.

The Slater determinantal form simply takes into account the fact that you deal with FERMIONS, by antisymmetrizing the product of single-particle orbitals, and NOTHING MORE.

The above remark makes your third point irrelevant.

As to your last point: it is easy to show that, if your N-particle wave function

consists of a single Slater determinant of single-particle orbitals, the density of the system is simply given by the sum over the N single-particle densities.

Thanks again for your answer. I think I understand now the Slater question.

Back to the "second point" though, do you mean that the Roothaan method is not used in K-S DFT? Mind that i'm talking about calculations on molecules and atomic basis sets.

Although most of literature explain the Roothaan method as applied to H-F, almost none of the books I have give much details as to how the K-S equations are solved by the software. I have Postnikov's lecture 4 ( www.home.uni-osnabrueck.de/apostnik/lectures.html ) that applies Roothaan to K-S, I have Nogueira, Castro and Marques' "a tutorial on density functional theory" (Perdew "density functional theory section 6.2.1), that says that the K-S equations are solved by "diagonalizing the hamiltonian matrix", and some more, so I understand it is used in K-S as well.

Specifically, I would like to know if Gaussian uses this method or what method does it use.

From a technical point of view, the KS equations and the HF equations are solved in exactly the same way. Both are using the approach of Roothaan and Hall (RH), i.e. using a basis set expansion. The reason why the RH equations were originally formulated only for HF is that KS-DFT was not formulated at that time. But all DFT codes I know first build the HF or KS matrix, i.e. the first derivative of the energy w.r.t. the orbitals or the electron density, respectively, using the procedure of RH (but not necessarily the original algorithm for solving the equations, but most likely some improvements like DIIS or quadratic convergent SCF, etc.). Regarding the theory, the original KS and HF equations just differ in the way how the exchange-correlation (XC) energy is calculated. In HF the exact-exchange energy, given by the Slater determinant, is used, while a density functional is employed for that purpose in KS-DFT, e.g. BLYP. However, strictly speaking, many of the DFT calculations nowadays are no 'real' KS-DFT calculations anymore. In particular, often hybrid or meta-GGA functionals are used, which include the kinetic energy density and exact exchange, which both are only implicit functionals of the density, but explicit functionals of the orbitals and the density matrix. Accordingly, derivatives w.r.t. the electron density as required for 'real' KS-DFT are not simply achievable. Therefore, one would need to perform an optimized effective potential (OEP) calculation, which is not the standard procedure in most DFT codes. Instead, the KS matrix is usually build like in HF by taking the first derivative w.r.t. the density matrix. Such a treatment is legitimized within the generalized KS framework (in contrast to the original KS), in which the non-interacting reference system is replaced by a partially interacting reference system which can be expressed by a single Slater determinant (SD), thus allowing the direct HF-like treatment of orbital-dependent XC functionals. In that sense HF is simply an XC functional within generalized KS, which is also known as HF-KS. In summary, generalized KS-DFT and HF are more or the less equivalent. That's why techniques to solve the corresponding equations are exactly the same in most QC programs (such as Gaussian, Turbomole and many others, but the list would be too long). Only the way the integrals are calculated usually differ. DFT requires a grid integration for the determination of the XC integrals, while exact exchange integrals in HF can be calculated using solely 'analytical' methods (or better to say, numerical methods almost equivalent to analytical results). However, also exact exchange can be calculated on a numerical grid, which nicely highlights similarities between the methods.

I hope, I didn't enhance confusion by adding too many details.

Dear Federico,

I am glad that the conceptual problems have been solved and that we are now down at the technical level.

What I meant in my "second point" was, that Rothaan originally wanted to solve the Hartre-Fock equations for molecules, a difficult problem due to to the non-local nature of the exchange potential. In KS one deals with a the much simpler case of a local potential.

The diagonalization of the Hamiltonian in a given basis set is a fairly standard procedure, and I would not necessarily associate it with the

name of Rothaan, although it can be done.

I have no experience with Gaussian, but I suppose that the name specifies the form of the basis set, and then the answer would be "yes".

Thanks mr. Mayer for your very detailed response. The first half of it is exactly what I figured. The second part, about the derivatives respect to the density and on, I had no idea and none of my bibliography said that. Mainly, they say, a kind of a "chain rule" is used to obtain the derivatives analitically. From what I gathered from the Siesta manual and the Gaussian manual, in DFT, the first derivatives are calculated analitically. The second derivatives are not usually calculated analitically in a geometry optimization, but they are in a vibrational calculation.

Thanks also mr. Monnier for your responses.

Dear Frederico,

one should be carful not to mix the different types of derivatives. For the single-point DFT calculation, i.e. the self-consistent field calculation, one needs to calculate the first derivative with respect to the density matrix (or orbitals or density). And as the manuals state, one simply uses a "kind of chain rule". In more detail, one constructs kind of a complete differential for the mentioned derivative (density matrix or orbital) employing the functional chain rule. In this way you can express your derivatives with respect to the density matrix (or orbitals) in terms of the ingredients of the exchange-correlation functional. Unfortunately, respective formulas are often not shown in the literature.

The second kind of derivatives you mentioned (which is usually called only 'derivatives') are derivatives of the total energy with respect to the molecular coordinates. These are not the derivatives needed to converge the orbital coefficients within the self-consistent-field calculation (the single-point calculation). As you stated, derivatives with respect to the molecular coordinates are required e.g. for structure optimizations. Here, usually quasi-Newton-Raphson methods are used, which need a start guess for the model Hessian at the initial structure and then update the model Hessian for the subsequent subsequent structures (during optimization) using only gradients. This is more efficient than calculating the real Hessian for each structure. However, the model Hessian is only suited for structure optimizations, since individual elements are not accurate enough.

For minimum structure optimizations, one usually uses the BFGS algorithm for the update of the Hessian and a diagonal Hessian (either with values of 1 or a fixed guess based on some parameters) is usually sufficient as initial guess. For transition-state optimizations, a more accurate Hessian is needed and one has to use the real Hessian (second derivatives) as initial guess. Also other update procedures are used. For the calculation of 'classical' vibrational frequencies, you also need the real Hessian (not a model Hessian).

However, such a Hessian does not need to be calculated analytically. Although analytical second derivatives are usually used for DFT, since respective modules are usually available in most quantum-chemical programs and it is still cheaper than the alternative (although it is far more expensive than analytical first derivatives), one can also perform a numerical determination of the Hessian using using finite displacements of the molecular coordinates by just using energies (and analytical gradients, which would reduce the required number of displacements).

Sorry for the long remark, but it seemed to me that there were some unclear points regarding the term 'derivatives'.

Best regards.

Thanks Mr. Maier for your detailed response. I think this belongs in another question of mine, though.

Do you mean that analytical second derivatives are computationally more expensive than numerical second derivatives? (respect to the nuclei coordinates). Thanks!

Dear Federico,

The answer to your latest question strongly depends a) on the size of your system, b) the employed method and c) the available computational resources.

In general, one has three options of calculating second derivatives:

1. fully numerical (energies only)

2. semi-numerical (analytical gradients)

3. analytical

For 1., the number of displacements scales as N^2, which, unless one deals with small molecules, is more expensive than the other two options. For expensive electronic methods without analytical gradients this, however, is the only option to do frequency calculations.

For 3., one only needs the molecular structure (so, no displacements). However, evaluation of analytical second derivatives is usually more expensive than the single-point calculation and often not available for high-level methods. Furthermore, memory requirements are usually increased, which might be problematic for large molecules (and already memory consuming methods). In most cases, especially for HF and KS-DFT this is a reasonable option.

For 2. , the number of displacements scales only as N. Hence, also relatively large molecules can be calculated without an exploding calculation cost. Furthermore, analytical gradients are often relatively cheap and do not need much more memory than the single-point calculation, especially for HF/KS-DFT (but perhaps not for more sophisticated methods such as CASPT2). Hence, this is a good compromise between 1. and 3.. Nonetheless, this option is usually worse than 3. regarding used CPU time (except perhaps for small molecules). In contrast to 3., option 2., however, is easily parallelizable, so that in terms of used life time this sometimes might be the better option (although some CPU time might be burned), especially if one has access to a lot of CPUs with low memory, e.g. on a computing cluster.

At the end, the actual performance strongly depends on your problem (molecule, method) and the used resources (program, computers). Hence, in my opinion there is no unique answer to your question (although there might be some more guidelines than the ones I mentioned here). I would say, this is also related to some kind of experience.

Best regards