Efficient HF exchange evaluation through Fourier convolution in Cartesian grid for orbital-dependent density functionals

TL;DR

This paper introduces a fast, purely numerical Cartesian grid method for calculating Hartree-Fock exchange in density functional theory, utilizing Fourier convolution and range-separated Coulomb kernels for improved efficiency and scalability.

Contribution

It presents a novel Fourier convolution-based numerical approach for HF exchange evaluation that scales logarithmically with molecular size, extending to global hybrid functionals within pseudopotential Kohn-Sham theory.

Findings

Excellent agreement with semi-numerical methods in total and orbital energies

Logarithmic scaling with molecular size demonstrated

Potential for developing advanced range-separated hybrid functionals

Abstract

We present a purely numerical approach in Cartesian grid, for efficient computation of Hartree-Fock (HF) exchange contribution in the HF and density functional theory models. This takes inspiration from a recently developed algorithm [Liu \emph{et al.}, J.~Chem.~Theor.~Comput.~ \textbf{13}, 2571 (2017)]. A key component is the accurate evaluation of electrostatic potential integral, which is the rate-determining step. This introduces the Fourier convolution theorem in conjunction with a range-separated Coulomb interaction kernel. The latter is efficiently mapped into real grid through a simple optimization procedure, giving rise to a constraint in the range-separated parameter. The overall process offers logarithmic scaling with respect to molecular size. It is then extended towards global hybrid functionals such as B3LYP, PBE0 and BHLYP within pseudopotential Kohn-Sham theory, through an LCAO-MO ansatz in Cartesian grid, developed earlier in our laboratory. For sake of comparison, a parallel semi-numerical approach has also been worked out that exploits the familiar Obara-Saika recursion algorithm. An excellent agreement between these two routes is demonstrated through total energy and orbital energy in a series of atoms and molecules (including 10 $\pi$-electron molecules), employing an LANL2DZ-type basis function. A critical analysis of these two algorithms reveals that the proposed numerical scheme could lead to very attractive and competitive scaling. The success of our approach also enables us for further development of optimally tuned range-separated hybrid and hyper functionals.

Tables

Table 1. Table 1: Convergence of HF energy †,‡ of HCl ( R = 1.275 ​ Å 𝑅 1.275 Å R=1.275\textup{\AA} ) in the grid ( h r = 0.3 subscript ℎ 𝑟 0.3 h_{r}=0.3 ). All results are in a.u.

Set I					Set II
$N_x$	$N_y$	$N_z$	$⟨E⟩$		$N_x$	$N_y$	$N_z$	$⟨E⟩$
32	32	32	$-$15.24881		36	36	60	$-$15.27363
–	–	36	$-$15.26277		40	40	–	$-$15.27496
–	–	40	$-$15.26604		44	44	–	$-$15.27521
–	–	44	$-$15.26684		46	46	–	$-$15.27525
–	–	48	$-$15.26703		48	48	–	$-$15.27526
–	–	52	$-$15.26707		50	50	–	$-$15.27527
–	–	56	$-$15.26708		52	52	–	$-$15.27527
–	–	60	$-$15.26708		54	54	–	$-$15.27527
–	–	64	$-$15.26708		56	56	–	$-$15.27527

Table 2. Table 2: Convergence of KS energy †,‡ of Cl 2 ( R = 4.2 𝑅 4.2 R=4.2 a.u.) in the grid ( h r = 0.3 subscript ℎ 𝑟 0.3 h_{r}=0.3 ) using B3LYP XC functional. All results are in a.u.

Set I					Set II
$N_x$	$N_y$	$N_z$	$⟨E⟩$		$N_x$	$N_y$	$N_z$	$⟨E⟩$
32	32	32	$-$28.49544		36	36	76	$-$29.78914
–	–	40	$-$29.59219		40	40	–	$-$29.79282
–	–	48	$-$29.75311		44	44	–	$-$29.79363
–	–	52	$-$29.76710		48	48	–	$-$29.79380
–	–	56	$-$29.77140		52	52	–	$-$29.79384
–	–	68	$-$29.77304		56	56	–	$-$29.79384
–	–	72	$-$29.77306		60	60	–	$-$29.79384
–	–	76	$-$29.77306		64	64	–	$-$29.79384
–	–	80	$-$29.77306		68	68	–	$-$29.79384

Table 3. Table 3: HF, B3LYP PBE0 energies (a.u.) of selected atoms. E diff = | E ​ NR − E SNR | subscript 𝐸 diff 𝐸 NR subscript 𝐸 SNR E_{\mathrm{diff}}=|E{\mathrm{NR}}-E_{\mathrm{SNR}}| .

Atom	$-⟨E⟩$
	HF			B3LYP			PBE0
	NR	SNR	Ediff	NR	SNR	Ediff	NR	SNR	Ediff
Be	0.96019	0.96019	0.00000	0.99386	0.99386	0.00000	0.99598	0.99598	0.00000
B	2.53426	2.53426	0.00000	2.59552	2.59552	0.00000	2.59751	2.59751	0.00000
C	5.30354	5.30354	0.00000	5.39301	5.39301	0.00000	5.39603	5.39603	0.00000
N	9.61972	9.61969	0.00003	9.73283	9.73282	0.00001	9.74190	9.74189	0.00001
O	15.61720	15.61681	0.00039	15.80556	15.80549	0.00007	15.80351	15.80341	0.00010
Al	1.86947	1.86947	0.00000	1.92379	1.92379	0.00000	1.93075	1.93075	0.00000
Si	3.67570	3.67570	0.00000	3.75077	3.75077	0.00000	3.76281	3.76281	0.00000
P	6.31523	6.31523	0.00000	6.40580	6.40580	0.00000	6.42550	6.42550	0.00000
S	9.87464	9.87464	0.00000	10.01827	10.01827	0.00000	10.03519	10.03519	0.00000
Cl	14.68130	14.68129	0.00001	14.87170	14.87170	0.00000	14.88873	14.88873	0.00000
Ga	1.94526	1.94526	0.00000	2.00195	2.00195	0.00000	2.00831	2.00831	0.00000
Ge	3.59814	3.59814	0.00000	3.67466	3.67466	0.00000	3.68693	3.68693	0.00000
As	5.95615	5.95615	0.00000	6.04808	6.04808	0.00000	6.06823	6.06823	0.00000
Se	9.01108	9.01108	0.00000	9.15439	9.15439	0.00000	9.17235	9.17235	0.00000
Br	12.91872	12.91872	0.00000	13.10652	13.10652	0.00000	13.12535	13.12535	0.00000

Table 4. Table 4: HF, B3LYP PBE0 energies (a.u.) of molecules. E diff = | E ​ NR − E SNR | subscript 𝐸 diff 𝐸 NR subscript 𝐸 SNR E_{\mathrm{diff}}=|E{\mathrm{NR}}-E_{\mathrm{SNR}}| .

Molecule	$-⟨E⟩$
	HF			B3LYP			PBE0
	NR	SR	Ediff	NR	SR	Ediff	NR	SR	Ediff
Cl2	29.34810	29.34809	0.00001	29.79384	29.79384	0.00001	29.82896	29.82896	0.00000
Br2	25.82400	25.82400	0.00000	26.25227	26.25227	0.00000	26.29481	26.29481	0.00000
I2	22.30720	22.30720	0.00000	22.71498	22.71498	0.00000	22.76592	22.76592	0.00000
HCl	15.27527	15.27527	0.00000	15.50844	15.50843	0.00001	15.52803	15.52803	0.00000
HBr	13.49506	13.49506	0.00000	13.72537	13.72537	0.00000	13.74603	13.74603	0.00000
HI	11.72279	11.72279	0.00000	11.94619	11.94619	0.00000	11.96939	11.96939	0.00000
H2S	11.02605	11.02605	0.00000	11.25511	11.25511	0.00000	11.27394	11.27394	0.00000
H2Se	10.14548	10.14548	0.00000	10.37371	10.37371	0.00000	10.39215	10.39215	0.00000
PH3	8.01636	8.01636	0.00000	8.23350	8.23350	0.00000	8.24942	8.24942	0.00000
CH3Cl	21.87306	21.87306	0.00000	22.29768	22.29768	0.00000	22.33525	22.33525	0.00000
CH4	7.78878	7.78888	0.00010	8.00843	8.00846	0.00003	8.02684	8.02686	0.00002
SiH3Cl	20.19863	20.19862	0.00001	20.58442	20.58441	0.00001	20.61885	20.61885	0.00000
Si2H6	10.93377	10.93377	0.00000	11.28249	11.28249	0.00000	11.31207	11.31207	0.00000
P4	25.13843	25.13843	0.00000	25.81974	25.81974	0.00000	25.91452	25.91452	0.00000
SiH4	6.03289	6.03289	0.00000	6.22621	6.22621	0.00000	6.23634	6.23634	0.00000

Table 5. Table 5: Negative HOMO energies, − ϵ HOMO subscript italic-ϵ HOMO -\epsilon_{\mathrm{HOMO}} (in a.u.) for selected atoms and molecules using HF, B3LYP, PBE0 and BHLYP XC functionals.

Atom	$-{ϵ}_{\mathrm{HOMO}}$(a.u.)					Molecule	$-{ϵ}_{\mathrm{HOMO}}$(a.u.)
	HF	B3LYP	PBE0	BHLYP	Expt.†		HF	B3LYP	PBE0	BHLYP	Expt.‡
Be	0.3090	0.2291	0.2387	0.2650	0.3426	Cl2	0.4786	0.3274	0.3433	0.3947	0.4219
B	0.3080	0.1830	0.1929	0.2348	0.3050	Br2	0.4139	0.2889	0.3034	0.3456	0.3864
C	0.4313	0.2615	0.2771	0.3290	0.4138	I2	0.3695	0.2658	0.2793	0.3138	0.3429
N	0.5657	0.3497	0.3707	0.4338	0.5341	HCl	0.4772	0.3274	0.3425	0.3937	0.4700
O	0.5098	0.3247	0.3379	0.4090	0.5005	HBr	0.4302	0.3030	0.3170	0.3602	0.4303
Al	0.2095	0.1291	0.1408	0.1621	0.2200	HI	0.3875	0.2804	0.2932	0.3295	0.3815
Si	0.3029	0.1937	0.2088	0.2374	0.3000	H2S	0.3895	0.2617	0.2746	0.3191	0.3851
P	0.3901	0.2545	0.2732	0.3076	0.3854	H2Se	0.3613	0.2466	0.2593	0.2987	0.3649
S	0.3631	0.2506	0.2625	0.3043	0.3807	PH3	0.3849	0.2675	0.2781	0.3200	0.3626
Cl	0.4731	0.3294	0.3442	0.3944	0.4766	CH3Cl	0.4340	0.2946	0.3084	0.3567	0.4149
Ga	0.2058	0.1263	0.1381	0.1590	0.2205	CH4	0.5416	0.3882	0.4013	0.4555	0.4998
Ge	0.2844	0.1825	0.1974	0.2232	0.2903	SiH3Cl	0.4509	0.3149	0.3288	0.3754	0.4281
As	0.3665	0.2426	0.2605	0.2910	0.3607	Si2H6	0.4068	0.3043	0.3152	0.3516	0.3870
Se	0.3319	0.2337	0.2451	0.2815	0.3584	P4	0.3844	0.2921	0.3075	0.3362	0.3381
Br	0.4206	0.3010	0.3144	0.3563	0.4341	SiH4	0.4871	0.3576	0.3687	0.4149	0.4520

Table 6. Table 6: Negative HOMO energies, − ϵ HOMO subscript italic-ϵ HOMO -\epsilon_{\mathrm{HOMO}} (in a.u.) for selected π 𝜋 \pi -electron molecules using HF, B3LYP, PBE0 and BHLYP XC functionals.

Molecule	$-{ϵ}_{\mathrm{HOMO}}$(a.u.)
	HF	B3LYP	PBE0	BHLYP	Theory65	Expt.†
Ethylene	0.3686	0.2649	0.2796	0.3140	0.376a	0.3859
Propene	0.3544	0.2503	0.2645	0.2994	0.354a	0.3565
1,3-Butadiene (E)	0.3188	0.2308	0.2444	0.2734	0.332a	0.3333
1,3-Pentadiene (E)	0.3053	0.2175	0.2309	0.2600	0.306b	0.3157
1,3,5-Hexatriene (E)	0.2862	0.2090	0.2218	0.2471	0.294c	0.3050
1,3-Cyclo-pentadiene	0.3051	0.2142	0.2283	0.2584	0.301d	0.3135
Benzene	0.3362	0.2491	0.2637	0.2920	0.329e	0.3399
Thiophene	0.3338	0.2400	0.2544	0.2856	0.327f	0.3271
Acetylene	0.4110	0.2928	0.3073	0.3473	0.409a	0.4189
Acetaldehyde	0.4628	0.3072	0.3193	0.3764	0.425a	0.3758

Equations

$$F_{\mu \nu }^{\mathrm{x}\sigma }=\sum _{\lambda \eta }{P}_{\lambda \eta }^{\sigma }(\mu \lambda |\eta \nu )=\sum _{\lambda \eta }{P}_{\lambda \eta }^{\sigma }\int \int \frac{{\chi }_{\mu ,\sigma }(\mathbf{r}){\chi }_{\lambda ,\sigma }(\mathbf{r}){\chi }_{\nu ,\sigma }({\mathbf{r}}^{\prime }){\chi }_{\eta ,\sigma }({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime }.$$

$$e_{\sigma }^{\mathrm{x}}(\mathbf{r})=-\sum _i^{occ}\sum _j^{occ}\int \frac{{\varphi }_{i,\sigma }^{\star }(\mathbf{r}){\varphi }_{j,\sigma }^{\star }(\mathbf{r}){\varphi }_{i,\sigma }({\mathbf{r}}^{\prime }){\varphi }_{j,\sigma }({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }$$

$$=-\sum _{\mu \nu }\sum _{\lambda \eta }{P}_{\mu \nu }^{\sigma }{P}_{\lambda \eta }^{\sigma }\int \frac{{\chi }_{\mu ,\sigma }(\mathbf{r}){\chi }_{\lambda ,\sigma }(\mathbf{r}){\chi }_{\nu ,\sigma }({\mathbf{r}}^{\prime }){\chi }_{\eta ,\sigma }({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }.$$

$$v_{\nu \eta }(\mathbf{r})=\int \frac{{\chi }_{\nu }({\mathbf{r}}^{\prime }){\chi }_{\eta }({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }=\sum _p\sum _q\int \frac{{\phi }_{\nu }^p({\mathbf{r}}^{\prime }){\phi }_{\eta }^q({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }.$$

$$E_{\sigma }^{\mathrm{x}}=\frac{1}{2}\int e_{\sigma }^{\mathrm{x}}(\mathbf{r})d\mathbf{r}.$$

$$e_{\sigma }^{\mathrm{x}}(\mathbf{r})=-\sum _{\nu }{Q}_{\nu }^{\sigma }(\mathbf{r})M_{\nu }^{\sigma }(\mathbf{r}),$$

$$Q_{\nu }^{\sigma }(\mathbf{r})=\sum _{\mu }{\chi }_{\mu ,\sigma }(\mathbf{r})P_{\mu \nu }^{\sigma },$$

$$M_{\nu }^{\sigma }(\mathbf{r})=\sum _{\eta }{Q}_{\eta }^{\sigma }(\mathbf{r})v_{\nu \eta }^{\sigma }(\mathbf{r}).$$

$$\frac{\partial E_{\sigma }^x}{\partial P_{\lambda \eta }^{\sigma }}=F_{\mu \nu }^{\mathrm{x}\sigma }=-{\int }_{\mathbf{r}}{\chi }_{\mu }(\mathbf{r})M_{\nu }^{\sigma }(\mathbf{r})d\mathbf{r}.$$

$$\bigg{[}-\frac{1}{2}\nabla^{2}+v^{\mathrm{p}}_{\mathrm{ion}}(\mathrm{\mathbf{r}})+v_{\mathrm{ext}}(\mathrm{\mathbf{r}})+v_{\mathrm{h}}[\rho(\mathrm{\mathbf{r}})]+v_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]\bigg{]}\psi_{i}^{\sigma}(\mathrm{\mathbf{r}})=\epsilon_{i}\psi_{i}^{\sigma}(\mathrm{\mathbf{r}}),$$

$$v_{\mathrm{ion}}^{\mathrm{p}}(\mathbf{r})=\sum _{{\mathbf{R}}_a}{v}_{\mathrm{ion},\mathrm{a}}^{\mathrm{p}}(\mathbf{r}-{\mathbf{R}}_a).$$

$$\sum _{\nu }{F}_{\mu \nu }^{\sigma }{C}_{\nu i}^{\sigma }={ϵ}_i^{\sigma }\sum _{\nu }^{\sigma }{S}_{\mu \nu }{C}_{\nu i}^{\sigma },$$

$$F_{\mu \nu }^{\sigma }=H_{\mu \nu }^{\mathrm{core}}+J_{\mu \nu }+F_{\mu \nu }^{\mathrm{xc}\sigma }.$$

$$r_i=r_0+(i-1)h_r,i=1,2,3,....,N_r,r_0=-\frac{N_r{h}_r}{2},r\in \{x,y,z\},$$

$$\rho ({\mathbf{r}}_g)=\sum _{\mu ,\nu }{P}_{\mu \nu }{\chi }_{\mu }({\mathbf{r}}_g){\chi }_{\nu }({\mathbf{r}}_g).$$

$$F_{\mu \nu }^{\mathrm{hxc}}=⟨{\chi }_{\mu }({\mathbf{r}}_g)|v_{\mathrm{hxc}}({\mathbf{r}}_g)|{\chi }_{\nu }({\mathbf{r}}_g)⟩=h_x{h}_y{h}_z\sum _g{\chi }_{\mu }({\mathbf{r}}_g)v_{\mathrm{hxc}}({\mathbf{r}}_g){\chi }_{\nu }({\mathbf{r}}_g).$$

$$E_{\sigma }^{\mathrm{x}}=\frac{1}{2}{h}_x{h}_y{h}_z\sum _g{e}_{\sigma }^{\mathrm{x}}({\mathbf{r}}_g),$$

$$F_{\mu \nu }^{\mathrm{x}\sigma }=-h_x{h}_y{h}_z\sum _g{\chi }_{\mu }({\mathbf{r}}_g)M_{\nu }^{\sigma }({\mathbf{r}}_g).$$

$$f(\mathbf{r})$$

$$F(\mathbf{k})$$

$$f\star g={\int }_{-\infty }^{\infty }f({\mathbf{r}}^{\prime })g(\mathbf{r}-{\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime }.$$

$$\mathcal{F}(f\star g)$$

$$f\star g$$

$$v_{\nu \eta }(\mathbf{r})=\int \frac{{\chi }_{\nu }({\mathbf{r}}^{\prime }){\chi }_{\eta }({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }=\int \frac{{\chi }_{\mu \nu }({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }={\chi }_{\nu \eta }(\mathbf{r})\star v^c(\mathbf{r}).$$

$$v_{\nu \eta }(\mathbf{r})={\mathcal{F}}^{-1}(v^c\{\mathbf{k}){\chi }_{\nu \eta }(\mathbf{k})\}\mathrm{where}{\chi }_{\nu \eta }(\mathbf{k})=\mathcal{F}\{{\chi }_{\nu \eta }(\mathbf{r})\}.$$

$$v^c(\mathbf{r})=\frac{\mathrm{erf}(\zeta \mathbf{r})}{\mathbf{r}}+\frac{\mathrm{erfc}(\zeta \mathbf{r})}{\mathbf{r}}$$

$$v^c({\mathbf{r}}_g)=v_{long}^c({\mathbf{r}}_g)+v_{short}^c({\mathbf{r}}_g).$$

$${\zeta }_{\mathrm{opt}}\equiv \underset{\zeta }{\mathrm{opt}}{E}_{\sigma }^{\mathrm{x}}\equiv \underset{\zeta }{\mathrm{opt}}{E}_{\mathrm{tot}}=\underset{N_x,N_y,N_z}{\mathrm{opt}}{E}_{\text{tot}},\mathrm{at}\mathrm{fixed}{h}_{\mathbf{r}},$$

$$E_{\mathrm{xc}}^{B3LYP}=(1-a_0)E_{\mathrm{LSDA}}^{\mathrm{x}}+a_0{E}_{\mathrm{HF}}^{\mathrm{x}}+a_{\mathrm{x}}{E}_{B88}^{\mathrm{x}}+a_{\mathrm{c}}{E}_{\mathrm{LYP}}^{\mathrm{c}}+(1-a_{\mathrm{c}})E_{\mathrm{VWN}}^{\mathrm{c}},$$

$$E_{\mathrm{xc}}^{PBE0}=b_0{E}_{\mathrm{HF}}^{\mathrm{x}}+(1-b_0)E_{\mathrm{PBE}}^{\mathrm{x}}+E_{\mathrm{PBE}}^{\mathrm{c}},$$

$$E_{\mathrm{xc}}^{\mathrm{BHLYP}}=(1-c_0)E_{\mathrm{LSDA}}^{\mathrm{x}}+c_0{E}_{\mathrm{HF}}^{\mathrm{x}}+E_{\mathrm{LYP}}^{\mathrm{c}}.$$

Full text

1 1. Introduction

Kohn-Sham density functional theory (KS-DFT) 1, 2 has been found to be an indispensable tool for determining structure and properties of many-electron systems like atoms, molecules, clusters, nano-materials and periodic systems, for past several decades3, 4, 5, 6. The KS mapping is exact, in principle, and has the ability to capture many-body correlation effects at a computational cost comparable to Hartree-Fock (HF) theory. But a major bottleneck lies in the exchange-correlation (XC) functional that uses electron density, $\rho(\mathrm{\mathbf{r}})$, to describe non-classical effects under the mean-field formalism of KS scheme. As its exact form is unknown as yet, it has to be approximated and accordingly, $E_{\mathrm{xc}}^{\mathrm{approx}}[\rho(\mathrm{\mathbf{r}})]=\int d^{3}\mathrm{\mathbf{r}}f^{\mathrm{approx}}_{\mathrm{xc}}(\rho,\nabla\rho,\tau,e^{\mathrm{x}},\dots)$. Commonly used functionals can be hierarchically characterized using these variables. Thus starting from a simpler local density approximation (LDA) (containing $\rho$ only), to general gradient approximation (GGA) (with addition of gradient of electron density, $\nabla\rho$), to meta-GGA (further addition of kinetic energy density, $\tau$), one gets a progressively involved and sophisticated range of functionals. With further inclusion of exact (also called HF) exchange energy, one gets the hybrid functionals, while incorporation of exact exchange energy density, $e^{\mathrm{x}}$, leads to a hyper-GGA type XC functionals, residing in the fourth rung of Jacob’s ladder. Now we also have functionals which go beyond the fourth rung (and include virtual orbitals), thus requiring higher computational cost, and have been successfully implemented in literature in the context of various physio-chemical properties7, 8.

In general the approximate local or semi-local functionals experience certain issues regarding (i) correct asymptotic behavior of exact exchange potential at long range and (ii) non-cancellation of spurious Coulomb self-repulsion energy–the so-called self-interaction error (SIE)9, 10. These two points are connected to each other and may be dealt with by bringing exact exchange in to the picture. This component may be combined with semi-local functionals, empirically11 or non-empirically 12, 13, thus improving the asymptotic nature, and consequently reducing the SIE to a greater extent. Accordingly, the global hybrid functionals 11 (B3LYP being a prototypical model) were proposed in the literature14, 15, 16, which improved the chemical accuracy in a large number of physico-chemical properties17, 18. On the other hand, non-dynamic (or static) correlation seems to be a serious problem for current DFT functionals5 which arises when an electronic state cannot be properly described by a single Slater determinant. It requires a linear combination of degenerate or nearly degenerate Slater determinants. Naturally it is very much demanding to model this correlation accurately within a single determinantal KS-DFT framework, and some methods have been proposed in the past decade under the name of hyper-GGA functionals19, 20, 21. More recently, general single-determinant model functionals have also been put forth to treat this correlation up to the dissociation limit of a co-valent bond22, 23. All these emerging functionals belonging to hyper-GGA category, require $e^{\mathrm{x}}$ as a fundamental variable and in general, their computational cost is much higher compared to the global hybrid functionals.

The HF exchange energy and corresponding matrix elements are usually evaluated through four-center electron repulsion integrals (ERI) analytically, within a Gaussian basis expansion of molecular orbitals (MO), providing following contributions to KS-Fock matrix,

[TABLE]

Here $(\mu\lambda|\eta\nu)$ is the four center ERI with $\mu,\nu,\lambda,\eta$ representing the atomic orbitals (AO) or contracted basis, and $P_{\lambda\eta}^{\sigma}$ is an element of the single-particle spin density matrix, $P^{\sigma}$ having spin $\sigma$. Throughout the past few decades, there has been extensive methodological developments to optimize its evaluation using a range of approaches varying in complexity, sophistication and accuracy. It is worthwhile mentioning some of the notable ones, such as pseudo-spectral scheme using Gaussian basis functions, multipole accelerated algorithm, LinK, truncated or short-range exchange kernels with/without the use of resolution-of-the-identity (RI) approximation, chain-of-sphere exchange based on quadrature, auxiliary density matrix, plane-wave basis functions, metropolis stochastic, pair-atomic RI, adaptively compressed exchange operator, density fitting with local domains, real-space with projection operators, rigorous integral screening, occupied-orbital fast multipole 24, 25, 26, 27, 28, 29, 30, 31, 32, etc. The HF exchange energy density can be represented as,

[TABLE]

The two expressions are defined in terms of KS occupied MO ($\phi_{\sigma}$) and AO ($\chi_{\sigma}$) respectively. Complex conjugate sign is omitted here since density matrix and basis are generally in real form. At first glance, it seems to be computationally more expensive than normal exchange energy calculation due to the fact that it needs to be evaluated at each grid point with four AO indices. But some recent developments proposed in the literature show substantial computational cost reduction via pair-atomic RI approximation33, 34, 35 or using a semi-numerical (SNR) scheme36, 37, 38, 39. Here we mainly focus on the latter implementation, where one of the integrations involved in ERI is carried out numerically on the effective spatial grid.

The main motivation of this work is thus to calculate HF exchange density, energy and its matrix elements which are essential ingredients for development of some of the new-generation (in particular, orbital-dependent) XC functionals. This is performed accurately in real-space Cartesian coordinate grid (CCG) following the procedure 39, where a key step is the identification and evaluation of an intermediate quantity, namely, two-center electrostatic potential (ESP), defined as,

[TABLE]

The first expression is defined in terms of AO, whereas the second one in terms of primitive basis functions, $\varphi_{\mu}^{p}$ for a given $\chi_{\mu}$ (for simplicity we remove the spin indices). Note that symbols $\mu,\nu$ signify contracted basis functions, whereas $p,q$ label primitive functions. We propose a direct, efficient numerical (NR) scheme based on Fourier convolution theorem (FCT) for accurate estimation of this ESP integral (first expression) using a range-separation (RS) technique, corresponding to long-and short-range interactions for Coulomb interaction kernel. A major concern is how one can define the optimal RS parameter for a successful mapping of this interaction kernel in CCG from first principles rather than empirically. Towards this direction, we have employed a recently developed grid optimization procedure 40, 41 with respect to total energy in non-uniform CCG (originally used for evaluation of classical Hartree potential) to optimize the RS parameter, to a very good level of accuracy, using a well-defined constraint, suggested in the literature. The successful implementation of this approach also paves the way for further development of the so-called range-separated hybrid functionals in connection with generalized KS theorem42. The underlying algorithm of the SNR method and subsequently a detailed account of our proposed NR scheme, denoted as RS-FCT, for estimation of this ESP integral, is provided in Sec. 2. Then its implementation is carried out using a pseudopotential KS-DFT in CCG, developed in our lab43, 44, 45, 40, 41. The relevant cost analysis between these two routes is also sketched upon. Section 3 offers the necessary computational and technical details. Finally, the feasibility, performance and accuracy of our results are critically assessed through quantities such as total energy and orbital energy, in Sec. 4, employing three-types of global hybrid functionals B3LYP”46, PBE012 and BHLYP11 along with the HF calculation for a decent (total 40) number of atoms and molecules. Some concluding remarks as well as the future prospect is given in Sec. 4.

2 2. Methodology

2.1 2.1. HF exchange energy, density and matrix elements

Let us begin with the SNR approach, where HF exchange energy is computed numerically by integrating the corresponding density, $e^{\mathrm{x}}_{\sigma}(\mathrm{\mathbf{r}})$ as given by,

[TABLE]

Now recasting Eq. (2) in the following form,

[TABLE]

one may envisage the construction of $e^{\mathrm{x}}_{\sigma}(\mathrm{\mathbf{r}})$ in three steps of comparable computational cost. The first quantity, $Q_{\nu}^{\sigma}(\mathrm{\mathbf{r}})$ may be represented as follows,

[TABLE]

in which the density matrix is combined with AOs through a simple matrix multiplication. The computational cost of this step scales as $\mathcal{O}(N_{g}N_{B}^{2})$, with $N_{g},N_{B}$ denoting total number of grid points, and number of AO basis. Next crucial (rate-determining) step is to evaluate the ESP integral (embedded in $M_{\nu}^{\sigma}(\mathrm{\mathbf{r}})$), which according to Eq. (3) (first expression), also scales as $\mathcal{O}(N_{g}N_{B}^{2})$. Alternately one may perform this integral analytically (corresponding to second expression of Eq. (3)) using primitive functions and standard recursion algorithm 47, 48. In that scenario, evidently each ESP integral scales as $\mathcal{O}(N_{g}N_{P}^{2})$, with $N_{P}$ identifying average number of primitive functions. The final step consists of computation of the quantity, $M_{\nu}^{\sigma}(\mathrm{\mathbf{r}})$ in accordance with the following expression,

[TABLE]

The scaling of this step is also same as that of ESP integral evaluation, but requires fewer steps than the latter; needing only one multiplication and one addition at innermost loop. This effectively provides an NR way to compute the HF exchange matrix, which according to Eq. (1), can be rewritten as,

[TABLE]

2.2 2.2. Pseudopotential KS-DFT in CCG

For a many-electron system, one can write the single-particle KS equation in presence of pseudopotential as (atomic unit employed unless stated otherwise),

[TABLE]

where $v^{\mathrm{p}}_{\mathrm{ion}}$ denotes the ionic pseudopotential, written as below,

[TABLE]

In the above equation, $v^{\mathrm{p}}_{\mathrm{ion},\mathrm{a}}$ signifies ion-core pseudopotential associated with atom A, situated at $\mathrm{\mathbf{R}}_{a}$. The classical Coulomb (Hartree) term, $v_{\mathrm{h}}[\rho(\mathrm{\mathbf{r}})]$ describes usual electrostatic interaction amongst valence electrons whereas $v_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]$ signifies the non-classical XC part of latter, and $\{\psi^{\sigma}_{i},\sigma=\alpha\quad\mathrm{or}\quad\beta\}$ corresponds to a set of $N$ occupied orthonormal spin molecular orbitals (spin-MOs). Within LCAO-MO approximation, the coefficients for expansion of spin-MOs satisfy a set of equations, very similar to that in HF theory,

[TABLE]

satisfying the orthonormality condition, $(\mathrm{\mathbf{C}}^{\sigma})^{\dagger}\mathrm{\mathbf{S}}\mathrm{\mathbf{C}}^{\sigma}=\mathrm{\mathbf{I}}$. Here $\mathrm{\mathbf{C}}^{\sigma}$ contains the respective spin-MO coefficients $\{C_{\nu i}^{\sigma}\}$ for a given spin-MO $\psi_{i}^{\sigma}(\mathrm{\mathbf{r}})$, $\mathrm{\mathbf{S}}$ is the usual overlap matrix corresponding to elements $S_{\mu\nu}$, $\bm{\epsilon}^{\sigma}$ refers to diagonal matrix of spin-MO eigenvalues $\{\epsilon_{i}^{\sigma}\}$. The KS-Fock matrix has elements $F_{\mu\nu}^{\sigma}$, constituting of following contributions,

[TABLE]

In this equation, all one-electron contributions, such as kinetic energy, nuclear-electron attraction and pseudopotential matrix elements are included in first term, whereas $J_{\mu\nu}$ and $F_{\mu\nu}^{\mathrm{xc}\sigma}$ account for classical Hartree and XC potentials respectively.

Now we discretize various quantities like localized basis function, electron density, MO as well as two-electron potentials directly on a 3D cubic box having $x,y,z$ axes,

[TABLE]

where $h_{r},N_{r}$ denote grid spacing and total number of points along each directions respectively. The electron density $\rho(\mathrm{\mathbf{r}})$ in this grid may be simply written as (“g” symbolizes discretized grid),

[TABLE]

At this stage, the two-electron contributions of KS matrix are computed through direct numerical integration in the grid,

[TABLE]

where $v_{hxc}$ refers to the combined Hartree and XC potential. The detailed construction of various potentials in CCG has been well documented in our earlier work43, 44, 45, 40, 41; hence not repeated here. Similarly, the HF exchange energy and its contributions towards KS-Fock matrix can be computed numerically in CCG as,

[TABLE]

2.3 2.3. RS-FCT for computation of ESP integral in CCG

Let us start with two functions $f(\mathrm{\mathbf{r}})$ and $F(\mathrm{\mathbf{k}})$ in $\mathrm{\mathbf{r}}$ and $\mathrm{\mathbf{k}}$ spaces, which are Fourier transforms (FT) (denoted by $\mathcal{F}$) of each other, as given below,

[TABLE]

With two functions $f(\mathrm{\mathbf{r}}),g(\mathrm{\mathbf{r}})$, one can form the convolution, defined by,

[TABLE]

The convolution theorem states that FT of convolution is just the product of individual FTs,

[TABLE]

One can make use of this theorem to rewrite the ESP integral as,

[TABLE]

The spin indices have been removed for simplicity. The last expression is in terms of convolution integral, where $\chi_{\nu\eta}$ denotes simple multiplication of two AO basis and $v^{c}(\mathrm{\mathbf{r}})$ represents Coulomb interaction kernel. Now one can invoke FCT to further simplify this integral,

[TABLE]

Here $v^{c}(\mathrm{\mathbf{k}})$ and $\chi(\mathrm{\mathbf{k}})$ stand for Fourier integrals of Coulomb kernel and AO basis respectively. The main concern lies in an accurate mapping of the former, which has singularity at $\mathrm{\mathbf{r}}\to 0$. In order to alleviate this problem, we apply a simple RS technique from the works of 49, 50, expanding the kernel into long- and short-range components with a suitably chosen RS parameter ($\zeta$). This gives rise to the following expression,

[TABLE]

In the above equation, $\mathrm{erf}(x)$ and $\mathrm{erfc}(x)$ denote error function and its complement respectively, while the second expression is written in real-space CCG. The Fourier integral of Coulomb kernel can be separated out as follows: (i) FT of short-range part is treated analytically and (ii) long-range portion is computed directly from FFT of corresponding real-space CCG values. Then the remaining problem lies in finding an optimum value of parameter $\zeta_{\mathrm{opt}}$ for successful mapping of Coulomb kernel in CCG from first principles. In this regard we propose a simple procedure which is found to be sufficiently accurate for all practical purposes (as exemplified in results that follows). This prompts us to write,

[TABLE]

using a suitably defined constraint 51 $(\zeta\times L=7)$, where $L\ (=N_{\mathrm{\mathbf{r}}}h_{\mathrm{\mathbf{r}}};\ \mathrm{\mathbf{r}}\in\{x,y,z\})$ is the smallest side of simulating box. This technique was already implemented in grid optimization procedure with respect to total energy in a non-uniform CCG 40, 41.

From the foregoing discussion it is clear that, each ESP integral involves only a combination of FFT (two forward and one backward transformation simultaneously) and hence scales as $\mathcal{O}(N_{g}\log N_{g})$. It is to be noted that this prescription gives a logarithmic scaling with grid size, in contrast to the quadratic scaling (with respect to primitive basis functions, $N_{P}$) in SNR scheme.

3 3. Computational details

The above described strategy for HF exchange component is implemented in case of three global hybrid functionals, namely, B3LYP, PBE0 and BHLYP, containing a variable amount of former with conventional DFT XC functional. As such, the XC energy corresponding to these three functionals are defined as follows,

[TABLE]

[TABLE]

and

[TABLE]

Following the work of Stephens et al., 52 values of $a_{0},a_{\mathrm{x}},a_{\mathrm{c}}$ are $0.2$, $0.72$, $0.81$ for B3LYP, whereas in case of PBE012, $b_{0}=0.25$. Note that the amount of HF exchange contribution in PBE0 is slighter higher than that of B3LYP, but even larger amount ($c_{0}=0.5$) is attributed in BHLYP. These are computed using KS orbitals obtained from solution of Eq. (9) in real-space CCG through the procedure delineated in previous subsections.

The present calculations employ following effective core potential (ECP) basis sets: SBKJC 54 for species containing Group-II elements and LANL2DZ 55 for Group-III or higher group elements. These are adopted from EMSL Basis Set Library56. All one-electron integrals are generated by standard recursion relations 47 using Cartesian Gaussian-type orbitals as primitive basis functions. The norm-conserving pseudopotential matrix elements in contracted basis are imported from GAMESS53 suite of program package. The relevant LDA- and GGA-type of functionals in connection with B3LYP, BHLYP and PBE0 are: (i) Vosko-Wilk-Nusair (VWN)–with the homogeneous electron gas correlation proposed in parametrization formula V 57 (ii) B88–incorporating Becke58 semi-local exchange (iii) Lee-Yang-Parr (LYP) 59 semi-local correlation (iv) Perdew-Burke-Ernzerhof (PBE)60 functional for semi-local exchange and correlation. All these are adopted from density functional repository program61 except LDA. The convergence criteria imposed in this communication are slightly tighter than our earlier work43, 44, 62, 45; this is to generate a more accurate orbital energies. Changes in following quantities were followed during the self-consistent field process, viz., (i) orbital energy difference between two successive iterations and (ii) absolute deviation in a density matrix element. They both were required to remain below a prescribed threshold set to $10^{-8}$ a.u.; this ensured that total energy maintained a convergence of at least this much, in general. In order to perform discrete Fourier transform, standard FFTW3 package63 was invoked. The resulting generalized matrix-eigenvalue problem is solved through standard LAPACK routine64 accurately and efficiently. All molecular calculations are performed in their experimental geometries, taken from NIST database65, excepting those in Table VI (see discussion later). Other details including scaling properties may be found in references 43, 44, 62, 45, 40, 41.

4 4. Result and Discussion

Before proceeding for main results, at first it may be appropriate to discuss the grid optimization (in CCG) in terms of the RS parameter, $\zeta_{\mathrm{opt}}$. This is illustrated through convergence of total energy, as delineated in Eq. (23). Table 1 shows this for a diatomic molecule HCl at its experimental geometry (bond length $1.275$Å). The HF energies are provided in non-uniform grid with respect to sparsity (regulated by $N_{x},N_{y},N_{z}$) for a fixed grid spacing (determined by $h_{r}$, chosen as 0.3), employing the NR and SNR schemes. The presentation strategy is similar to that in our previous works 40, 41; accordingly we first vary $N_{z}$, the number of grid points along inter-nuclear axis, keeping the same along $xy$ plane static at certain reasonable value, say $N_{x}=N_{y}=32$. As $N_{z}$ is gradually increased from 32 to 64 with an increment of 4, there is a smooth convergence in energy at around $N_{z}=56$ with a difference in total energy (we term it as grid accuracy) of about $10^{-6}$ a.u., between two successive steps. In the beginning, when $N_{z}$ goes through 40-44-48-52, one notices slow improvement in energy; after that it eventually attains the convergence for $N_{z}$ at around 60. Then in Set II in right-hand side, we vary $N_{x},N_{y}$ along $xy$ plane keeping $N_{z}$ fixed at its previously determined value of $60$. Now we can see that convergence in energy takes place at $N_{x}=N_{y}=50$ with same grid accuracy of Set I. Next, we have repeated the same calculation through the SNR scheme (ESP integral was approached analytically using recursion relations 47 containing expensive incomplete Gamma functions); the corresponding electronic energy value is supplied in footnote of this table. For sake of completeness, energy is also quoted from GAMESS 53 in footnote. It is quite encouraging to note that the two routes provide completely identical results as that in the reference. A similar exercise was undertaken for Cl2 ($R_{eq}=4.2$ a.u.), and the conclusion remains same, i.e., converged total energies in this case were found to be as follows: $-$29.34810 (NR), $-$29.34809 (SNR), $-$29.34813 (GAMESS) a.u. Thus it is evident that our direct, NR approach using RS-FCT in CCG, can correctly estimate the optimized RS-parameter, $\zeta_{\mathrm{opt}}$, and consequently, HF energy.

Next we move towards implementation of hybrid functional, B3LYP using above procedure. For this, we again consider a closed shell molecule, such as chlorine (Cl2) at its experimental bond length of $4.2$ a.u., along $z$ axis as a specimen case. The formal convergence and stability of solution is illustrated in Table 2 at a grid spacing of $h_{r}=0.3$. Same convergence criteria of previous table were imposed here as well. Once again the NR, SNR (provided in footnote) and GAMESS (also in footnote) energies excellently match with each other. An analogous study on HCl produced same result in energy, i.e., converged values in this case were $-$15.50844 (NR), $-$15.50843 (SNR), $-$15.50845 (GAMESS) a.u. This therefore establishes that our simple RS-FCT scheme in CCG can produce correct, meaningful results for XC functionals containing HF exchange.

In order to extend the scope and applicability of proposed scheme, we now report total energies for a set of 15 atoms and 15 molecules (along with test molecular systems) in Tables 3 and 4 respectively. To put things in perspective, we compare the NR and SNR results for three sets of calculations, namely HF, B3LYP and PBE0. Both tables engage identical basis set, pseudopotential and convergence (both grid and SCF) criteria. In all occasions, the maximum absolute deviation in energy (labeled $\mathrm{E}_{\mathrm{diff}}$) between NR and SNR, are produced side by side for easy comparison. The overall agreement between these two sets of results are in general excellent. A closer observation at Table 3 reveals that, the two energies are virtually indistinguishable for all the species with exception of O atom, where absolute deviation remains well below 0.0004 a.u. Similarly in Table 4 again, these two results coincide for 13 molecules; in this case the largest discrepancy (0.0001 a.u.) occurs for CH4. Note that, since the previous tables have established that NR and SNR energies are very close to reference GAMESS results (which have also been verified in these cases), we do not reproduce those energies any more.

In the last part of this study, in Table 5, we investigate the highest occupied molecular orbital (HOMO) energies of same set of atoms and molecules of Tables 3 and 4. As expected from previous tables, NR and SNR schemes produce practically identical results in these cases too. Hence now onwards, we present only NR results to save space. The available theoretical and experimental IPs are recorded for easy comparison. It is revealed that even though the theory does not account for electron correlation and relaxation effects, HF HOMO energies still compare very favorably with experiment than any of the three DFT functionals considered. This is apparent from the respective deviations; one notices underestimation of 14-41%, 9-37%, 1-28% for B3LYP, PBE0 and BHLYP functionals, whereas HF results show absolute deviation by 0-14%. However, it is interesting to note that, amongst the three functionals, the deviations fall from B3LYP to BHLYP, mainly because the fractional content of HF exchange (which has a pivotal role in determining correct asymptotic nature at long range) gradually increases in that order. This observation along with a consideration of total energies (from Tables 3 and 4), suggests that hybrid functionals provide a systematic improvement over the HF result. These observations are further complemented in Table 6, where a similar analysis is made for some $\pi$-electron molecules (simple, conjugated as well as aromatic). A very important aspect of such species is the so-called fundamental gap, which is the difference in energy between HOMO and LUMO. A satisfactory estimation of this gaps requires accurate knowledge of these two orbital energies. The current table is an attempt towards this direction. In this case, we have taken the calculated geometry (using B3LYP XC functional and cc-PVTZ basis set) from NIST 65 data base. As expected, we witness similar kind of pattern in calculated IP values as in previous table; with respect to experimental/theoretical values, the range of underestimation for hybrid functionals are: 18-32%, 15-28% and 11-19% respectively, for B3LYP, PBE0 and BHLYP.

A few remarks may now be made before passing. During the calculation, it is realized that the computational burden of our NR route is quite less (about eight times) than the SNR one. It occurs because of the obvious reason that the former route obviates the necessity of two primitive loops in Eq. (3); instead the ESP integral is directly calculated in CCG using contracted Gaussian functions. In order to make the SNR route attractive and fruitful for practical applications, an elegant solution was proposed 39 in the form of some optimal recurrence relations for analytical evaluation of ESP integral 39. We have not attempted to incorporate this procedure in this communication, as this is aside the main theme and purpose of this work. Here, we just wanted to establish the validity, efficacy and feasibility of a simple direct NR scheme in the broader context of certain DFT functionals for large-scale calculations. As demonstrated in this work, the NR route is quite efficient, and may be favorable, especially for larger basis set. The accuracy and ease of implementation augurs well for its further application in the development of range-separated hybrid functional from first principles. This leads to a more complete representations of HF exchange in asymptotic region which closes the gap between theoretical and experimental result, keeping the computational overhead at same level as in commonly used hybrid functionals66, and currently we are working on it.

5 4. Future and Outlook

We have demonstrated the feasibility and practicability of a direct NR scheme for computation of HF exchange energy density, energy and matrix in real-space CCG. This was applied for a host of atoms and molecules; properties such as total energy and orbital energies were offered, within a pseudopotential KS-DFT. Besides HF components, these were also presented for B3LYP, PBE0 and BHLYP hybrid functionals. The success of this approach relies on accurate estimation of ESP integral, which in turn depends on optimization of the RS parameter in Coulomb kernel. The obtained results for all these species from both NR and SNR scheme, are in excellent agreement with each other. A consideration of scaling properties suggests that this approach may be quite helpful in large-scale DFT calculations involving orbital-dependent functionals. Application to range-separated hybrid, hyper as well as local hybrid XC functionals would further enhance its scope in a wide range of systems. It may also be desirable to examine its performance in various important physio-chemical properties of many-electron systems.

6 Acknowledgement

AG is grateful to UGC for a Senior Research Fellowship (SRF). TM acknowledges IISER Kolkata for a Junior Research fellowship (JRF). Financial support from DST SERB, New Delhi, India (sanction order number EMR/2014/000838) is gratefully acknowledged.