Exchange-correlation potential in pyscf

[AlexanderFreeman]

Hello! Could you please explain me how does pyscf calculate exchange-correlation potential in GGA case? I analyzed code and it seems to me that pyscf uses only first-order gradients of xc energy density w.r.t. rho from libxc (vrho, vgamma). I found in file "/dft/numint.py" a function "nr_rks" and in case of GGA, for example, it slices the result obtained from "eval_xc" function:
exc, vxc = ni.eval_xc(xc_code, rho, spin=0, relativity=relativity, deriv=1, verbose=verbose)[:2]
And it uses only exc, vxc, but discards kxc and fxc. Does it mean that pyscf somehow calculates second-order gradients needed for exchange-correlation potential calculation? Or it uses another way for calculating exchange-correlation potential?

[gkc1000]

The standard definition of GGA only requires first derivatives of rho. Meta-gga's requires the kinetic energy density, and TDDFT requires fxc.

…

On Thu, Jul 2, 2020 at 7:50 AM AlexanderFreeman ***@***.***> wrote: Hello! Could you please explain me how does pyscf calculate exchange-correlation potential? I analyzed code and it seems to me that pyscf uses only first-order gradients from libxc (vrho, vgamma). I found in file "/dft/numint.py" a function "nr_rks" and in case of GGA, for example, it slices the result obtained from "eval_xc" function: exc, vxc = ni.eval_xc(xc_code, rho, spin=0, relativity=relativity, deriv=1, verbose=verbose)[:2] And it uses only exc, vxc, but discards kxc and fxc. Does it mean that pyscf somehow calculates second-order gradients needed for exchange-correlation potential calculation? Or it uses another way for calculating exchange-correlation potential? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#650>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRDSPSWTZI5WVPJ2PS3RZSNEHANCNFSM4OO62ZMA> .

[AlexanderFreeman]

The standard definition of GGA only requires first derivatives of rho. Meta-gga's requires the kinetic energy density, and TDDFT requires fxc.
…
On Thu, Jul 2, 2020 at 7:50 AM AlexanderFreeman @.***> wrote: Hello! Could you please explain me how does pyscf calculate exchange-correlation potential? I analyzed code and it seems to me that pyscf uses only first-order gradients from libxc (vrho, vgamma). I found in file "/dft/numint.py" a function "nr_rks" and in case of GGA, for example, it slices the result obtained from "eval_xc" function: exc, vxc = ni.eval_xc(xc_code, rho, spin=0, relativity=relativity, deriv=1, verbose=verbose)[:2] And it uses only exc, vxc, but discards kxc and fxc. Does it mean that pyscf somehow calculates second-order gradients needed for exchange-correlation potential calculation? Or it uses another way for calculating exchange-correlation potential? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#650>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRDSPSWTZI5WVPJ2PS3RZSNEHANCNFSM4OO62ZMA .

I mean not gradients of rho, but gradients of xc energy density w.r.t. rho. In gga case vxc (exchange-correlation potential) depends not only on dexc/dp, dexc/dsigma, but also on d2exc/dpdsigma, d2exc/dsigma^2. These gradients are included in fxc.

[AlexanderFreeman]

I used the GGA Vxc definition from article : https://www.researchgate.net/publication/1826446_Evaluation_of_Exchange-Correlation_Energy_Potential_and_Stress . You can find it on page 2 (Vxc gga case).

[gkc1000]

Do you mean, in the functional derivative, the term that looks like vxc~ = \nabla. dE/d(|\nabla \rho|)? If I recall correctly, for the vxc matrix this is done by integration by parts, e.g. (a | vxc~ | b) --> (\nabla a | dE/d(\nabla \rho) | b). You can see something to that effect on line 508. On Thu, Jul 2, 2020 at 9:10 AM AlexanderFreeman <notifications@github.com> wrote:

…

The standard definition of GGA only requires first derivatives of rho. Meta-gga's requires the kinetic energy density, and TDDFT requires fxc. … <#m_-1143627517613182789_> On Thu, Jul 2, 2020 at 7:50 AM AlexanderFreeman *@*.***> wrote: Hello! Could you please explain me how does pyscf calculate exchange-correlation potential? I analyzed code and it seems to me that pyscf uses only first-order gradients from libxc (vrho, vgamma). I found in file "/dft/numint.py" a function "nr_rks" and in case of GGA, for example, it slices the result obtained from "eval_xc" function: exc, vxc = ni.eval_xc(xc_code, rho, spin=0, relativity=relativity, deriv=1, verbose=verbose)[:2] And it uses only exc, vxc, but discards kxc and fxc. Does it mean that pyscf somehow calculates second-order gradients needed for exchange-correlation potential calculation? Or it uses another way for calculating exchange-correlation potential? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#650 <#650>>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRDSPSWTZI5WVPJ2PS3RZSNEHANCNFSM4OO62ZMA . I mean not gradients of rho, but gradients of xc energy density w.r.t. rho. In gga case vxc (exchange-correlation potential) depends not only on dexc/dp, dexc/dsigma, but also on d2exc/dpdsigma, d2exc/dsigma^2. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#650 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRC4UABQPPY6EF3W3U3RZSWQZANCNFSM4OO62ZMA> .

[gkc1000]

Ok, so if you are looking for the value of vxc on a grid specifically, this quantity is not used to evaluate (a|vxc|b) for GGA's. Although you could reconstruct it from the libxc interface.

…

On Thu, Jul 2, 2020 at 9:19 AM AlexanderFreeman ***@***.***> wrote: I used the GGA definition from article : https://www.researchgate.net/publication/1826446_Evaluation_of_Exchange-Correlation_Energy_Potential_and_Stress <http://url> — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#650 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRD4ZDCL5CMDA3YXQRDRZSXQ5ANCNFSM4OO62ZMA> .

[AlexanderFreeman]

Ok, so if you are looking for the value of vxc on a grid specifically, this quantity is not used to evaluate (a|vxc|b) for GGA's. Although you could reconstruct it from the libxc interface.
…
On Thu, Jul 2, 2020 at 9:19 AM AlexanderFreeman @.***> wrote: I used the GGA definition from article : https://www.researchgate.net/publication/1826446_Evaluation_of_Exchange-Correlation_Energy_Potential_and_Stress http://url — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#650 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRD4ZDCL5CMDA3YXQRDRZSXQ5ANCNFSM4OO62ZMA .

Thank you for your answers! Everything becomes much clearer.

[AlexanderFreeman]

Ok, so if you are looking for the value of vxc on a grid specifically, this quantity is not used to evaluate (a|vxc|b) for GGA's. Although you could reconstruct it from the libxc interface.
…
On Thu, Jul 2, 2020 at 9:19 AM AlexanderFreeman @.***> wrote: I used the GGA definition from article : https://www.researchgate.net/publication/1826446_Evaluation_of_Exchange-Correlation_Energy_Potential_and_Stress http://url — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#650 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRD4ZDCL5CMDA3YXQRDRZSXQ5ANCNFSM4OO62ZMA .

Can pyscf be used for vxc potential calculation on specific grid? Is procedure for that implemented in pyscf?

[gkc1000]

To my knowledge, not for GGA. The integration by parts using AO derivatives, and use of vxc(r) itself on the grid are not the same numerically (due to grid integration error) and as far as I recall, using AO derivatives was found to be better in the early 90's in the early implementations in Gaussian and CADPAC. So, vxc(r) on the grid is mainly useful for real-space analysis of DFT functionals. On Thu, Jul 2, 2020 at 10:41 AM AlexanderFreeman <notifications@github.com> wrote:

…

Ok, so if you are looking for the value of vxc on a grid specifically, this quantity is not used to evaluate (a|vxc|b) for GGA's. Although you could reconstruct it from the libxc interface. … <#m_-2845581592048752796_> On Thu, Jul 2, 2020 at 9:19 AM AlexanderFreeman *@*.***> wrote: I used the GGA definition from article : https://www.researchgate.net/publication/1826446_Evaluation_of_Exchange-Correlation_Energy_Potential_and_Stress http://url — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#650 (comment) <#650 (comment)>>, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAN5NRD4ZDCL5CMDA3YXQRDRZSXQ5ANCNFSM4OO62ZMA . Can pyscf be used for vxc potential calculation on specific grid? Is procedure for that implemented in pyscf? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#650 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRERSRQOQTSQWPQVO43RZTBC3ANCNFSM4OO62ZMA> .

[sunqm]

The xc potential and the deriviatives of xc potential on grids are provided by libxc library or xcfun library. You can first compute the electron density then pass to the eval_xc https://github.com/pyscf/pyscf/blob/master/pyscf/dft/libxc.py to evaluate the xc functional on real space grids.

…

On Thu, Jul 02, 2020 at 09:33:04AM -0700, AlexanderFreeman wrote: > Ok, so if you are looking for the value of vxc on a grid specifically, this quantity is not used to evaluate (a|vxc|b) for GGA's. Although you could reconstruct it from the libxc interface. > […](#) > On Thu, Jul 2, 2020 at 9:19 AM AlexanderFreeman ***@***.***> wrote: I used the GGA definition from article : https://www.researchgate.net/publication/1826446_Evaluation_of_Exchange-Correlation_Energy_Potential_and_Stress <http://url> — You are receiving this because you commented. Reply to this email directly, view it on GitHub <[#650 (comment)](#650 (comment))>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAN5NRD4ZDCL5CMDA3YXQRDRZSXQ5ANCNFSM4OO62ZMA> . Thank you for your answers! Everything becomes much clearer. Does pyscf calculate only vxc matrix or can be used for vxc potential calculation on specific grid? -- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: #650 (comment)

--

[ajz34]

I've probably reproduced GGA vxc(r) on specific grid, using fxc from PySCF's LibXC interface (NumInt).
https://ajz34.readthedocs.io/zh_CN/latest/QC_Notes/XC_Potential_Series/PBE_xc_potential.html
Code's in the blog. It's written in Chinese, but mostly equations and code. Hope not be confused by language (0.0)
It's default grid on Ne atom, 6-311G. Restricted PBE.
Figure 5 of V. N. Staroverov JCTC 2013 was reproduced, though this code is not reproducing concept of that article directly.
I think equation (2) of L. C. Balbas PRB 2001 is very similar to expression of vxc(r) in my blog.
Expression of vxc(r) in my blog was derived independently, so no reference for that.

[ajz34]

To my understanding, doing PBE SCF process does need vxc matrix in atomic orbital form, and does not require fxc; but generating vxc(r) does need fxc in GGA's case.
For LDA, only vxc (i.e. xc kernel's first derivative) is required for vxc(r). Also, vxc(r) is actually functional of rho.
It's really tricky for me 😹

[zhendongli2008]

Hi Zhenyu, As answered by others before, vxc(r) in GGA is usually not computed for SCF, since an equivalent simple expression for the contribution to the Fock matrix can be obtained using integration by part. If you really need it, the chain rule can be used to derive an expression for vxc(r), which for example can be found in Eqs. (A1-A6) [ https://aip.scitation.org/doi/10.1063/1.3479401], where we used it for the XC kernel in spin-flip TD-DFT. Best, Zhendong Li Zhenyu Zhu ajz34 <notifications@github.com> 于2020年8月12日周三 下午8:20写道：

…

To my understanding, doing PBE SCF process does need vxc matrix in atomic orbital form, and does not require fxc; but generating vxc(r) does need fxc in GGA's case. For LDA, only vxc (i.e. xc kernel's first derivative) is required for vxc(r). Also, vxc(r) is actually functional of rho. It's really tricky for me 😹 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#650 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVP65PRFHPL4NQLYKORVDTSAKCJHANCNFSM4OO62ZMA> .

[zhendongli2008]

Sorry, the correct reference should be https://doi.org/10.1063/1.3676736. Zhendong Li <zhendongli2008@gmail.com> 于2020年8月13日周四 上午11:58写道：

…

Hi Zhenyu, As answered by others before, vxc(r) in GGA is usually not computed for SCF, since an equivalent simple expression for the contribution to the Fock matrix can be obtained using integration by part. If you really need it, the chain rule can be used to derive an expression for vxc(r), which for example can be found in Eqs. (A1-A6) [ https://aip.scitation.org/doi/10.1063/1.3479401], where we used it for the XC kernel in spin-flip TD-DFT. Best, Zhendong Li Zhenyu Zhu ajz34 ***@***.***> 于2020年8月12日周三 下午8:20写道： > To my understanding, doing PBE SCF process does need vxc matrix in atomic > orbital form, and does not require fxc; but generating vxc(r) does need > fxc in GGA's case. > For LDA, only vxc (i.e. xc kernel's first derivative) is required for > vxc(r). Also, vxc(r) is actually functional of rho. > It's really tricky for me 😹 > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub > <#650 (comment)>, or > unsubscribe > <https://github.com/notifications/unsubscribe-auth/ABVP65PRFHPL4NQLYKORVDTSAKCJHANCNFSM4OO62ZMA> > . >