DFT calculations of atoms and molecules in Cartesian grids

TL;DR

This paper develops a Cartesian grid-based variational DFT method for atoms and molecules, offering accurate results and a practical alternative to traditional atom-centered grid approaches, with initial steps towards non-uniform grids to reduce computational costs.

Contribution

The paper introduces a Cartesian grid-based variational DFT approach for atoms and molecules, including a Fourier convolution method for Coulomb potentials and initial non-uniform grid implementation.

Findings

Accurate calculation of energies and properties using Cartesian grids.

Effective Coulomb potential computation via Fourier convolution and Ewald decomposition.

Potential reduction in computational overhead with non-uniform grids.

Abstract

Density functional theory (DFT) has emerged as one of the most versatile and lucrative approaches in electronic structure calculations of many-electron systems in past four decades. Here we give an account of the development of a variational DFT method for atoms and molecules \emph{completely} in a Cartesian grid. The non-relativistic Kohn-Sham equation is solved by using an LCAO-MO ansatz. Atom-centered localized basis set, electron density, molecular orbitals, two-body potentials are directly constructed on the grid. We adopt a Fourier convolution method for classical Coulomb potentials by making an Ewald-type decomposition technique in terms of short- and long-range interactions. It produces quite accurate and competitive results for various properties of interest, such as component energy, total energy, ionization energy, potential energy curve, atomization energy, etc. Both local and non-local functionals are employed for pseudopotential as well as full calculations. While most results are offered in a \emph{uniform} grid, initial exploratory attempts are made in a \emph{non-uniform} grid, which can significantly reduce the computational overhead. This offers a practical, viable alternative to atom-centered grid-based implementations, currently exploited by the majority of programs available world-wide.

Tables

Table 1. Table 1: Energy components and total number of electrons, N , 𝑁 N, in various grids, for Cl 2 , at R = 4.2 𝑅 4.2 R=4.2 a.u., along with reference results. All quantities are in a.u. It is adopted from roy08 .

Set	A	B	C	D	E	F	G	H	Ref. schmidt93
$h_r$	0.3	0.4	0.2	0.3	0.4	0.1	0.2	0.1
$N_r$	32	32	64	64	64	128	128	256
$⟨T⟩$	11.00750	11.17919	11.18733	11.07195	11.06448	11.18701	11.07244	11.07244	11.07320
$⟨V_t^{\mathrm{ne}}⟩$	$-$83.43381	$-$83.68501	$-$83.70054	$-$83.45722	$-$83.44290	$-$83.69988	$-$83.45810	$-$83.45810	$-$83.45964
$⟨V_t^{\mathrm{ee}}⟩$	32.34338	31.22265	31.22885	31.00981	31.00306	31.22832	31.01000	31.01000	31.01078
$⟨E_{\mathrm{nu}}⟩$	11.66667	11.66667	11.66667	11.66667	11.66667	11.66667	11.66667	11.66667	11.66667
$⟨V⟩$	$-$39.42376	$-$40.79570	$-$40.80503	$-$40.78074	$-$40.77317	$-$40.80489	$-$40.78144	$-$40.78144	$-$40.78219
$⟨E_{\mathrm{el}}⟩$	$-$40.08293	$-$41.28318	$-$41.28437	$-$41.37545	$-$41.37535	$-$41.28455	$-$41.37566	$-$41.37566	$-$41.37566
$⟨E⟩$	$-$28.41626	$-$29.61651	$-$29.61770	$-$29.70878	$-$29.70868	$-$29.61789	$-$29.70900	$-$29.70900	$-$29.70899111Refers to the grid-DFT calculation; corresponding grid-free DFT value is $-$29.71530 a.u.
$N$	13.89834	13.99939	13.99865	14.00002	14.00003	13.99864	14.00000	13.99999	13.99998

Table 2. Table 2: Variation of the energy components and N 𝑁 N with respect to the grid parameters for Cl 2 and HCl with reference values. BLYP results in a.u.

	Cl2 ($R=4.2$ a.u.)			HCl ($R=2.4$ a.u.)
Set	A	B	Ref. schmidt93	A	B	Ref. schmidt93
$N_r$	64	128		64	128
$h_r$	0.3	0.2		0.3	0.2
$⟨T⟩$	11.21504	11.21577	11.21570	6.25431	6.25464	6.25458
$⟨V_t^{\mathrm{ne}}⟩$	$-$83.72582	$-$83.72695	$-$83.72685	$-$37.29933	$-$37.29987	$-$37.29979
$⟨V_t^{\mathrm{ee}}⟩$	31.07572	31.07594	31.07594	12.63884	12.63903	12.63901
$⟨E_{\mathrm{nu}}⟩$	11.66667	11.66667	11.66667	2.91667	2.91667	2.91667
$⟨V⟩$	$-$40.98344	$-$40.98434	$-$40.98424	$-$21.74382	$-$21.74417	$-$21.74411
$⟨E_{\mathrm{el}}⟩$	$-$41.43506	$-$41.43524	$-$41.43522	$-$18.40618	$-$18.40620	$-$18.40620
$⟨E⟩$	$-$29.76840	$-$29.76857	$-$29.76855111Corresponding grid-free DFT value is $-$29.74755 a.u. schmidt93 .	$-$15.48951	$-$15.48953	$-$15.48953222Corresponding grid-free DFT value is $-$15.48083 a.u. schmidt93 .
$N$	14.00006	14.00000	13.99998	8.00002	8.00000	8.00000

Table 3. Table 3: Comparison of calculated negative eigenvalues of Cl 2 and HCl with reference values. BLYP results are given in a.u. See text for more details.

MO	Cl2 ($R=4.2$ a.u.)			MO	HCl ($R=2.4$ a.u.)
Set	A	B	Ref. schmidt93		A	B	Ref. schmidt93
$N_r$	64	128			64	128
$h_r$	0.3	0.2			0.3	0.2
$2{\sigma }_g$	0.8143	0.8143	0.8143	$2\sigma$	0.7707	0.7707	0.7707
$2{\sigma }_u$	0.7094	0.7094	0.7094	$3\sigma$	0.4168	0.4167	0.4167
$3{\sigma }_g$	0.4170	0.4171	0.4171	$1{\pi }_x$	0.2786	0.2786	0.2786
$1{\pi }_{xu}$	0.3405	0.3405	0.3405	$1{\pi }_y$	0.2786	0.2786	0.2786
$1{\pi }_{yu}$	0.3405	0.3405	0.3405
$1{\pi }_{xg}$	0.2778	0.2778	0.2778
$1{\pi }_{yg}$	0.2778	0.2778	0.2778

Table 4. Table 4: Comparison of various energies and N 𝑁 N for several atoms and molecules for different XC functionals. For more details, see text.

	$⟨T⟩$				$-⟨V⟩$
System	LDA	Ref. (schmidt93, )	BLYP	Ref. (schmidt93, )		LDA	Ref. (schmidt93, )	BLYP	Ref.  (schmidt93, )
P	2.35430	2.35334	2.38891	2.38890		8.73501	8.73404	8.78249	8.78248
H2S	4.90204	4.90197	4.98071	4.98066		16.10707	16.10698	16.21919	16.21913
MgCl2	11.62114	11.62208	11.75947	11.75999		42.34513	42.34621	42.54049	42.54103
SiH2Cl2	13.95036	13.94989	14.14948	14.14945		48.78729	48.78685	49.04463	49.04461
	$-⟨E⟩$				N
System	LDA	Ref. (schmidt93, )	BLYP	Ref. (schmidt93, )		LDA	Ref. (schmidt93, )	BLYP	Ref. (schmidt93, )
P	6.38070	6.38071	6.39358	6.39358		5.00000	4.99999	4.99999	4.99999
H2S	11.20503	11.20501	11.23848	11.23846		8.00000	7.99999	8.00000	7.99999
MgCl2	30.72399	30.72413	30.78102	30.78104		16.00004	15.99957	16.00004	15.99999
SiH2Cl2	34.83693	34.83696	34.89515	34.89516		19.99999	20.00015	19.99999	20.00000

Table 5. Table 5: Comparison of kinetic ( ⟨ T ⟩ delimited-⟨⟩ 𝑇 \left\langle T\right\rangle ), potential ( ⟨ V ⟩ delimited-⟨⟩ 𝑉 \left\langle V\right\rangle ), total (E) energies and N 𝑁 N for several atoms and molecules for different XC functionals. See text for details.

	$⟨T⟩$				$-⟨V⟩$
System	LDA	BLYP	PBE	FT97		LDA	BLYP	PBE	FT97
Br	4.22038	4.27022	4.17153	4.22828		17.28157	17.36122	17.29168	17.37615
Al2	1.35409	1.37642	1.32292	1.34896		5.21624	5.24844	5.20904	5.25881
SiH4	3.65832	3.76989	3.53302	3.60707		9.82333	9.97905	9.75741	9.86882
S2	7.58306	7.68138	7.51586	7.59269		27.62272	27.77409	27.64063	27.76097
	$-⟨E⟩$				$N$
System	LDA	BLYP	PBE	FT97		LDA	BLYP	PBE	FT97
Br	13.06119	13.09100	13.12015	13.14787		7.00000	7.00000	6.99999	6.99999
Al2	3.86215	3.87203	3.88612	3.90985		5.99999	5.99999	5.99999	5.99999
SiH4	6.16500	6.20916	6.22438	6.26175		7.99999	8.00000	7.99999	7.99999
S2	20.03966	20.09271	20.12477	20.16828		12.00000	12.00000	12.00000	12.00000

Table 6. Table 6: Comparison of − ϵ HOMO subscript italic-ϵ HOMO -\epsilon_{\mathrm{HOMO}} (in a.u.) and atomization energies (in kcal/mole) for several molecules for different XC functionals with literature experimental data.

	$-{ϵ}_{\mathrm{HOMO}}$ (a.u.)						Atomization energies (kcal/mol)
System	LDA	BLYP	PBE	FT97	LBVWN	Expt. (afeefy05, )	LDA	BLYP	Expt. (afeefy05, )
SiH4	0.3188	0.3156	0.2919	0.2702	0.4624	0.4042	339.43	312.02	302.6
PH2	0.2170	0.2111	0.1970	0.1741	0.3504	0.3610	152.77	139.92	149.2
P4	0.2712	0.2575	0.2525	0.2369	0.3964	0.3432	200.77	142.99	285.9
Al2	0.1407	0.1400	0.1318	0.1161	0.2371	0.1984	22.92	21.42	37.0

Table 7. Table 7: Energy components and orbital energies for two atoms (N, O) using Cartesian grid. All-electron calculations with LDA XC functionals using STO-3G basis set are given along with reference data. PR implies Present Result.

Energy	N		O		Orbital	N		O
	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )		PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )
$⟨T⟩$	53.66407	53.66407	73.44497	73.44497	${ϵ}_{1s}^{\alpha }$	$-$13.6312	$-$13.6311	$-$18.3331	$-$18.1330
$⟨V^{\mathrm{ne}}⟩$	$-$127.32649	$-$127.32649	$-$176.32432	$-$176.32432	${ϵ}_{2s}^{\alpha }$	$-$0.6152	$-$0.6153	$-$0.7538	$-$0.7537
$⟨V^{\mathrm{ee}}⟩$	20.25538	20.25536	29.42801	29.42799	${ϵ}_{2p_x}^{\alpha }$	$-$0.1671	$-$0.1672	$-$0.1941	$-$0.1942
$⟨V⟩$	$-$107.07111	$-$107.07113	$-$146.89631	$-$146.89634	${ϵ}_{2p_y}^{\alpha }$	$-$0.1671	$-$0.1672	$-$0.1941	$-$0.1942
$⟨E⟩$	$-$53.40704	$-$53.40701	$-$73.45134	$-$73.45137	${ϵ}_{2p_z}^{\alpha }$	$-$0.1671	$-$0.1672	$-$0.1085	$-$0.1085
$N$	6.99999	6.99999	7.99999	7.99999	${ϵ}_{1s}^{\beta }$	$-$13.5837	$-$13.5836	$-$18.2972	$-$18.2971
					${ϵ}_{2s}^{\beta }$	$-$0.4449	$-$0.4450	$-$0.6301	$-$0.6302
					${ϵ}_{2p_x}^{\beta }$			$-$0.0378	$-$0.0379

Table 8. Table 8: Electronic energy a of Cl 2 in non-uniform grid ( h r = 0.3 subscript ℎ 𝑟 0.3 h_{r}=0.3 ). LDA results in a.u.

Set I					Set II
$N_x$	$N_y$	$N_z$	$⟨E⟩$		$N_x$	$N_y$	$N_z$	$⟨E⟩$
30	30	30	$-$25.32542		30	30	70	$-$29.66504
-	-	40	$-$29.48297		32	32	-	$-$29.68722
-	-	50	$-$29.65381		34	34	-	$-$29.69838
-	-	60	$-$29.66461		36	36	-	$-$29.70387
-	-	70	$-$29.66504		38	38	-	$-$29.70651
-	-	80	$-$29.66505		40	40	-	$-$29.70777
-	-	90	$-$29.66505		42	42	-	$-$29.70837
36	36	40	$-$29.52104		44	44	-	$-$29.70865
-	-	60	$-$29.70343		46	46	-	$-$29.70877
-	-	70	$-$29.70387		48	48	-	$-$29.70884
-	-	80	$-$29.70388		50	50	-	$-$29.70886
-	-	90	$-$29.70388		60	60	-	$-$29.70888

Equations

$${\rho }_1({\mathbf{x}}_1^{{}^{\prime }}|{\mathbf{x}}_1)=N\int \cdots \int \psi ({\mathbf{x}}_1^{{}^{\prime }},{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N){\psi }^{\star }({\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N)d{\mathbf{x}}_2\cdots d{\mathbf{x}}_N$$

$${\tau }_2({\mathbf{x}}_1^{{}^{\prime }},{\mathbf{x}}_2^{{}^{\prime }}|{\mathbf{x}}_1,{\mathbf{x}}_2)=\frac{N(N-1)}{2}\int \cdots \int \psi ({\mathbf{x}}_1^{{}^{\prime }},{\mathbf{x}}_2^{{}^{\prime }},{\mathbf{x}}_3,\cdots ,{\mathbf{x}}_N){\psi }^{\star }({\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,\cdots ,{\mathbf{x}}_N)d{\mathbf{x}}_3\cdots d{\mathbf{x}}_N$$

$$\rho (\mathbf{r})=N\int \cdots \int \psi (\mathbf{r}\sigma ,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N){\psi }^{\star }(\mathbf{r}\sigma ,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N)d\sigma d{\mathbf{x}}_2\cdots d{\mathbf{x}}_N.$$

$${\mu }_{\mathrm{TF}}=\frac{\delta E_{\mathrm{TF}}[\rho (\mathbf{r})]}{\delta \rho (\mathbf{r})}=\frac{5}{3}{C}_{\mathrm{F}}\rho (\mathbf{r}{)}^{\frac{2}{3}}+\int d{\mathbf{r}}^{\prime }\frac{\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+v_{\mathrm{ext}}(\mathbf{r}),$$

$$E_{\mathrm{TFD}}[\rho (\mathbf{r})]=C_{\mathrm{F}}\int \rho (\mathbf{r}{)}^{\frac{5}{3}}d\mathbf{r}+\int \rho (\mathbf{r})v_{\mathrm{ext}}(\mathbf{r})d\mathbf{r}-C_{\mathrm{x}}\int \rho (\mathbf{r}{)}^{\frac{4}{3}}d\mathbf{r}+\frac{1}{2}\int \int \frac{\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime },$$

$$\delta\bigg{\{}E[\rho(\mathrm{\mathbf{r}})]-\mu\left[\int\rho(\mathrm{\mathbf{r}})d\mathrm{\mathbf{r}}-N\right]\bigg{\}}=0,$$

$$E[\rho ]$$

$$F[\rho ]$$

$$E[\rho (\mathbf{r})]\ge E[{\rho }_0(\mathbf{r})],$$

$$\frac{\partial \rho (\mathbf{r},t)}{\partial t}+\nabla \cdot \mathbf{j}(\mathbf{r},t)=0,$$

$$\frac{\partial \chi (\mathbf{r},t)}{\partial t}+\frac{1}{2}(\nabla \chi {)}^2+\frac{\delta \mathrm{G}[\rho ]}{\delta \rho }+\frac{\delta {\mathrm{E}}_{\mathrm{el}-\mathrm{el}}[\rho ]}{\delta \rho }+v_{\mathrm{ext}}(\mathbf{r},t)=0,$$

$$\Psi (\mathbf{r},t)=\rho (\mathbf{r},t{)}^{1/2}{e}^{i\chi (\mathbf{r},t)},$$

$$\left[-\frac{1}{2}{\nabla }^2+v_{\mathrm{eff}}([\rho ];\mathbf{r},t)\right]\Psi (\mathbf{r},t)=i\frac{\partial \Psi (\mathbf{r},t)}{\partial t},$$

$$H_{\mathrm{s}}=-\frac{1}{2}\sum _i^N{\nabla }_i^2+\sum _i^N{v}_{\mathrm{s}}({\mathbf{r}}_i).$$

$${\hat{f}}^{\mathrm{KS}}{\psi }_i={ϵ}_i{\psi }_i,$$

$${\hat{f}}^{\mathrm{KS}}=-\frac{1}{2}{\nabla }^2+v_{\mathrm{s}}(\mathbf{r}).$$

$${\rho }_{\mathrm{s}}(\mathbf{r})=\sum _i^N\sum _{\sigma }|{\psi }_i(\mathbf{r},\sigma ){|}^2={\rho }_0(\mathbf{r}).$$

$$T_{\mathrm{s}}=-\frac{1}{2}\sum _i^N{\nabla }_i^2.$$

$$F[\rho ]$$

$$E_{\mathrm{xc}}[\rho ]$$

$$E[\rho ]=\int v_{\mathrm{ext}}(\mathbf{r})\rho (\mathbf{r})d\mathbf{r}+J[\rho ]+T_{\mathrm{s}}[\rho ]+E_{\mathrm{xc}}[\rho ]$$

$$\bigg{[}-\frac{1}{2}\nabla^{2}+v_{\mathrm{eff}}(\mathrm{\mathbf{r}})\bigg{]}\psi_{i}(\mathrm{\mathbf{r}})=\epsilon_{i}\psi_{i}(\mathrm{\mathbf{r}})$$

$$v_{\mathrm{eff}}(\mathbf{r})=v_{\mathrm{ext}}(\mathbf{r})+\int \frac{\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }+v_{\mathrm{xc}}(\mathbf{r}).$$

$$I=\int F(\mathbf{r})d\mathbf{r}=\sum _j{I}_j=\sum _j\int F_j(\mathbf{r})d\mathbf{r},$$

$$F_j(\mathbf{r})=w_j(\mathbf{r})F(\mathbf{r})$$

$$\sum _j{w}_j(\mathbf{r})=1\forall \mathbf{r},$$

$${\mu }_{ij}=\frac{{\mathbf{r}}_i-{\mathbf{r}}_j}{R_{ij}}$$

$$w_i(\mathbf{r})=\frac{W_i(\mathbf{r})}{{\sum }_j{W}_j(\mathbf{r})};W_i(\mathbf{r})=\prod _{j\ne i}s({\mu }_{ij}).$$

$$s_k(\mu )=\frac{1}{2}\left[1-f_k(\mu )\right],$$

$$I_j={\int }_o^{\infty }{\int }_0^{\pi }{\int }_0^{2\pi }{F}_j(r_j,{\theta }_j,{\varphi }_j)r_j^{2}sin{\theta }_{j}d{r}_{j}d{\theta }_{j}d{\varphi }_j.$$

$$F_j(r_j)={\int }_0^{\pi }{\int }_0^{2\pi }{F}_j(r_j,{\theta }_j,{\varphi }_j)sin{\theta }_{j}d{\theta }_{j}d{\varphi }_j\approx \sum _{i=1}^{n_j^{\Omega }}{\omega }_i^{\Omega }{F}_j(r_j,{\Omega }_i),$$

Full text

Keywords: Density functional theory, many-electron system, exchange-correlation functional, Cartesian grid.

I Introduction

The past several decades have witnessed a proliferation of ab initio quantum mechanical methods for elucidation of electronic structure, properties and dynamics of many-electron systems such as atoms, molecules, clusters, solids. With rapid advances in mathematical formalisms and smart numerical algorithms, as well as computer architecture and resources, lately, some of these have been extended to the case of materials science, nano-science and biological science, etc., with unprecedented success. The electronic Schrödinger equation (SE) of such an $N$-electron system is essentially a many-body problem involving space, spin and time coordinates. This usually causes enormous complexity and hence, exact analytical solutions are unavailable in vast majority cases, leaving aside a few countable idealized model problems which offer such solutions. Thus almost invariably recourse must be taken to approximate methods. So a well-defined theoretical framework of approximation is indeed necessary to explore the electronic problem in a systematic way as accurate as possible. Nowadays, it is not a mere contemplation to obtain an optimal solution for a system containing 100 or more electrons, with physically and chemically meaningful accuracy. However, as the number of electrons increases, computational complexity grows up exponentially, which requires new paradigms of development regarding basic theory, methodology, mathematical algorithm, sophisticated computer code and practical computational implementation. The ultimate objective is to probe bigger and bigger systems with as much accuracy and efficiency, as possible. Significant strides have been made in last few decades in these aforementioned areas, and in the coming years, a lot of progress will be continued in order to reach the eventual goal of theory predicting experimental results in an explanatory level relevant to practical applications.

Besides ground-state electronic structure, a lot of important aspects are involved in atomic physics to explain the experimental observations in a satisfactory level. Some of these, e.g., regard excited states of discrete and continuum spectra, ionization dynamics of Rydberg states, spectroscopic levels, oscillator strengths for electronic transitions, fine-structures due to spin-orbit coupling, external perturbations due to static probes, multi-photon ionizations, above-threshold ionizations, electron scattering, relativistic corrections to the Hamiltonian and so on. But in case of molecules, calculation of physical and chemical properties are somehow less straightforward and clouded, simply because of the presence of more than one nucleus. At a first glance, it requires a rigorous understanding of the molecular potential energy surface through construction of the Hessian matrix, and all other molecular properties are inherently dependent on it. Such a surface consists of a large number of maxima and minima at unknown locations, and the global minimum corresponds to the ground-state structure, whereas the paths between the minima are the primary focus of chemical reactions. It is worthwhile to mention a few relevant molecular properties, such as dipole and multipole moments, vibrational frequencies, molecular dissociation energies, linear and non-linear responses to external static probes, molecular scattering through electron beams, molecular reaction dynamics, nonlinear reaction dynamics, photochemistry, etc. A significant amount of interest concerns with modern spectroscopic techniques, namely, electronic spectroscopy, photo-electron spectroscopy, ro-vibrational spectroscopy, nuclear magnetic resonance, nuclear quadrupole resonance, electron spin resonance, Mössbauer spectroscopy etc.

An integral consideration in modern quantum chemistry is the obligatory computational work, which has profound impact on future theoretical development. In order to perform the complex mathematical operations for larger systems within a suitable framework, it is imperative to have a systematic way of computation, and sufficient intricacy to capture the characteristic subtle interactions. Sometimes the computational output may be very crude (overestimating or underestimating), and other times it could be within the bound accuracy, but that is mostly dictated by the nature of the practical problem at hand. Now, in order for a given computational method to be successful and viable, the two deciding factors are its accuracy and cost effectiveness. In today’s computational chemistry, it is desirable to achieve the energy of a chemical reaction within the bounds of chemical accuracy, which is less than 1 kcal/mole. And the main challenge is to lower computational cost without much compromise on accuracy for increasingly bigger and complicated systems. Starting from relatively less accurate and dependable empirical or semi-empirical schemes to modern ab initio methods, all fall under the purview of current computational quantum chemistry repertoire. The latter covers a large array of methods including Hartree-Fock (HF), post-HF, multi-reference, quantum Monte-Carlo, density functional theory (DFT), etc., each providing a different computational cost depending on the particular problem. In order to broaden the domain of feasibility to larger systems, significant developments have been made towards the so-called linear-scaling approaches, which slash computational scaling very swiftly with molecular size. In this context, it is appropriate to mention some of the most popular and prominent electronic structure packages, which are routinely used by the scientific community, such as Gaussian, GAMESS, ACES, CADPAC, Dalton, Hyperchem, Jaguar, Molcas, Molpro, NWChem, PQS, Psi4, Q-Chem, Spartan, TURBOMOLE, UTChem, etc. On the other hand, one could take help of molecular mechanics which uses classical mechanics to model a molecular, as well as large biological system or material cluster, having much lower computational cost, but is also less sophisticated and reliable. These are usually applied in the field of molecular dynamics (in the context of thermodynamic and kinetic properties of a system) and potential energy minimization (e.g., to calculate binding constants, protein folding kinetics, protonation equilibria, active site coordinates). In recent years, computational chemistry has also been used vigorously in medicinal chemistry, leading to the discovery of new life-saving drugs. Essentially this helps to model the synthetic compounds in a systematic way, providing some critical checks on huge labor and chemical cost. A large number of beautiful and in-depth books, reviews and articles are available on the subject. Some of them are referred here (yarkony95, ; szabo96, ; simon97, ; kohn99, ; helgaker00, ; springborg00, ; young01, ; foulkes01, ; lewars03, ; cramer04, ; martin04, ; hoffman07, ; jensen07, ; kaisas07, ; borden11, ; springborg12a, ; springborg14, ).

In order to solve the electronic SE for an interacting system, two distinct well-established theoretical formalisms have gained popularity and credence in the last few decades, with contributions from physicists and chemists alike. We begin our discussion with a comparatively simple HF theory, which in a sense, is founded on a mean-field approximation. In essence, the instantaneous inter-electronic interaction is accounted for, in an average way, where each electron sees an average field because of the presence of the other electrons, leading to the optimization of a single Slater determinantal wave function consisting of individual spin orbitals. Though it offers a crucial first step towards the problem, its application is limited due to a complete neglect of electron correlation effects. Post HF formalisms are mainly concerned on how this subtle effects could be incorporated in an efficient manner. This family includes a variety of versatile formalisms, such as Møller-Plesset perturbation theory (MPn) (saebo93, ), configuration interaction (CI) (szabo96, ), coupled-cluster ansatz (they come in many flavours like CCSD, CCSD(T), CCSDT, CCSDTQ, etc.) (bartlett89, ), multi-reference perturbation theory (such as CASPT2) (roos96, ), etc. These methods offer potentially authentic and reliable results, but are quite difficult to be implemented computationally for large $N$, primarily because of their unfavorable scaling. They differ in scaling, cost, efficiency and accuracy; naturally a major research effort for such methods is directed towards betterment of these issues. They have their own space of applicability depending upon a given problem and often provide a benchmark of accuracy for the system under investigation.

The second approach relies on utilizing the limited information coming from electron density, density matrix or Green’s function in an optimized way with the help of the variational principle. Amongst them, however, DFT has appeared as the most versatile and outstanding method in electronic structure calculation. Many beautiful books and reviews are available on the subject; some recent ones are (parr89, ; chong95, ; seminario96, ; joubert98, ; dobson98, ; koch01, ; young01, ; fiolhais03, ; gidopoulos03, ; hu08, ; cramer09, ; cohen11, ; burke12, ; roy12, ; becke14, ; jones15, ). The basic variable, the single-particle electron density, $\rho(\mathrm{\mathbf{r}})$, is a real-valued and experimentally realizable quantity, in contrast to the traditional complex-valued $N$-electron wave function. Furthermore, the most attractive part is that it incorporates the complicated electron correlation effects in a transparent and elegant way. Being a single determinant approach, it is computationally more tractable; the burden remains at same level as HF. However, this apparent simplicity and novelty is somehow tainted and compensated by the hitherto unavailable, all-important, universal, exchange-correlation (XC) density functional. Although uncovering its exact form constitutes one of the most active and fertile area of research, unfortunately it still remains elusive; hence it must be approximated. Success of DFT lies heavily on the choice of these approximate XC functional. Historically, the journey started with the so-called local density approximation (LDA) in 1930 by Dirac (dirac30, ) using the form of exchange assuming a uniform electron gas model. Later on, a working LDA correlation part was developed (vosko80, ; perdew92, ). However, the LDA XC functional suffers from serious problems in atoms, molecules and solids, as the density of such systems are far apart from the uniform electron gas. Since it is essential to have a good knowledge of the gradient of the electron density, a whole family of so-called generalized gradient approximation (GGA)-based functionals were proposed over the years. A few prominent ones are: B88 (becke88a, ), LYP (lee88, ), PBE (perdew96, ) etc. In the next step, higher-order derivatives of the density including the kinetic energy density lead to the development of meta-GGA kind of functionals, such as B88C (becke88, ), Becke-Rousse (becke89, ), TPSS (tao03, ) etc. Next crucial step lead to hybrid functionals, which bring exact HF exchange into the picture. Depending on the co-factor of different kinds of XC functionals like LDA, GGA with latter, a variety of hybrid functionals have been published, which include one of the most popular and versatile candidates in quantum chemistry, namely, B3LYP (becke88a, ; lee88, ). The LDA functional typically does not have an impressive energetic performance, but interestingly it executes better performance than GGA functional regarding geometry optimization. Generally it is true that hybrid functionals work better than LDA due to incorporation of HF exchange. But in case of transition-metal chemistry and response properties, hybrid functionals show rather quite poor performance. Also, a majority of these approximate functionals has some sort of deficiency regarding activation energies of chemical reactions, which is somehow circumvented in the recently developed range-separated, density fitting, adiabatic-connection, local hybrid functionals or those involving unoccupied orbitals and eigenvalues, etc. These are found to be largely good for dispersion and van der Waals interactions. A vast amount of literature exists on the topic; the interested reader may look at the lucid reviews cohen11 ; becke14 ; jones15 and references therein.

Triumphs of DFT is now well-established in quantum chemistry and condensed matter physics; so much so that nowadays it is the most dominant and visible workforce. For large scale computation of materials, two distinct numerical approaches are engaged to solve the relevant many-electron equation. The condensed matter community chooses plane-wave basis sets, where periodic boundary conditions are appropriate. On the other hand, chemists prefer a localized atom-centred Gaussian basis set, because of being fondly attached with molecules (non-periodic systems), where the problematic multi-center integrals can be evaluated analytically. With advances in computational facilities, modern DFT can provide results with benchmark accuracy accompanied with experimental results at least for small molecules in gaseous phase. Throughout the past several decades, a wide variety of noteworthy applications were made regarding molecular properties (including structures, thermo-chemistry, various spectroscopic quantities, responses to external perturbations), bulk and surface properties of solids, band-gaps and optical properties of solids, interactions of small molecules with surfaces with focus on structures, binding energies and catalytic chemistry, modelling of photochemical reactions, etc. Of late, extensive applications were reported in some other areas, such as nano-technology, semiconductor quantum dots, quantum confinement imposed by various kinds of nano structures. Practical applications from biological fields include modelling of enzymatic catalysis and active-site chemistry, cooperative activity in backbone of hydrogen bonding, modelling of beta-sheet formation, enzyme functioning etc.

The present chapter gives an account of the work done in our laboratory for static, non-relativistic ground states in many-electron atoms/molecules. Within the Born-Oppenheimer approximation and the Hohenberg-Kohn-Sham framework, an implementation of DFT in a Cartesian Coordinate grid (CCG) is offered. By using a linear combination of Gaussian functions, molecular orbitals (MO) and quantities like basis functions, electron density as well as classical Hartree and non-classical XC potentials are constructed on a 3D real CCG directly roy08 ; roy08a ; roy09 ; roy10 ; roy11 . No additional auxiliary basis is used for the charge density. A Fourier convolution method, involving a combination of Fast Fourier transform (FFT) and inverse FFT martyna99 ; minary02 is used to obtain the Coulomb potential quite accurately. Analytical one-electron Hay-Wadt-type effective core potentials wadt85 ; hay85 , made of sums of Gaussian type functions, are used to represent the inner core electrons whereas energy-optimized truncated Gaussian bases are used for valence electrons. Viability and suitability of this simpler grid is compared and contrasted with a routinely used atom-centered grid (ACG) for a decent number of atoms, molecules. This is demonstrated by presenting the total energy, energy components, orbital energy, potential energy curve, atomization energy for a bunch of local and non-local XC functionals.

II The methodology

II.1 Electron density as basic variable

In the usual wave function-based approach, the difficulty to solve the many-electron SE increases rapidly with the number of electrons present in the system. In order to reduce the dimensionality of the problem, one can introduce the concept of reduced density matrices (lowdin95, ). One- and two-particle reduced density matrices, which are adequate for the evaluation of expectation values of one- and two-particle operators, are defined as,

[TABLE]

[TABLE]

where each of ($\{\mathrm{\mathbf{x}}_{i}\};\ i=1,2,3,\cdots,N$) consists of both spatial coordinates $\{\mathrm{\mathbf{r}}_{i}\}$ as well as spin coordinates ($\sigma_{i}$). Now, the spin-less, single-particle density $\rho(\mathrm{\mathbf{r}})$, which is the diagonal element of one-particle density matrix, is given as,

[TABLE]

In case of fermions, it actually describes a three-dimensional distribution of electrons in a system. In sharp contrast to the complex-valued wave function, $\rho(\mathrm{\mathbf{r}})$ is a function of only three coordinates, irrespective of the number of electrons present in it. As it is a fundamental physical observable and can be determined experimentally, this enables one to test directly the veracity of chemical calculations and approximations. The beauty of DFT lies in formulating a many-particle problem within a single-particle framework with the aid of density being considered as a central quantity.

The first seminal work in this direction was provided independently by Thomas (thomas27, ) and Fermi (fermi27, ; fermi28, ). They proposed a model for calculating atomic properties based on $\rho(\mathrm{\mathbf{r}})$ alone. For the kinetic energy, a local density approximation of a uniform, homogeneous electron gas model was suggested. On the other hand, the electron-electron repulsion energy was obtained from the classical Coulomb potential. This leads to the familiar Thomas-Fermi (TF) Euler-Lagrange equation for the density as,

[TABLE]

where $C_{\mathrm{F}}=\frac{3}{10}(3\pi^{2})^{\frac{2}{3}}$, $v_{\mathrm{ext}}(\mathrm{\mathbf{r}})$ signifies the external potential at point $\mathrm{\mathbf{r}}$ consisting of nuclear-electron attraction and $\mu_{\mathrm{TF}}$ is, as usual, the Lagrange multiplier, to be identified as chemical potential. Thus we have the first genuine density functional for the energy of an interacting system. Though TF theory approximately describes charge density, electrostatic potential and total energy (at least for atoms and periodic solids), there are serious deficiencies due to a lack of proper description of the charge density at outer regions of atoms. Further, it is unable to produce the periodic variation of many properties with changing atomic number Z, due to the lack of a shell structure in the atoms. Besides, it does not bind atoms (neutral or ionic) to form molecules and solids. In a way, its utility is very limited to atomic systems.

Within this approximate framework, Dirac (dirac30, ) first made a correction to the TF energy, incorporating exchange effects into it. He derived an exchange term from exchange energy of a homogeneous electron gas of density $\rho(\mathrm{\mathbf{r}})$ by reforming the HF theory solely in terms of density function. This gives the Thomas-Fermi-Dirac (TFD) equation as,

[TABLE]

where $C_{\mathrm{x}}=\frac{3}{4}(\frac{3}{\pi})^{\frac{1}{3}}$. The second term corresponds to the external potential energy containing the nuclear-electron attraction, and the last expression refers to the classical electrostatic Hartree repulsion energy. While the sheer simplicity of replacing the troublesome many-electron SE by a single equation in terms of $\rho(\mathrm{\mathbf{r}})$ alone is manifestly appealing, underlying approximations are crude and grossly inadequate to be of any practical use in quantum chemistry.

As overwhelming complications are met to advance beyond the primitive levels of the TFD model, this line of thought went into oblivion until 1964, when the future hopes were rekindled in a pioneering work by Hohenberg and Kohn (HK) hohenberg64 . It laid a rigorous theoretical foundation of modern-day DFT by asking a plain obvious question as, whether the information contained in $\rho(\mathrm{\mathbf{r}})$ is sufficient for elucidating a many-electron system completely. Following an amazingly simple proof, they answered the question in an affirmative way. In a landmark paper (hohenberg64, ), they first proved that a non-degenerate ground state of an $N$-particle system, moving under the influence of their mutual Coulomb repulsion and an external interacting potential $v_{\mathrm{ext}}(\mathrm{\mathbf{r}})$, is solely characterized by single-particle density $\rho(\mathrm{\mathbf{r}})$, i.e., $v_{\mathrm{ext}}(\mathrm{\mathbf{r}})$, $\psi(\mathrm{\mathbf{r}})$ and hence all other ground-state properties are uniquely determined by $\rho(\mathrm{\mathbf{r}})$ alone (or as unique functionals of $\rho(\mathrm{\mathbf{r}})$). They also showed that, for a given $v_{\mathrm{ext}}(\mathrm{\mathbf{r}})$, the energy functional is minimum corresponding to the true density. Further, the total energy and density of a given system can be determined variationally by minimizing the functional $E[\rho(\mathrm{\mathbf{r}})]$ subject to the normalization condition, $\int\rho(\mathrm{\mathbf{r}})d\mathrm{\mathbf{r}}=N$, as a constant, through the following equation,

[TABLE]

where $N$ is a measure of the total number of electrons present. This is the central equation of DFT providing a deterministic route for $\rho(\mathrm{\mathbf{r}})$. Later, it has been established that the HK theorem is equally valid for degenerate ground states and has been extended to excited states.

Thus the one-to-one mapping between density and energy, the basic preamble of the HK theorem, can be succinctly expressed by conveniently writing the energy functional as follows,

[TABLE]

Here $F[\rho]$ signifies the so-called universal energy density functional, independent of $N$ and $v_{\mathrm{ext}}(\mathrm{\mathbf{r}})$, while $T,V_{\mathrm{ne}}$, $V_{\mathrm{ee}}$ refer to the usual kinetic, electron-nuclear attraction and electron-electron repulsion energy operators in the Hamiltonian respectively. The original HK theorem implies that ground-state properties are functionals of $v$-representable densities only and not of any other trial densities. Later this restriction was lifted; thus widening its domain of applicability from $v$-representable to $N$-representable densities (one that can be obtained from an antisymmetric wave function).

In order to make sure that a given density is indeed the true ground-state density, now second HK theorem gives the crucial variational theorem,

[TABLE]

where $E[\rho_{0}]$ corresponds to the ground-state energy of a Hamiltonian with $v_{\mathrm{ext}}(\mathrm{\mathbf{r}})$ as external potential, and $\rho_{0}$ is the true ground-state density. In other words, the functional attains its minimum value with respect to all allowed densities if and only if our input density corresponds to the true ground-state density. These two theorems, although, validate previous efforts of TFD to utilize $\rho(\mathrm{\mathbf{r}})$ as a basic fundamental variable to describe many-electron systems, emphasis must be given to the fact that it is merely a proof of existence only. It confirms the much sought-after unique mapping between density and energy, in principle; however it remains absolutely quiet on any guidelines for construction of such a functional. Even though $\rho(\mathrm{\mathbf{r}})$ is sufficient, the relation is delicate, intricate, and engenders tremendous complexity. Another major discomfiture lies in the fact that minimization of $E[\rho]$ is, in general, an unusually challenging numerical task. Nature of actual calculations before and after HK theorem, does not change in any noticeable way; for they are as hard as before. No clear-cut clue is available to approximate the unknown functional. Since this does not provide any simplification over MO theory, the ultimate step still remains to be the solution of an excessively clumsy SE equation. It is of very little comfort for chemistry and physics, leading to virtually no actual progress for realistic calculation. Another unpleasant predicament is that the density variation principle holds true only for exact functionals. This implies that, in contrast to standard wave-function-based approaches (such as HF or CI), where wave functions are rigorously variational, within the realm of HK DFT, energy offered by a trial functional has no physical meaning whatsoever.

II.2 The single-equation approach

At this stage, we will briefly mention about the so-called single-equation approach in DFT, which functions exclusively in terms of $\rho(\mathrm{\mathbf{r}})$ alone precluding the consideration of orbitals. Thus this preserves the true spirit of DFT without sacrificing its associated physics. A real crux of problem for development of such methods lies exclusively on the unavailability of accurate (not to speak of exact) kinetic and XC functionals. Nevertheless two major attempts are noteworthy. The first one invokes a link between quantum fluid dynamics and DFT through $\rho(\mathrm{\mathbf{r}})$. The former was proposed long times ago by Madelung, de Broglie, and Bohm based on a hydrodynamic formulation of quantum mechanics. It is more complicated than the Schrödinger picture, in a sense that the former requires solving a set of nonlinear partial differential equations, in contrast to linear equations in the latter. However the charm is that this route has a conceptual appeal, for the electron cloud is treated as a “classical fluid” moving under the influence of classical Coulomb forces and an additional quantum potential. Moreover, amplitude and phase of the wave function (treated explicitly in the fluid formulation as independent variables) are more slowly varying in time than the wave function itself, which leads to computational advantages and offers additional physical insights of quantum dynamics or the quantum trajectory. Thus all the electrons, assumed to be distributed over a 3D space-like continuous fluid, being governed by two basic equations in terms of local variables, $\rho(\mathrm{\mathbf{r}},t)$ and the current density $\mathbf{j}(\mathrm{\mathbf{r}},t)$, can be eventually described by a set of equations, viz., (atomic units used unless otherwise mentioned):

i) an equation of continuity:

[TABLE]

and ii) an Euler-type equation of motion:

[TABLE]

where $\mathbf{j}(\mathrm{\mathbf{r}},t)=\rho\nabla\chi(\mathrm{\mathbf{r}},t)$, with $\chi(\mathrm{\mathbf{r}},t)$ being the velocity potential. $\mathrm{E_{el-el}}$ represents the inter-electronic Coulomb repulsion energy, $v_{\mathrm{ext}}(\mathrm{\mathbf{r}},t)$ is the external potential containing electron-nuclear attraction and any time-dependent (TD) interaction present in the system, while $\mathrm{G}[\rho]$ is a universal functional comprising of kinetic and XC energy contributions. Equations (9) and (10) can be combined into a single equation by defining a complex-valued TD hydrodynamical wave function for the entire time-evolving system as (in polar form):

[TABLE]

and eliminating $\chi(\mathrm{\mathbf{r}},t)$ from them. The result is a generalized nonlinear SE:

[TABLE]

where $\rho(\mathrm{\mathbf{r}},t)=|\Psi(\mathrm{\mathbf{r}},t)|^{2}$, and the effective potential $\mathit{v_{\mathrm{eff}}}([\rho];\mathrm{\mathbf{r}},t)$ contains both classical and quantum potentials. This procedure has been successfully applied to a variety of interesting static and dynamic many-electron situations. Some of these are electronic structure calculations (roy02b, ), ion-atom collisions (deb89, ; dey98a, ), atoms, molecules in presence of strong TD electric and/or magnetic fields (dey98, ; roy02a, ; sadhukhan11, ; vikas11, ; vikas13, ; sadhukhan16, ) including high-harmonic generation and multi-photon ionization, etc. Many more details could be found in the above references and therein.

The other one is the so-called orbital-free DFT (wang00, ), whose realization most decisively depends on the availability of correct, good-quality kinetic energy density functionals. Reportedly, it is possible to compute all components of energy functionals effectively in momentum space within a linear scaling process, which is certainly advantageous in terms of simplicity and computational cost. During the past several decades, an immense amount of work has been attempted towards designing suitable kinetic energy functionals (wang98, ) including a density dependent kernel (wang99, ) and non-local functionals (ho08, ), to mention a few; mostly on periodic system. Efforts were also made to supplement the density-dependent kernel to non-periodic systems, with the aid of finite element modelling and coarse graining (gavini07, ). Later, its scope was enlarged to covalent systems and semiconductors (huang10, ) as well. Up to now, applications are predominantly limited to metallic system.

II.3 The Kohn-Sham method

While the single-equation route of the previous subsection provides an elegant and attractive direction to deal with an interacting system and has witnessed numerous successes, insurmountable difficulties (as delineated in II.A) need to be overcome for its routine practical application. Unsurprisingly, a massive amount of today’s work actually follows an alternative formalism, originally suggested by Kohn-Sham (KS). The prevailing unyielding state of affairs after publication of the HK theorem dramatically changed a year later, in a path-breaking contribution kohn65 , that suggested a manageable way to approach the hitherto unknown universal functional. In order to mitigate the problem, they introduced the clever concept of a fictitious, non-interacting system built from a set of orbitals (KS orbitals) such that the major part of the kinetic energy can be computed exactly. The remaining fairly small portion is absorbed in our non-classical contribution of the electron-electron repulsion, which is also unknown. However, the advantage is that electrons now move in an effective KS single-particle potential, while the mapped auxiliary system yields the same ground-state density as our real interacting system, but this simplifies the actual calculation tremendously.

This is accomplished by setting up a non-interacting reference system (designated with suffix ‘s’) with the Hamiltonian containing an effective local potential $v_{\mathrm{s}}(\mathrm{\mathbf{r}})$ expressed as a sum of one-electron operators,

[TABLE]

Note that there is no electron-electron interaction in the Hamiltonian. It is well known that the exact wave function of a system of non-interacting fermions can be represented as Slater determinants of individual one-electron eigenfunctions, and energies are sums of one-electron eigenvalues. In complete analogy to the HF equation, thus, one can write an eigenvalue equation for KS spin-orbitals as,

[TABLE]

where the one-electron KS operator $\hat{f}^{KS}$ is defined by,

[TABLE]

Hence the crucial step is to map our artificial non-interacting system to the interacting, real one through an effective potential $v_{\mathrm{s}}(\mathrm{\mathbf{r}})$ which plays a pivotal role. This is ensured by choosing $v_{\mathrm{s}}(\mathrm{\mathbf{r}})$ such that the density resulting from a summation of moduli of square of orbitals (\big{\{}\psi_{i}\big{\}}), is exactly the same as the ground-state density of our real system ($\sigma$ signifies spin), i.e.,

[TABLE]

So the leading part of the kinetic energy of the real system can be restored adequately as a sum of individual kinetic energies of the reference system with same density as real system,

[TABLE]

However, since the non-interacting kinetic energy is not equal to the true kinetic energy, they suggested that the residual, kinetic energy ($T_{\mathrm{c}}$), often a small contribution, is submerged to the unknown, non-classical component of electron-electron repulsion, namely,

[TABLE]

Associated terms have following meanings: $J[\rho]$ is the known classical part of $E_{\mathrm{ee}}[\rho]$, whereas $E_{\mathrm{xc}}[\rho]$ contains everything that is unknown, i.e., non-classical electrostatic effects of $E_{\mathrm{ee}}[\rho]$ as well as the difference between the true kinetic energy $T[\rho]$ and $T_{\mathrm{s}}[\rho]$. Then the expression of total energy $E[\rho]$ of our real system can be cast in the following way,

[TABLE]

This allows us to write the celebrated KS orbital equation in canonical form (henceforth atomic units employed unless otherwise stated),

[TABLE]

where the “effective” potential contains following terms,

[TABLE]

In the above equation, $v_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]$ signifies the functional derivative, $\frac{\delta E_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]}{\delta\rho(\mathrm{\mathbf{r}})},$ with respect to charge density. So the KS equation, although exact in principle, is structurally similar to the HF equation, but less complicated than the former. Both equations must be solved iteratively by self-consistent field. It has a profound effect on the quantum community because of its ability to account for XC effects in a rigorous, quantitative and transparent manner with a cost level of HF theory. However, it may be noted that, although being local in nature, the KS potential may be formally innocuous and less complicated than its counterpart in HF approximation, it may have a rather not-so-straightforward and non-local dependence on $\rho(\mathrm{\mathbf{r}})$.

Thus, as explained above, for practical purposes, minimization of our explicit energy functional, although possible, in principle, is hardly ever recommended. A far more viable and feasible route invokes solution of the KS Eq. (20), which, however, interestingly brings back an orbital picture, and strictly speaking, does not function solely in terms of density. But, the good thing is that we have a mechanism to incorporate the intractable many-body effects formally exactly, within a single-particle theory.

II.4 Grid consideration

For realistic solution of the KS equation, it is necessary to deal with mathematically non-trivial integrals that cannot be evaluated analytically and pose certain amount of challenge. Even in finite basis-set expansion methods, there does not exist any explicit analytical formulas for requisite XC integrals. In absence of such procedures, one is left with no option but for numerical calculation. It is well acknowledged that such a discrete procedure for multi-center integrals in 3D space is not straightforward. The task becomes all the more formidable from a consideration of the fact that, in order to achieve a satisfactory level of chemical accuracy, a generous grid size and a longer computation time is often mandatory, as a prohibitively huge number of operations is needed. A vital problem for constructing such integration route arises due to cusps in the density and the singular nature of the Coulombic potential. Naturally, in order to pursue high-accuracy calculations within a reasonable number of quadrature grids, one needs efficient and sophisticated numerical integrators which can capture the forms of density at a satisfactory level. This paves the way for a considerable number of integrators having varied degrees of performance. Amongst them, two distinct, well recognized partitioning schemes have shown significant promise. The Voronoi cellular approach was originally proposed by te Velde and Baerends (te92, ), in which the space is divided into non-overlapping regions of simple geometry. On the other hand, fuzzy cells avenue, or commonly known as ACG, was initially recommended by Becke (becke88, ); later championed by several others groups (gill93, ; treutler95, ; delley90, ; lindh01, ). Though the basic idea was implemented for completely numerical, non-basis set calculations of polyatomic molecules, now it has reached the status of standard for linear combination of atomic orbitals (LCAO)-MO DFT approach, which is convenient while using spherical grids.

The basic strategy is to divide the space into fuzzy (atom-centered overlapping analytically continuous regions) cells. In order to integrate over each cell, a spherical grid (corresponding to radial and angular quadrature grid) centered on each atom is the natural choice. The integral ($I$) of an arbitrary integrand $F(\mathrm{\mathbf{r}})$ over whole region is then converted to discrete numerical summations over individual atomic regions $j$:

[TABLE]

where

[TABLE]

i.e., a multi-center integration therefore translates to a sum of single-center integrations over each atomic nucleus, whereas the $w_{j}$’s correspond to relative weight functions over each atomic center respectively. In order to make the discrete cell fuzzy, certain restrictions are imposed on $w_{j}$, so that,

[TABLE]

and each $w_{j}(\mathrm{\mathbf{r}})$ takes value unity in the neighborhood of its own nucleus, but vanishes near any other nucleus in a well-behaved manner.

Traditional cellular separation of molecular space is best realized by using confocal elliptical coordinates ($\lambda,\mu,\phi$). Taking center $i$ as reference, one considers other centers $j\neq i$ and establishes $\lambda_{ij}$, $\mu_{ij}$, $\phi_{ij}$ on the foci $i,j$. Of particular interest is,

[TABLE]

where $r_{i},r_{j}$ identify distances to nuclei $i,j$ respectively, and $R_{ij}$ the internuclear distance. In terms of an auxiliary function, $s(\mu_{ij}$), normalized weight functions, $w_{i}(\mathrm{\mathbf{r}})$, are defined as:

[TABLE]

The most recommended form of $s(\mu_{ij}$) is a cut-off function as given below,

[TABLE]

where $f_{k}(\mu)$ indicates a $k$-th order iterative polynomial. The forms of polynomial, $f_{k}(\mu)$ and $k$ are chosen such that $w_{i}(\mathrm{\mathbf{r}})$ satisfies the restriction imposed by Eq. (24). There is no hard and fast rule to evaluate the value of $k$; it is chiefly decided through experience ($k=3$ is used in some earlier works (becke88, )). Instead of polynomial, if one chooses $s(\mu_{ij})$ as a step function, the fuzzy scheme transforms into Voronoi scheme. Moreover, a size adjustment procedure has been advocated for hetero-nuclear molecules. The above prescription is most suited for molecules; an alternative approach towards periodic system is nicely presented in (franchini13, ).

The single-center integral, $I_{j}$ in Eq. (22) is expressed in a spherical coordinate system as:

[TABLE]

Now we can distribute the 3D integration into 2D integration of variables $\Omega(\theta,\phi)$ and 1D integration of radial variable $r$ in the following fashion, by applying some appropriate $n$-point quadrature formula having weight factor $\omega$ for both angular and radial parts,

[TABLE]

and

[TABLE]

Angular integration is eventually carried out on numerous sizes of a Lebedev grid (lebedev92, ; lebedev99, ) (also called octahedral grid), which is distributed on the surface of a unit sphere in such a way that the resulting point distribution remains invariant with respect to inversion. This grid is a consequence of various quadrature formulas (mclaren63, ) of remarkably high order, where all the spherical harmonics ($Y_{l,m}$) and its squares are integrated accurately. For optimal performance, it is best that the total number of grid points satisfies a certain condition as:

[TABLE]

where L is the maximum degree of the spherical harmonics ($0\leq l\leq L$). Open-ended quadrature schemes are generally advisable; this is achieved via product type of grids quite easily in polar coordinates ($\Omega$) (murray93, ). Efficiency can be further improved using a Lobatto grid, which uses a different criterion for the total number of points $\left(\approx\frac{L(L+1)^{2}}{2}\right)$, for optimal realization treutler95 . However, the majority of modern KS-DFT programs (within basis-set approach) like Gaussian, NWChem, Q-Chem, GAMESS etc., exploit the Lebedev grid.

For radial integration, a fairly decent number of implementations are conceivable. These are typically carried out with the aid of some suitable quadrature rules. Typically, a transformation is applied to map radial grid onto the quadrature grid varying the number of sampling points on the grid. Generally, a grid-pruning procedure is also administered in conjunction, to curtail the number of sampling points (el04, ; chien06, ). Some commonly used quadrature grids are: Gauss-Chebyshev integration of second kind, Euler-Maclaurin summation (extended trapezoidal) formula, log-squared quadrature, Gauss-Legendre, Gauss-Laguerre quadrature etc., (becke88, ; murray93, ; treutler95, ; mura96, ; gill03, ; kakhiani09, ). The Gauss-Chebyshev integration formula is given as,

[TABLE]

whereas the Euler-Maclaurin scheme is executed through,

[TABLE]

Subsequently, a double exponential formula has been proposed, and later refined by several other groups (mori85, ; mori01, ; muhammad03, ) to generate a radial grid. For large $N$, it becomes beneficial in the way of providing fast convergence rate. This is achieved by adopting an equally meshed trapezoidal rule for numerical integration over an arbitrary interval that transforms into an infinite integral by variable transformation. It has several variants depending on the mapping transformation and was invoked for finite, semi-finite, and infinite integrals. Another promising approach involves construction of an adaptive grid, pioneered by (krack98, ; perez94, ; termath96, ) which is automatically generated for a given accuracy corresponding to a chosen basis function. Some of the prominent DFT programs like NRLMOL, deMon2k routinely use such grids.

It is well-known that the current enviable status of DFT that it enjoys largely depends on basis-set calculations. While such studies are heavily dominated by ACG, real-space grid has been invoked for fully numerical, basis-set free DFT methods. Apart from ACG, some scattered works exist for other grids in literature, e.g., an adaptive Cartesian grid with a hierarchical cubature method (challacombe00, ), a transformed sparse-tensor product grid rodriguez08 , a Fourier Transform Coulomb (fustimolnar02, ; fustimolnar03, ) method interpolating density from ACG to a more regular grid that helps to enhance its efficiency. In the latter, a difficulty arises due to computation of XC matrix elements, which is circumvented successfully with help of a multi-resolution method (brown06, ). Here, in this work, we report a simple fruitful DFT implementation roy08 ; roy08a ; roy09 ; roy10 ; roy11 within the LCAO-MO framework which solely uses CCG. So far, our results have been reported in a uniform grid. All relevant quantities such as basis functions, electron densities, MOs as well as various two-electron potentials (Hartree and XC) are directly set up on a real 3D Cartesian grid simulating a cubic box as,

[TABLE]

where $h_{r}$,$N_{r}$ signify grid spacing and total number of grid points respectively ($r_{0}=-\frac{N_{r}h_{r}}{2}$).

A major concern in the grid-based approach constitutes an accurate estimation of the classical electrostatic repulsion potential. For finite systems, the simplest and crudest way to calculate $v_{\mathrm{h}}(\mathrm{\mathbf{r}})$ is through direct numerical integration on the grid. For smaller systems, this is a feasible option; in all other cases, it is generally tedious and cumbersome. However, the most rewarding and widespread approach is through a solution of the corresponding Poisson equation,

[TABLE]

The usual way to solve this is by conjugate gradients (saad03, ) or through multi-grid solvers (brandt77, ). As an alternative, the current work exploits a conventional Fourier convolution method originally suggested by (martyna99, ; minary02, ; skylaris02, ) and adapted in the context of molecular modelling (hine11, ; chang12, ). The basic principle can be formulated as:

[TABLE]

where $v_{\mathrm{h}}^{\mathrm{c}}(\mathbf{k})$ and $\rho(\mathbf{k})$ stand for Fourier integrals of the Coulomb interaction kernel and density respectively, in the cubic box. The quantity $\rho(\mathbf{k})$ can be easily computed by using a discrete Fourier transformation of its real space value. Thus our primary concern here lies in calculation of the Coulomb interaction kernel which has a singularity in real space. For this, we utilize an Ewald summation-type approach (chang12, ), expanding the Hartree kernel into long-range and short-range components,

[TABLE]

where erf(x) and erfc(x) signify the error function and its complementary function respectively. The Fourier transform of the short-range part can be treated analytically, whereas the long-range portion needs to be computed directly from FFT of real-space values. A convergence parameter $\alpha$ is used to adjust the range of $v^{\mathrm{c}}_{\mathrm{h}_{\mathrm{short}}}(\mathrm{\mathbf{r}})$, such that the error is minimized. Following the conjecture of martyna99 , here we employ $\alpha\times L=7$ ($L$ denotes length of cubic box), which produces quite accurate results. For a non-uniform grid, $L$ is chosen as the smallest side of the simulating box. Some other routes which calculate the long-range part efficiently are fast multipole (beatson97, ), multi-level summation (skeel02, ), fast Fourier-Poisson (york94, ) method, etc.

II.5 Basis-set free methods

This sub-section gives an overview of fully numerical, real-space methods, used for solving the pertinent KS equation without needing any explicit basis set. Typically, the wave function is directly sampled on a discrete grid through a variety of representations; each having its own pros and cons. Most of these attempts also invoke pseudopotentials to freeze the inner core. One central problem is the formation of an appropriate finer mesh, which normally calls for an enormous number of grid points, making calculations of the electronic structure (especially for large systems) very expensive in terms of both time and computer memory. Apparently, there are ways to minimize discretization error systematically by controlling size, shape and spacing of the 3D mesh. It is worth mentioning two notable works, viz., a higher-order real-space pseudopotential method in uniform CCG chelikowsky94a ; chelikowsky94b , as well as real-time propagation of KS orbitals mundt07 ; marques12 . Often, they are easily amenable to linear scaling schemes. Recently, extensions have been made to Graphical Processing Unit interface (andrade13, ) and domain-decomposition techniques chelikowsky09 for parallel computation. Some excellent reviews are available on the subject (see, for example, beck00 ; natan08 ; bowler12 ; bernholc08 ).

Broadly speaking, real-space techniques can be categorized into three major classes, namely, finite difference (FD), finite element (FE) and wavelets (mohr14, ; genovese08, ). In all cases, the target discretized differential equation produces structured and highly banded matrices, which can be solved readily using efficient multi-scale techniques. The potential operator is diagonal in coordinate space and the Laplacian operator is nearly non-local; that makes them particularly suited for parallel computing. In a uniform orthogonal 3D CCG of spacing $h$, e.g., the Laplacian operator, within an FD approximation, can be expressed approximately in 1D ($x$ direction) as a part of 3D:

[TABLE]

The last two terms account for truncation errors; the first contribution being second order in $h$ with a pre-factor involving the fourth derivative of $\phi(x)$. Here grid spacing, $h$, is one of the key parameters that controls convergence. This approximation is valid in the limit of $h\rightarrow 0$, but for practical purposes, it is unfeasible to go below a certain limit, which leads to the evolution of higher-order FD schemes chelikowsky94a . Most commonly, regular CCG is engaged and FD methods result from a Taylor series expansion of our desired function around the grid points. However, it is non-variational, because the error can be of either sign depending on the derivatives and value of $\phi(x)$. It is often quite difficult to achieve convergence due to this lack of variationality. While its chief advantage lies in a simplicity of representation and easy implementation, sometimes it is rather troublesome to construct flexible meshes which can reproduce physical geometry satisfactorily. There are reports of a Mehrstellen discretization technique as an alternative route to high-order FD methods (briggs96, ; fattebert00, ). Several workers (alemany07, ; Iwata10, ; fujimoto10, ; heiskanen01, ), have suggested use of FD techniques with in amalgamation of projector augmented wave methods (mortensen05, ). Applications of the FD method is quite impressive; mention may be made about some of the notable ones for clusters and other finite systems modine95 ; lee00 ; chelikowsky94a ; chelikowsky94b .

On the other hand, FE-based methods are well known from the engineering field for a long time, and have been implemented on applications of electronic structure theory (pask05, ). Customarily, one uses some suitable non-orthogonal basis functions, such as piecewise polynomials which are non-zero over a local region of space, and also divide the unit cell into elements. For solution of Poisson equation, one can expand the potential, $v(x)$, in a basis:

[TABLE]

where $u_{i}$’s are expansion coefficients corresponding to actual function values on mesh points and $\chi_{i}$ are basis functions. A similar kind of treatment can be performed on the charge density, and the resulting matrix equation will be similar to a one-dimensional Poisson equation. At least for uniform meshes, there is a close correspondence between FD and FE real-space representation. Unlike the FD method however, it has a variational foundation, and provides greater flexibility regarding the unit cell construction. A multitude of efficient basis sets have been suggested for 3D electronic structure calculations for diverse problems (tsuchida95, ), like cubic-polynomial basis white89 , tetrahedral discretization with orders $p=1-5$ (ackermann94, ), cubic functions (pask99, ), Lobatto-Gauss basis set with orders $p=5-7$ (yu94, ), B-spline basis (goringe97, ) etc.

Standard iterative processes suffer from so-called critical slow down beck00 ; a common complication that occurs when more grid points are used to obtain an enhanced accuracy on a fixed domain and thus makes it less efficient on finer meshes. In order to facilitate better convergence of self-consistent iteration procedure, a few advancements have been made. Most prominent of them is the family of multi-grid iteration techniques briggs95 ; briggs96 ; beck00 , which ideally utilize information from multiple lengths scales to alleviate this problem. Another promising route involves engaging an adaptive grid in conjunction with multi-grid (gygi95, ; modine95, ). This has found applications within FD-type representations, but maximum effort was devoted to develop efficient solvers for FE representations (braess90, ; brenner94, ).

Another useful source of real-space techniques is based on the wavelet basis method having a more complicated matrix structure than either the FD or FE representation (arias99, ). Such basis sets are semi-cardinal, local, and not only that, they are conceived with multi-resolution properties. One crucial advantage is that calculations can be performed with different kinds of boundary conditions. Daubechies wavelet (genovese08, ) is one of the most favored and frequently used, which is local in both real and Fourier space. These are easily adaptable to parallel programming and are employed in connection with linear scaling methods mohr14 for electronic structure calculations. Several DFT codes like ABINIT, ONETEP, CONQUEST, CP2K, and SIESTA etc., successfully implement wavelet schemes.

II.6 Basis set and LCAO-MO ansatz

As emphasized earlier, the basis-set approach remains by far the most convenient and pragmatic line of approach towards the solution of KS DFT. It dominates in every field of science whether being chemical, physical or biological applications employed to atoms, molecules or solid etc. The inspiration mainly comes from the success of basis-set related methodologies in traditional wave function theory, such as HF and post-HF. This is a coupled integro-differential equation where kinetic and potential energies are defined by a differential and integral operator. It must be solved numerically iteratively leading to a self-consistent set of orbitals $\{\psi_{i}(\mathrm{\mathbf{r}})\}$. Essentially, the unknown KS MOs are expanded in terms of $K$ suitably chosen, known basis functions $\{\chi_{\mu}(\mathrm{\mathbf{r}});\mu=1,2,3,....,K\}$, conventionally called atomic orbitals, in a manner analogous to that in the Roothaan-HF method, such as,

[TABLE]

The electron density then takes the following expression in this basis,

[TABLE]

In principle, one requires a complete basis set ($K=\infty$) in order to get an exact expansion of MOs, but in reality it is not feasible. Therefore, suitable truncation is required for realistic computational purposes; it suffices to work with a finite basis set usually.

It is of paramount importance to choose appropriate functions which reproduce KS MOs as precisely as possible. The numerical accuracy of KS DFT is quite sensitive to the choice and design of the basis set for a particular problem, as an incomplete basis set inducts certain constraints on the relaxation of the density through KS orbitals. Considerable developments have taken place towards their construction as well as effects on assorted physico-chemical properties. There is no restriction, as such, for the effective implementation of a specific basis set, but it is required to fulfill certain conditions which may serve as some guidelines: (a) it should have proper behavior according to the physics of the problem (b) it should be friendly in terms of ease and efficiency of computation of requisite integrals (c) number of functions in basis-set expansion should be scaled down as much as possible without sacrificing any accuracy. A sizeable number of elegant, flexible, versatile basis sets has been proposed over the past several decades, from different perspectives such as: exponential (slater30, ), Gaussian (boys50, ), polynomial, cube functions, wavelets (cho93, ), numerical functions (delley90, ), plane waves (ihm79, ) and many others. However, the problem still now remains open, as there is no universal basis set applicable for all methods or chemical systems. The interested reader may consult many nice reviews and books available on the topic (szabo96, ; helgaker00, ; cramer04, ; jensen07, ).

For a periodic system (e.g., solid), the preferred choice is a plane wave basis whereas for non-periodic systems (e.g., molecules, clusters) atom-centred localized basis functions, such as Slater Type Orbitals (STO) or Gaussian Type Orbitals (GTO) stand out to be two most attractive options. STOs assume the following functional form:

[TABLE]

where $N_{c}$ is the normalization constant and $Y_{l,m}$ defines spherical harmonics. It is a well acknowledged fact that multi-center (3c or 4c) two-electron integrals can not be evaluated analytically with STOs, making them somehow less favorable for molecular system. Thus their application principally remains confined to atomic or diatomic systems or in semi-empirical approximations, where two- or three-center integrals are generally neglected.

On the other hand, GTOs are expressed in Cartesian coordinates as:

[TABLE]

where $n_{x}$, $n_{y}$ and $n_{z}$ define the angular momentum such that $\lambda=n_{x}+n_{y}+n_{z}=0,1,2,\cdots,$ ($s,p,d,\cdots$ functions), while the orbital exponent $\zeta$ governs its compactness (large $\zeta$) or diffuseness (small $\zeta$). In spite of the fact that GTOs face problems to describe the proper behavior near and far from the nucleus than STOs, they are almost universally accepted for quantum chemistry calculations. To a great measure, this is due to their ability to provide closed, analytic formulas for multi-center integrals through Gaussian product theorem. In principle, both STOs and GTOs qualify to form a complete basis, but to reach a certain accuracy, ordinarily more GTOs are needed (roughly at least three times) than STOs. In practical purpose, instead of using individual GTOs as basis function, it is customary to use fixed linear combination of GTOs, called contracted GTOs, defined as:

[TABLE]

Here $d_{p\mu}$ are called contraction coefficients, $L$ signifies length of contraction, whereas “p” stands for individual primitive Gaussians from which contracted Gaussian functions are formed. The normal procedure to find out primitives is to optimize exponents and contraction coefficients so as to obtain variational lowest-energy states through an atomic calculation (within HF, CI or so). However there are instances in literature where they are explicitly optimized through KS procedure using LDA XC functionals. Historically, basis sets were designed for applications in wave-function based methodology, and fortunately it turns out that, for most common molecular property calculations (like equilibrium geometry, energy), results are pleasantly insusceptible towards the manner through which these exponents and coefficients have been determined. That offers a welcome possibility to take advantage of such basis sets, initially constructed for traditional methods, to be directly usable in DFT calculations with confidence.

Now inserting Eq. (40) in Eq. (20), multiplying the left side of the resulting equation with $\chi_{\mu}^{\star}(\mathrm{\mathbf{r}})$, then integrating over the whole space, followed by some algebraic manipulation generates the following KS matrix equation, in parallel to the HF case,

[TABLE]

where $\mathbf{F}$ and $\mathbf{S}$ imply the $K\times K$ real, symmetric total KS and overlap matrices respectively. The eigenvector matrix $\mathbf{C}$ contains the basis-set expansion coefficients $\mathbf{C}_{\mu i}$ and the diagonal matrix $\mathbf{\epsilon}$ holds the orbital energies $\epsilon_{i}$. It could be readily solved by standard numerical techniques of linear algebra. The individual elements of KS matrix can be written as:

[TABLE]

where $H_{\mu\nu}^{\mathrm{core}}$ represents the core bare-nucleus Hamiltonian matrix element consisting of the kinetic energy and nuclear-electron attraction, thus accounting for one-electron energies. The one-electron matrix elements can be evaluated analytically with the help of well-established recursion relations (obara86, ) for Gaussian bases. Second term $v_{\mathrm{hxc}}(\mathrm{\mathbf{r}})$ contains all two-electron interactions involving classical Coulomb repulsion and the non-classical XC potential. $J_{\mu\nu}$ denotes the matrix element of the classical Hartree potential defined in Eq. (21), while the remaining term, $V_{\mu\nu}^{\mathrm{XC}}$ supplies the XC contribution into two-body matrix element, whose construction remains one of the most vital steps in the whole KS DFT process. In absence of any analytical method, two-body matrix elements can be either calculated numerically or fitted by an auxiliary set of Gaussian functions, as suggested by (sambe75, ; dunlap79a, ; dunlap79b, ). In our work, we employed direct numerical integration on the CCG to obtain these matrix elements:

[TABLE]

For gradient-corrected functionals, the non-local XC contribution of the KS matrix is implemented by a finite-orbital basis expansion, without requiring to evaluate the density Hessians. Thus, in such cases, the XC contribution is written in a convenient working form, as suggested in pople92 ,

[TABLE]

where $\gamma_{\alpha\alpha}=|\nabla\rho_{\alpha}|^{2}$, $\gamma_{\alpha\beta}=\nabla\rho_{\alpha}\cdot\nabla\rho_{\beta}$, $\gamma_{\beta\beta}=|\nabla\rho_{\beta}|^{2}$ and $f$ is a function only of local quantities $\rho_{\alpha}$, $\rho_{\beta}$ and their gradients.

II.7 Pseudopotential

In today’s repertoire of computational chemistry/physics electronic structure theory, a pseudopotential approximation is a very useful and powerful concept which is extensively used for atoms (especially with heavy nuclei), molecules (containing one/more transition metal/lanthanide/actinide), solid state, etc. (huang57, ; cohen75, ; leggett01, ). The idea was originally introduced a long time ago by Hellmann, Fermi and Gombas independently; thereafter applied and popularized by several other groups (cohen84, ; schwerdtfeger11, ). It simplifies by dividing the space of many electrons into core and valence categories. To a very good approximation, because of their strong binding with the nucleus, inner-shell electrons form a passive “inert” core (retaining an atomic-like configuration) which supposedly plays a less significant role in the understanding of chemical binding. Consequently, a surprising majority of chemical properties can be followed very satisfactorily by considering only valence electrons. In effective terms, this implies to replace the strong Coulomb potential between nucleus and tightly bound core electrons (as a part of the exact Hamiltonian) by a non-local, smoother potential acting on the valence electrons. That means, the spectrum of the resulting pseudo-Hamiltonian is as close to the exact one and this is also true for the pseudo eigenfunction. It diminishes computational cost in larger molecular systems dramatically, simply by cutting down the basis set, and thus avoiding any integral arising from core electrons. A lot of developments have been recorded regarding the suitability and effectiveness of this approach, compared to the more complete and rigorous, all-electron calculations (maron98, ; schwerdtfeger95, ; schwerdtfeger00, ). Moreover, relativistic effects can be brought into a many-electron picture efficiently through pseudopotential schemes (batyrev01, ; schwerdtfeger11, ).

Thus a molecular pseudo-Hamiltonian $H_{v}$ can be factored as:

[TABLE]

where $n_{v}$ denotes number of valence electrons, $N_{c}$ the number of core (nuclei) ones and $Q_{a}$ signifies residual charge of core $a$. Indices $(i,j)$ and $(a,b)$ run over all valence electrons and all cores respectively; while $v^{\mathrm{pp}}$, the pseudopotential operator takes account of the fact that valence space is repulsive in short range and attractive in long range. Therefore, the basic challenge is to find out a decent and computer-proficient approximation to $v^{\mathrm{pp}}$ comparable with all-electron calculations. It is to be noted that the distinction between core and valence shells is essential for the success of the pseudopotential approximation and hence, it is beneficial to work within an orbital-based theory. Under this situation, our desired KS equation could be rewritten as:

[TABLE]

where $v_{\mathrm{ion}}^{\mathrm{pp}}$ designates ionic pseudopotential for the system,

[TABLE]

In above equation, $v_{\mathrm{ion},a}^{\mathrm{pp}}$ represents the ion-core pseudopotential associated with atom $a$, situated at $\mathrm{\mathbf{R}}_{a}$; whereas $v_{\mathrm{h}}[\rho(\mathrm{\mathbf{r}})]$ denotes the classical electrostatic potential among the valence electrons and $v_{\mathrm{xc}}[\rho(\mathrm{\mathbf{r}})]$, as usual defines the non-classical part of the many-electron Hamiltonian.

At present, two possible approaches are predominant for pseudopotential approximations in molecular applications, namely (i) model core potential (huzinaga71, ) and (ii) semi-local approximation. In the former, $v^{\mathrm{pp}}$ is formulated as,

[TABLE]

where $A_{k}$, $n_{k}$ and $\alpha_{k}$ are all adjustable parameters and $B_{c}$ is a chosen parameter such that $B_{c}$ = $-2\epsilon_{c}$, the core energies ($k$ runs over Gaussian expansion). It is possible to get HF orbital energies and corresponding radial functions (zeng09, ) by adjusting above equation. An advantage of this representation is that the inner-shell structure of the pseudo-valence orbitals is conserved and closely resembles the all-electron valence orbitals; as a result, relativistic and spin-orbit effects can be efficiently included into it (fedorov02, ). Incorporation of non-local core-valence exchange of the all-electron Fock operator into this approximation leads to ab initio model potential approximation (andzelm85, ), which has been profitably implemented into molecular program packages such as Molcas, GAMESS-US etc.

On the other hand, within a semi-local approximation (schwarz71, ; kahn72, ; redondo77, ; preuss81, ; schwerdtfeger82, ; kahn84, ), $v^{\mathrm{pp}}$ is expressed as,

[TABLE]

where $A_{k},B_{lk},n_{k},n_{lk},\alpha_{lk}$ and $\beta_{lk}$ are all adjustable parameters. The last term in above equation contains a projection operator $P_{l}$ which projects onto the Hilbert sub-space with angular momentum $l$. Unlike Eq. (52), it does not contain any core function and hence, uses orbitals that are smooth and node-less in short-range of the radial function. This approximation has found profound applications, especially in relativistic and spin-orbit effects (mcmurchie81, ; fuentealba83, ; lee09, ) and, has been implemented in several molecular program packages like Gaussian90, Molpro, Turbomole and also in solid-state program CRYSTAL. Depending on the manner in which adjustable parameters in Eq. (53) are ascertained, there are two leading categories in this case, namely, (i) energy-consistent (thierfelder10, ) and (ii) shape-consistent (hay85a, ). In the former scenario, these parameters are settled down through a fitting procedure (often least square) figgen09 to a large number of pre-calculated transitions in valence space. In latter plan, this is done in such way that the valence orbitals of different symmetries resemble the shape of the corresponding all-electron orbitals very closely after a certain cut-off radius and also match the orbital energies with high accuracy. This scheme was formulated and implemented by a number of research groups (kahn72, ; redondo77, ; kahn84, ; hay85a, ; hay85b, ; wadt85, ). Additionally it is usually faster and more efficient than the former, though it requires a large number of Gaussian fitting functions for inversion of the Fock equation involved. Moreover, the first procedure requires a vast number of valence spectra to fulfill the shape-consistent requirement. On the other hand, there is a family of non-conserving pseudopotentials (hamann79, ) (slightly different from shape-consistent), which has been found to be most suitable for plane-wave codes within DFT. In our current development, we adopt the latter form as proposed by (hay85a, ; wadt85, ).

III Results and discussion

Before we proceed for a discussion of results, mention may be made about some of the technical aspects. Our desired KS equation is solved following the usual self-consistent procedure by imposing three convergence criteria namely, (i) electronic energy differences between two successive iterations lies below $10^{-6}$ a.u. (ii) maximum absolute deviation in the potential is less than $10^{-5}$ a.u. and (iii) standard deviation in the density matrix remains below $10^{-5}$ a.u. The generalized matrix-eigen value problem is solved using standard LAPACK routine (anderson99, ) efficiently and accurately. Computation of FFT is done by means of the standard FFTW3 package fftw05 .

We begin by scrutinizing the stability and convergence for a representative homo-nuclear diatomic molecule Cl2, at an internuclear distance of 4.2 a.u. For this, the non-relativistic ground-state total energy as well as various energy components are reported in Table I for eight chosen grid Sets A-H, varying in $N_{r}$ and $h_{r}$ ($r\in\{x,y,z\}$). The LDA XC potential, corresponding to homogeneous electron-gas correlation of Vosko-Wilk-Nusair vosko80 , is used in this table, which is quoted from roy08 . All results in this and the following tables compare these with corresponding reference values obtained from the standard GAMESS suite of quantum chemistry program (schmidt93, ) incorporating the same XC functional, basis set and effective core potential. All pseudopotential calculations in this work, utilize the Hay-Wadt-type hay85 valence basis set in which the orbital is split into inner and outer components (derived by two and one primitive Gaussian functions respectively). Reference values are quoted for both “grid” and “grid-free” calculations; former has a default “army” grade grid with Euler-MacLaurin quadratures for radial integration and Gauss-Legendre quadrature for angular integrations. The latter (glaesemann98, ) follows an auxiliary basis-set expansion procedure for molecular integrals instead of a quadrature grid; it has the ability to bypass any error associated with a finite grid, though it suffers from an apparently inherent weakness of incompleteness. Following expectation values of energy components are reported for comparison, viz., kinetic energy $\langle T\rangle$, total nucleus-electron potential energy $\langle V_{\mathrm{t}}^{\mathrm{ne}}\rangle$, total two-electron potential energy $\langle V_{\mathrm{t}}^{\mathrm{ee}}\rangle$, total potential energy $\langle V\rangle$($=\langle V_{\mathrm{t}}^{\mathrm{ne}}\rangle+\langle V_{\mathrm{t}}^{\mathrm{ee}}\rangle+\langle E_{\mathrm{nu}}\rangle$) where $\langle E_{\mathrm{nu}}\rangle$ is nuclear repulsion energy, plus total integrated electron density $N$. Individual components of $\langle E_{\mathrm{hxc}}\rangle$ are not produced as they can not be directly compared (unavailable from output). Clearly, Set A energies deviate maximum from reference values, presumably because our box is too small to account for all important interactions; this is also reflected in a poor $N$ value. As grid spacing is increased from 0.3 of Set A to 0.4 in Set B keeping $N_{r}$ fixed at 32, with a corresponding enlargement of box size, all quantities improve considerably. Sets C, F also furnish similar kind of accuracy as the box length is same; however both are inadequate to reach the correct value. It is interesting to note that Sets B, C, F all correspond to the same box size; but the quality of results is slightly better in case of Set F, as it offers a finer grid structure ($h_{r}=0.1$). A similar conclusion can be drawn for Sets E, G too. Results from Sets D, E are quite comparable to reference values. In order to prove stability, some further calculations are performed by elongating the box in Sets F-H in columns 7-9. This discussion evidently establishes that either of four Sets D, E, G, H offer the desired accuracy; while D, E are capable of providing sufficient accuracy for all practical purposes. A similar analysis was done for HCl roy08 (at $R=2.4$ a.u.) in various grids, giving analogous pattern. Set B seems to be quite reasonable for reproducing literature results up to third decimal place, although it was inadequate for Cl2. All the sets produce results of similar accuracy; Set B amongst them does it slightly poorly in terms of component energies and $N$. As one passes from B-C-D, a gradual improvement in results is recorded, as expected. Sets C and F indicate agreement with each other as in Cl2. For Sets D and E, $N$ remains practically unchanged and energy increases by $\approx$0.0002 a.u. Additional calculations with $N_{r}=128,h_{r}=0.3$ suggests complete agreement with Set G. As in Cl2, the best three sets for HCl are again D, E and G; with the first two being fairly accurate for most calculations.

Up to this point, we focused on results with LDA functionals only. As mentioned in the Introduction, such calculations are fraught with some severe limitations. Hence in order to make practical applications, more refined and accurate functionals are essential. In this context, let us recall that $E_{\mathrm{xc}}$ consists of two components, each arising as a difference: (a) between classical and quantum mechanical inter-electronic repulsion (b) between kinetic energy of the hypothetical, non-interacting and the real, interacting system concerned. Usually, the second portion is not explicitly incorporated in the current family of functionals. In this sub-section, we present a cross-section of CCG results within the familiar non-local Becke exchange becke88a and LYP lee88 correlation, who converted the Colle-Salvetti formula for the correlation energy from HF second order density matrix into a density functional form. This will extend the range and scope of applicability of the present scheme in the context of more realistic non-local functionals which will be crucial for future chemical applications. The latter is implemented using an alternate equivalent form containing only the first derivative through a procedure due to miehlich89 . Table II presents a comparison of our BLYP energy components and $N$ at some selective grid sets at same $R$ as in Table I. The same quantities as in the LDA case are reported. Several test calculations were performed in several grids to ensure convergence of results; which offered conclusions quite similar to those found in LDA. From this, two grid sets are given for each of them. Once again, our CCG results exhibit excellent agreement with literature values, as observed in LDA scenario. In Set B, results are slightly better than Set A for obvious reasons of finer grid and an extended length of box. In Cl2, this effect is more pronounced than in case of HCl. Absolute deviations in Set B energies are 0.00002 and 0.00000 a.u. respectively for Cl2, HCl. Grid-free and grid DFT results are substantially different for both molecules. More details could be found in roy08a .

In order to justify the dependability and usefulness of CCG, Table III now samples calculated negative eigenvalues for Cl2 and HCl using BLYP XC functional at same geometries of previous tables. These are provided in two selected grid sets. Comparison with reference values again indicates excellent agreement for both molecules for all orbital energies except the lone case of $3\sigma_{g}$ for Cl2 (Set A) and $3\sigma$ for HCl (Set B); even then the absolute deviation remains well within 0.0001 a.u. More details could be found in roy08 .

Now, Table IV offers a comparative chart of $\left\langle T\right\rangle$, $\left\langle V\right\rangle$, total (E) energies and $N$ for few chosen atoms and molecules using both LDA and XC functionals. LDA results employ Set E grid ($N_{r}=64,h_{r}=0.4$), while for BLYP calculations, except P and H2S (using Set E), grid Set I ($N_{r}=128,h_{r}=0.3$) was used. Overall, the present calculated results agree with literature values more or less up to two decimal places and $N$ completely matches with the reference values with few exceptions. Results for more atoms and molecules can be found in (roy08, ; roy08a, ).

For further examination, Fig. (1) displays the changes in above LDA energies in Cl2 (relative to $-$29 a.u.) for four selected Sets D, E, G, I and HCl (relative to $-$15 a.u.) for Sets B, C, D respectively. They cover a span of 3.5-5.0 and 1.6-3.1 a.u., for Cl2 and HCl. Respective reference plots are also given for easy comparison. It is quite gratifying that all four sets replicate the shape of the potential energy curve quite well for the entire range of $R$, in Cl2. A thorough analysis of results (not given here and the reader may be referred to roy08 ) reveals that, until $R=4$ a.u., Set D energies are quite close to the reference (higher by only 0.0001 a.u.), while implying a gradual tendency to deviate thereafter. However, the maximum deviation still remains quite small overall (only 0.0007 a.u.) occurring for $R=5$ a.u. Sets G, I computed energies either completely match with reference or provide a maximum difference of 0.0001 a.u. Further in all sets, our CCG energies have always been above reference except in two occasions ($R=$4 and 4.3 for Set G). Now we turn to the bottom panel in the left side for HCl. Once again, we observe excellent qualitative agreement for whole range of $R$, as evident from their nearly identical shapes. Careful examination admits that maximum absolute deviation for three sets B, C, D remains rather small, viz., 0.0014, 0.0007, 0.0002 for $R=$3, 2.3, 2.3 respectively. It is possible to fine-tune these results further by enforcing more rigorous convergence criteria. For more details, see roy08 . Similar plots for BLYP functionals are given in right side for Cl2, HCl in top, bottom panels, for two grids A, B. In both cases, again respective grid-DFT plots are provided for comparison, which practically coincide. A careful analysis (roy08a, ) shows that maximum disagreements are 0.0002 a.u. (only for 2 instances) and 0.0001 a.u. for Sets A, B respectively in case of Cl2. However, for HCl, the same for both sets are 0.0001 a.u. (only for 2 instances). It is a matter of fact that the two sets confirm to each other more for HCl than Cl2.

Now results are given for four selected species in Table V, for two other functionals which are used in applications of many-electron systems with considerable success, namely, PBE perdew96 and Filatov-Thiel (FT97) filatov97a ; filatov97b . These non-local functionals are implemented with the help of Density Functional Repository program repository . For sake of completeness, respective LDA and BLYP values are also provided. Experimental geometries are taken from NIST computational chemistry database johnson06 . Set E grid was used for this table. Since correctness of our results has already been established in previous tables, here we omit the reference values; all entries are very close to them. Total energies decrease as one passes from LDA to non-local functionals. Within the latter category, there is some agreement amongst the three, with FT97 offering the lowest. $N$ maintains same accuracy for all functionals.

Finally Table VI, provides a comparison of calculated $-\epsilon_{\mathrm{HOMO}}$ and atomization energies with experimental results (afeefy05, ) for 4 selected molecules at their experimental geometries (johnson06, ), with LDA and BLYP XC functionals. For both quantities, reference theoretical results (schmidt93, ) are practically identical to those found from CCG, and so omitted here to avoid crowding. It is observed that calculated ionization energies generally signal large deviations (typically off by 30-50%) from experimental results for both functionals. Estimated ionization energies from PBE and FT97 are also reported, which also deliver similar performances. This leads to a general trend, which some of the frequently used functionals suffer from, i.e., even though energies are reproduced properly, ionization energies are described rather poorly. Now, it is a crucial requirement for real-time electron dynamics to have a proper description of orbital energies, especially the higher-lying ones. For this, we present the computed HOMO energies with modified Leeuwen-Baerends (LB) exchange potential (leeuwen94, ; schipper00, ) along with the other four functionals. It is known to give a more satisfactory representation of the long-range nature of the XC potential, as reflected from a better matching with experimental results. This is formulated as:

[TABLE]

where $\alpha=1.19$, $\beta=0.01$ are two empirical parameters, $\sigma$ signifies up, down spins and last term containing gradient correction is reminiscent of the exchange functional of (becke88a, ), while $x_{\sigma}(\mathrm{\mathbf{r}})=|\nabla\rho_{\sigma}(\mathrm{\mathbf{r}})|[\rho_{\sigma}(\mathrm{\mathbf{r}})]^{-4/3}$ is a dimensionless quantity. This ensures the required long-range nature, i.e., $v_{\mathrm{xc\sigma}}^{\mathrm{LB\alpha}}(\mathrm{\mathbf{r}})\rightarrow-1/r,r\rightarrow\infty.$ A combination of LB exchange with the VWN correlation functional produces the HOMO energies of column 6, which abundantly establishes its superiority over all the other four functionals considered here. These are substantially lower and more accurate than four functionals. In case of atomization energy, the present results reveal a considerable deviation from experimental values for both LDA and BLYP functionals. Note that latter contains zero-point vibrational corrections and relativistic effects. Interestingly, LDA atomization energies are apparently better than BLYP counterparts in some cases, but that should not be misconstrued to identify the former as a better choice. There could be some cancellation of errors. Moreover there are other important factors, such as using more complete accurate basis set or more appropriate pseudopotential (the current work uses one designed for HF calculations; those developed for DFT would be more suitable), which could circumvent some of these discrepancies. However, such deviations are not uncommon and found, even in all-electron calculations with more sophisticated, extended basis sets. More discussion on these issues could be found in roy08a . Nevertheless, this is aside the primary objective of this chapter and does not directly impact the main conclusions drawn herein.

So far, we were mainly concerned with valence-shell electrons using some sort of effective core potentials to incorporate the effects of frozen core electrons. Undoubtedly, while the pseudopotential approximation is indispensable for larger systems (especially those containing heavy atoms), “full” calculations provide a more detailed and accurate result at the expense of a heavy computational cost. Nowadays, it is all about a balance between cost and accuracy. If the cost accuracy ratio permits, unquestionably, the latter remains the preferred choice for most chemical and physical problems. Thus, as a further extension of our method, we present some sample results on “full” calculations in Table VII. The electronic energy including the relevant components and $N$, are produced for two representative atoms, namely, N and O in their open-shell ground-state configuration. We employed an STO-3G basis set, LDA XC potential and grid Set I ($N_{r}=128,h_{r}=0.3$). In all cases computed quantities, once again, establish excellent agreement with reference values. This validates the effectiveness and level of performance of our method in the context of all-electron calculations. It also gives a comparison of orbital energies for both $\alpha$ and $\beta$ spins for the same two atoms, which again proves the usefulness. It is a matter of fact that we use very crude basis functions and LDA functional for the purpose of illustration; involvement of better basis sets and non-local XC functionals should enhance the applicability and accuracy of full calculations.

Before passing, a few comments are in order. The computational scaling goes as follows: (a) construction of the localized basis set scales as $N_{g}$, where $N_{g}$ is total number of grid points (b) formation of the basis set in grid scales as $N_{b}\times N_{g}$, with $N_{b}$ giving the number of basis functions (c) generation of the electron density in grid scales as $N_{b}^{2}\times N_{g}$ (d) one- and two-body matrix elements of the KS matrix scale as $N_{b}^{2}$ and $N_{b}^{2}\times N_{g}$ respectively. One obvious way to achieve better scaling is by reducing $N_{b}$ and $N_{g}$. Now, up to this point, all our calculations were performed in a uniform Cartesian grid; further there was a restriction of $2^{2N}$ grid points due to the FFT algorithm employed. To lift these restrictions, we pursue some calculations in a non-uniform, unequal box of arbitrary size other than a cubic shape. For treatment of larger systems, it is highly desirable, as it offers much flexibility. At this stage, we provide some preliminary results on a general non-uniform grid with arbitrary number of points. For this, Cl2 is chosen as a representative in Table VIII. These are performed with LDA XC functional, using a grid spacing of $h_{r}=0.3$, keeping the same pseudopotential. A snapshot of computed energies, from two sets (labelled I, II) is offered for illustration. In Set I in the left-hand side, we first vary $N_{z}$, the number of grid points along internuclear axis ($z$) keeping same along $xy$ plane static at certain reasonable value ($N_{x}=N_{y}=30$). This can be appreciated from a glance at Table I; the calculated total electronic energy of Cl2 is $-$28.41626 a.u. in case of Set A using the LDA XC functional. As $N_{z}$ is increased from 30 to 40, there is a huge jump in energy value ($\approx 4.16$ a.u); with any further increase, it slowly converges to $-$29.66505 a.u., at around $N_{z}=70$. Applying the same procedure for fixed $N_{x},N_{y}$ at 36, offers an energy value of $-$29.70388, converging again nearly at $N_{z}=70$. There is no need to show all these details. Rather it suffices to mention that for each such $(N_{x},N_{y})$ pair, energies attain convergence for a certain $N_{z}$, which in this case happens to be around 70. Then Set II, in the right side, varies the number of points, $N_{x},N_{y}$ along $xy$ plane keeping $N_{z}$ constant at 70. It is apparent that formal convergence of energy takes place for $N_{x}=N_{y}=48$. A re-look at Table I reveals that the comparable convergence in a uniform grid occurs for Set D corresponding to $N_{x}=N_{y}=N_{z}=64$ and $h_{r}=0.3$. This leads to a drastic reduction in the total number of points, which eventually helps to achieve an enhanced computational scaling compared to the uniform grid. This will have a tremendous bearing on routine calculations of larger systems. A more thorough analysis of these aspects will be undertaken later.

IV Future and outlook

We have demonstrated the validity and feasibility of a Gaussian-based LCAO-MO approach to DFT using CCG in the context of atomic and molecular properties. This produces practically identical results with those obtained from other grid-based/grid-free quantum chemistry programs available in literature. The classical Coulomb potential is obtained by means of a Fourier convolution technique, accurately and efficiently in the real-space grid. The matrix elements of all two-body potentials are computed through numerical integration. Quantities such as HOMO energies, potential energy curves, orbital energies, atomization energies are reproduced very well. The method performs quite decently for both local and non-local XC density functionals. While we focused mainly on pseudopotential situations, it is equally applicable for all-electron calculations as well. So far, these works were restricted to a uniform grid; here we report some exploratory calculations for Cl2 in a more general non-uniform grid, for the first time. A more thorough analysis is in progress. Energies are variationally well-founded.

The success of the CCG method, as delineated above, encourages us to employ (i) more appropriate effective core potentials suitable for DFT calculations (ii) elaborate, extended and sophisticated basis sets (iii) superior quality density functionals in our future works. It would be worthwhile to assess its merit and suitability for geometry optimization of molecules using CCG. In this regard, a non-uniform grid, could be quite useful, for which some early results have been presented here and needs further inspection. An important concern would be to lessen the computational cost by incorporating a linear scaling approach; apparently the most promising will be real-space multi-grid techniques. A useful application of this methodology could be in the real-time dynamical studies, particularly laser-atom/molecule interactions in intense/super-intense regimes within a TDDFT framework. Some of these works are currently being investigated by us.

V Acknowledgement

We express our gratitude to the Editors Prof. Michael Springborg and Prof. Jan-Ole Joswig for their kind invitation to contribute in this issue. We also thank them for extending the deadline of submission. Mr. Siladitya Jana provided help in getting some of the references. Computer support from Mr. Suman Chakraborty is duly thanked. AKR acknowledges financial support from DST (Grant No. EMR/2014/000838) and AG from UGC, New Delhi, for a Junior Research Fellowship.