Density functional electric response properties of molecules in Cartesian grid

TL;DR

This paper presents a Cartesian grid-based DFT approach to calculate electric response properties of diatomic molecules, demonstrating accuracy and convergence comparable to traditional atom-centered methods.

Contribution

It develops and validates a Cartesian grid implementation for DFT electric response calculations, including grid convergence analysis and application to molecules at various geometries.

Findings

Excellent agreement with atom-centered-grid results

Accurate response properties for HCl, HBr, HI

Reliable results for molecules far from equilibrium

Abstract

Within the finite-field Kohn-Sham framework, static electric response properties of diatomic molecules are presented. The electronic energy, dipole moment ({\boldmath$\mu$}), static dipole polarizability ({\boldmath$\alpha$}) and first-hyperpolarizability ({\boldmath$\beta$}) are calculated through a pseudopotential-DFT implementation in Cartesian coordinate grid, developed in our laboratory earlier. We engage the Labello-Ferreira-Kurtz (LFK) basis set; while four local and non-local exchange-correlation (LDA, BLYP, PBE and LBVWN) functionals have been adopted. A detailed analysis of \emph{grid convergence} and its impact on obtained results, is presented. In each case the \emph{electric field optimization} was carefully monitored through a recently prescribed technique. For all three molecules (HCl, HBr, HI) considered, the agreement of all these quantities with widely successful and popular atom-centered-grid procedure, is excellent. To assess the efficacy and feasibility, companion calculations are performed for these on a representative molecule (HCl) at distorted geometries, far from equilibrium. Wherever possible, relevant comparison is made with available \emph{all-electron} data and experimental results. This demonstrates that Cartesian grid provides accurate, reliable results for such properties of many-electron systems within pseudopotential representation.

Tables

Table 1. Table 1: Convergence of electronic energy of HCl ( R = 1.275 ​ Å 𝑅 1.275 Å R=1.275\textup{\AA} ) in the grid ( h r = 0.3 subscript ℎ 𝑟 0.3 h_{r}=0.3 ) using LDA XC functional. The reference value is − - 15.433429 schmidt93 . All results are in a.u.

Set I					Set II
$N_x$	$N_y$	$N_z$	$⟨E⟩$		$N_x$	$N_y$	$N_z$	$⟨E⟩$
40	40	40	$-$15.405457		50	50	80	$-$15.432989
-	-	50	$-$15.414922		54	54	-	$-$15.433262
-	-	60	$-$15.415930		58	58	-	$-$15.433361
-	-	70	$-$15.416056		62	62	-	$-$15.433401
-	-	80	$-$15.416061		66	66	-	$-$15.433417
-	-	90	$-$15.416062		70	70	-	$-$15.433424
-	-	100	$-$15.416062		74	74	-	$-$15.433427
50	50	50	$-$15.432747		78	78	-	$-$15.433428
-	-	60	$-$15.432955		90	90	90	$-$15.433429
-	-	70	$-$15.432985		100	100	100	$-$15.433429
-	-	80	$-$15.432989		110	110	110	$-$15.433429
-	-	90	$-$15.432990		120	120	120	$-$15.433429

Table 2. Table 2: Optimal sparsity in grid, in case of HCl for four different XC functionals. The reference values schmidt93 are − - 15.43343, − - 15.47384, − - 15.50509, for LDA, BLYP, PBE respectively. Here h r ≡ h x = h y = h z subscript ℎ 𝑟 subscript ℎ 𝑥 subscript ℎ 𝑦 subscript ℎ 𝑧 h_{r}\equiv h_{x}=h_{y}=h_{z} , and R = 1.275 ​ Å 𝑅 1.275 Å R=1.275\textup{\AA} . All results in a.u.

	LDA				BLYP				PBE				LBVWN
$h_r$	$N_x$	$N_y$	$N_z$	$⟨E⟩$	$N_x$	$N_y$	$N_z$	$⟨E⟩$	$N_x$	$N_y$	$N_z$	$⟨E⟩$	$N_x$	$N_y$	$N_z$	$⟨E⟩$
0.3	74	74	80	$-$15.43343	76	76	80	$-$15.47383	76	76	80	$-$15.50513	68	68	80	$-$15.43287
0.2	106	106	120	$-$15.43343	110	110	120	$-$15.47383	108	108	120	$-$15.50509	98	98	110	$-$15.43287
0.1	196	196	230	$-$15.43342	200	200	230	$-$15.47383	198	198	210	$-$15.50508	184	184	210	$-$15.43286

Table 3. Table 3: Static dipole moment μ z subscript 𝜇 𝑧 \mu_{z} and Finite field 𝜶 𝜶 \bm{\alpha} , α ¯ ¯ 𝛼 \overline{\alpha} and 𝜷 𝜷 \bm{\beta} values (in a.u.) of three diatomic molecules for different XC functionals. PR implies Present Result.

	XC	${\mu }_z$		${\alpha }_{xx=yy}$		${\alpha }_{zz}$		${\overline{\alpha }}^{f,g}$	${\beta }_{xxz=yyz}$		${\beta }_{zzz}$
molecule	functionals	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )	PR	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )
HCl11footnotemark: 1,22footnotemark: 2,33footnotemark: 3	LDA	$-$0.43826	$-$0.43826	18.48	18.48	19.38	19.38	18.79	8.26	8.27	20.77	20.77
	BLYP	$-$0.42337	$-$0.42337	18.19	18.19	19.24	19.24	18.55	6.28	6.28	19.60	19.60
	PBE	$-$0.43420	$-$0.43425	18.05	18.04	19.01	19.01	18.37	7.19	7.19	18.91	19.89
	LBVWN	$-$0.45357	–	15.39	–	17.41	–	16.07	3.77	–	15.20	–
HBr44footnotemark: 4,55footnotemark: 5	LDA	$-$0.31611	$-$0.31612	25.33	25.32	26.58	26.58	25.74	7.27	7.25	23.13	23.12
	BLYP	$-$0.29881	$-$0.29882	24.71	24.70	26.16	26.16	25.19	4.16	4.16	20.90	20.90
	PBE	$-$0.31207	$-$0.31208	24.63	24.64	25.99	25.99	25.08	5.79	5.74	20.51	20.49
	LBVWN	$-$0.30761	–	21.31	–	24.04	–	22.22	2.19	–	15.60	–
HI	LDA	0.18225	0.18226	37.38	37.37	39.13	39.13	37.96	$-$3.11	$-$3.17	$-$16.47	$-$16.36
	BLYP	0.16093	0.16103	36.46	36.48	38.34	38.33	37.08	1.21	1.21	$-$12.66	$-$12.63
	PBE	0.17943	0.17945	36.35	36.33	38.18	38.18	39.24	$-$1.84	$-$1.73	$-$13.25	$-$13.19
	LBVWN	0.11947	–	32.15	–	35.74	–	33.34	$-$2.44	–	$-$9.18	–

Table 4. Table 4: Static dipole moment μ z subscript 𝜇 𝑧 \mu_{z} , FF 𝜶 𝜶 \bm{\alpha} , α ¯ ¯ 𝛼 \overline{\alpha} and 𝜷 𝜷 \bm{\beta} values (in a.u.) of HCl molecule at various distorted geometries. All quantities are in a.u. PR implies Present Result.

$R$	XC	${\mu }_z$		${\alpha }_{xx=yy}$		${\alpha }_{zz}$		${\beta }_{xxz=yyz}$		${\beta }_{zzz}$
	functional	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )	PR	Ref. (schmidt93, )
1.5	LDA	$-$0.31111	$-$0.31111	16.80	16.80	13.85	13.85	20.29	20.29	18.59	18.59
	BLYP	$-$0.31210	$-$0.31213	16.69	16.69	13.64	13.64	19.20	19.19	16.55	16.52
	PBE	$-$0.32008	$-$0.32013	16.46	16.46	13.52	13.52	19.37	19.35	17.05	16.89
2.0	LDA	$-$0.37343	$-$0.37344	17.71	17.71	16.38	16.38	15.24	15.25	20.66	20.67
	BLYP	$-$0.36844	$-$0.36844	17.50	17.50	16.20	16.20	13.68	13.67	19.14	19.15
	PBE	$-$0.37698	$-$0.37703	17.14	17.14	15.42	15.42	15.56	15.56	18.66	18.63
2.5	LDA	$-$0.45428	$-$0.45429	18.67	18.67	20.19	20.19	6.48	6.48	20.98	20.99
	BLYP	$-$0.43630	$-$0.43630	18.36	18.35	20.06	20.06	4.39	4.40	19.84	19.83
	PBE	$-$0.44800	$-$0.44804	18.22	18.21	19.81	19.81	5.39	5.39	19.07	19.07
3.0	LDA	$-$0.55506	$-$0.55506	19.63	19.63	25.44	25.44	$-$3.83	$-$3.82	26.01	26.01
	BLYP	$-$0.51364	$-$0.51361	19.22	19.21	25.42	25.42	$-$6.49	$-$6.47	24.69	24.68
	PBE	$-$0.53309	$-$0.53312	19.13	19.13	25.05	25.05	$-$4.98	$-$4.99	23.25	23.26

Equations

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

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

$$\rho (\mathbf{r})={\rho }^{\alpha }(\mathbf{r})+{\rho }^{\beta }(\mathbf{r})=\sum _i{f}_i^{\alpha }|{\psi }_i^{\alpha }(\mathbf{r}){|}^2+\sum _i{f}_i^{\beta }|{\psi }_i^{\beta }(\mathbf{r}){|}^2,$$

$$\rho (\mathbf{r})=\sum _{\mu ,\nu }^{{\mathrm{N}}_{\mathrm{c}}}{P}_{\mu \nu }{\chi }_{\mu }^{\star }(\mathbf{r}){\chi }_{\nu }(\mathbf{r}),$$

$$P_{\mu \nu }=P_{\mu \nu }^{\alpha }+P_{\mu \nu }^{\beta },P_{\mu \nu }^{\sigma }=\sum _i{f}_i^{\sigma }{C}_{\mu i}^{\star \sigma }{C}_{\nu i}^{\sigma },$$

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

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

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

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

$$v_h({\mathbf{r}}_g)=FF{T}^{-1}\{v_h({\mathbf{k}}_g)\rho ({\mathbf{k}}_g)\}\mathrm{and}\rho ({\mathbf{k}}_g)=FFT\{\rho ({\mathbf{r}}_g)\}.$$

$$v_h({\mathbf{r}}_g)=\frac{\mathrm{erf}(\zeta r)}{r}+\frac{\mathrm{erfc}(\zeta r)}{r}\equiv v_{h_{\mathrm{long}}}({\mathbf{r}}_g)+v_{h_{\mathrm{short}}}({\mathbf{r}}_g),$$

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

$${\mu }_i(\mathbf{F})={\mu }_i(0)+\sum _j{\alpha }_{ij}{F}_j+\frac{1}{2}\sum _{j,k}{\beta }_{ijk}{F}_j{F}_k+\cdots .$$

$$\left.\begin{aligned} \alpha_{ii}F_{i}=\frac{2}{3}\bigg{[}\mu_{i}(F_{i})-\mu_{i}(-F_{i})\bigg{]}-\frac{1}{2}\bigg{[}\mu_{i}(2F_{i})-\mu_{i}(-2F_{i})\bigg{]}\\ \alpha_{ij}F_{j}=\frac{2}{3}\bigg{[}\mu_{i}(F_{j})-\mu_{i}(-F_{j})\bigg{]}-\frac{1}{2}\bigg{[}\mu_{i}(2F_{j})-\mu_{i}(-2F_{j})\bigg{]}\\ \beta{iii}F_{i}^{2}=\frac{1}{3}\bigg{[}\mu_{i}(2F_{i})+\mu_{i}(-2F_{i})\bigg{]}-\frac{1}{3}\bigg{[}\mu_{i}(F_{i})+\mu_{i}(-F_{i})\bigg{]}\\ \beta{ijj}F_{j}^{2}=\frac{1}{3}\bigg{[}\mu_{i}(2F_{j})+\mu_{i}(-2F_{j})\bigg{]}-\frac{1}{3}\bigg{[}\mu_{i}(F_{j})+\mu_{i}(-F_{j})\bigg{]}\end{aligned}\right\}$$

$$\bar{\alpha }=\frac{1}{3}({\alpha }_{zz}+{\alpha }_{xx}+{\alpha }_{yy})$$

$$H_{\mu \nu }^{{}^{\prime }\mathrm{core}}=H_{\mu \nu }^{\mathrm{core}}+F_i⟨\mu |\mathbf{r}|\nu ⟩,i\in \{x,y,z\}.$$

$$\mathbf{\mu }\equiv {\mathbf{\mu }}_{el}=\sum _{\mu \nu }{P}_{\mu \nu }⟨\mu |\mathbf{r}|\nu ⟩+\sum _a{Z}_a{\mathbf{R}}_a,$$

$$F_n=F_0\times 2^{\frac{3\mathrm{n}}{100}}$$

$${\zeta }_{\mathbf{\beta }}=\frac{{\mathbf{\beta }}_{FF}}{{\mathbf{\beta }}_{FF}^{ref}}-1.$$

Full text

Keywords: Density functional theory, electric response, dipole moment, polarizability, first-hyperpolarizability, cartesian grid.

I Introduction

Density functional theory (DFT) (hohenberg64, ; kohn65, ) has been found to be an indispensable tool for determining structure and properties of many-electron systems like atoms, molecules, clusters, nano-materials and periodic systems over the past three decades parr89 ; joubert98 ; gidopoulos03 ; cohen11 ; engel11 ; burke12 ; becke14 ; jones15 . The theory is exact, in principle and has the ability to capture many-body correlation effects at a computational cost comparable to Hartree-Fock theory. Naturally, it has been widely applied in solid-state physics and material chemistry for understanding the physio-chemical phenomena. The linear and non-linear electric properties such as dipole moment, polarizability, hyperpolarizability are important for many applications, e.g., the development of nonlinear optical materials munn01 , Raman and infrared spectroscopy tanimura09 , structural identifications of atomic clusters (wang09, ), separations of molecular isomers (kinnibrugh06, ), etc. Throughout the past several decades, a lot of theoretical developments have taken place to determine the electric response properties of atoms and molecules using DFT as well as different ab initio methods. Several excellent reviews are available; here we refer to a selective set (pugh99, ; champagne10, ; helgaker12, ).

The above-mentioned properties may be computed using a number of different theoretical approaches within the Kohn-Sham (KS) DFT rubric. Broadly speaking, there exists two distinct routes. In the first case, one obtains the response of density matrix analytically by solving a set of coupled perturbed KS (CPKS) equations (fournier90, ; colwell93, ; komornicki93, ; lee94, ). The atomic-orbital basis along with the iterative nature of this procedure makes it usually somehow expensive in terms of computational overhead for larger systems. Moreover, this also requires information about analytical gradients or excited states. Two major prescriptions have been put forth in the literature, namely (i) a non-iterative approximation to CPKS (sophy03, ; sophy08, ) and (ii) an auxiliary density perturbation theory (moreno08, ; carmona12, ). The former is semi-numerical in nature where the derivative of KS matrix is estimated using finite-field (FF) approximation from which the response density matrix is calculated. The latter brings in a set of auxiliary functions to express the approximated density to calculate the perturbed density matrix. Further, time-dependent density functional theory (TDDFT) was invoked in the same spirit as CPKS method for static as well as dynamic (frequency dependent) electric responses of ground, and subsequently excited states of molecule (casida95, ; gorling99, ; hu02, ). The advantage lies in the fact that it has the ability to determine nonlinear optical properties such as, frequency-dependent hyperpolarizabilities, second harmonic generation, and etc gisbergen98 . Besides these methods, which involve variations of density matrix with respect to electric field, there exists another attractive approach, called perturbative sum-over states expression over all dipole-allowed electronic transitions. Here, time-independent/time-dependent perturbation theory is used to identify the expressions for static/dynamic response functions respectively, by expanding the energy in order of external perturbations (orr71, ; bishop94, ).

Another important direction is the FF method, which is a common technique in quantum chemistry for calculating these properties, as reflected from a large volume of works recorded (jasien90, ; guan93, ; calaminici98, ; kamada00, ; bulat05, ; touma11, ). Its real advantage lies in the ease of implementation over the analytical response theory. In general, the response properties are very sensitive towards basis set, electron correlation, relativistic effects and vibrational structure, in case of molecules. Note that, as the response largely results from valence electrons, sufficient diffuse functions must be included in basis set. A very careful consideration of polarization functions is an important requisite (miadokova97, ) to obtain accurate values. Furthermore, one could use pseudopotential approximation to reduce the computational burden substantially instead of full calculations. In that case, their determination depends on the accuracy and transferability of pseudopotential (filippetti95, ). Besides, valence basis sets which are used for pseudopotential calculations, also suffer from similar shortcomings as standard all-electron basis sets (labello05, ). Fortunately, in many cases response of core electrons is small compared to that of valence electrons, which makes the calculated results acceptable and realistic. In yet another effort, basis-set free methods offered a practical alternative, as they do not suffer from the incompleteness of basis set. For example, the real-space approach of (vasiliev10, ) provides quite comparable results for atomic and molecular polarizabilities within pseudopotential FF DFT framework. Another interesting development for polarizability was put forth by Talman (talman12, ), by advocating the use of Sternheimer method within the unoccupied HF and coupled perturbed HF approximations.

In this communication, we report the implementation of FF method within the frame work of LCAO-MO pseudopotential KS formalism in Cartesian Coordinate grid (CCG) (roy08, ; roy08a, ; roy10, ; roy11, ; ghosal16, ). The scope and applicability of the method was demonstrated for a decent number of atoms/molecules for various exchange-correlation (XC) functionals in terms of energy, ionization energy, atomization energy, orbital energy and potential energy curves. Some preliminary results on non-uniform grid was published as well ghosal16 . Here we investigate molecular properties, namely, dipole moment ($\mu$), dipole polarizability ($\alpha$) and first-hyperpolarizability ($\beta$), computed explicitly including the static electric field contribution in the time-independent KS Hamiltonian. These are done for three molecules (HCl, HBr, HI), taking four XC functional, viz., (i) local density approximation (LDA) with the Vosko-Wilk-Nusair (VWN) (vosko80, ) correlation (ii) Becke becke88 exchange and Lee-Yang-Parr (LYP) lee88 correlation–abbreviated as BLYP (iii) Perdew-Burke-Ernzerhof (PBE) (perdew96, ) functional (iv) asymptotically corrected Leeuwen-Baerends (LB94) (leeuwen94, ) exchange plus VWN correlation–abbreviated as LBVWN. The last one is considered because its superior asymptotic long-range behavior ensures the ionization energies (obtained through HOMO) as well as higher-lying states, quite satisfactorily. Throughout our presentation, Labello-Ferreira-Kurtz (LFK) basis set labello05 was used, which appears to be a quite appropriate in the literature for such pseudopotential studies. In each case, both spatial and field grid optimization was performed quite carefully. We also thoroughly examine the role of XC functionals to determine the optimal FF strength in case of a representative test molecule (HCl). The computed quantities are compared with available theoretical data in order to evaluate the effectiveness of CCG in this context. Additionally, the average polarizability for these molecules is also reported to facilitate comparison with existing experimental and theoretical results. Next, $\bm{\alpha}$ is probed in our test molecule HCl, at different internuclear separations ($R$), for which static correlation plays a crucial role. The suitability and applicability of our results is then correlated with respective all-electron calculations using Sadlej (sadlej92, ) basis set, in the range of $R$ covered in this work.

II Method of Calculation

The single-particle KS equation of a many-electron system under the influence of pseudopotential can be written as (atomic unit employed unless stated otherwise),

[TABLE]

where $v^{p}_{ion}$ denotes the ionic pseudopotential (denoted by a “p” superscript), expressed as,

[TABLE]

Here $v^{p}_{ion,a}(\mathrm{\mathbf{r}})$ represents the ion-core pseudopotential associated with atom A, located at $\mathrm{\mathbf{R}}_{a}$. The classical Coulomb (Hartree) potential, $v_{h}[\rho(\mathrm{\mathbf{r}})]$ describes the electrostatic interaction among the valence electrons whereas $v_{xc}[\rho(\mathrm{\mathbf{r}})]$ signifies XC potential, the non-classical part of effective Hamiltonian. The single-particle charge density is then given by,

[TABLE]

where $\{\psi^{\sigma}_{i},\sigma=\alpha\quad\mathrm{or}\quad\beta\}$ corresponds to a set of $N$ occupied orthonormal spin molecular orbitals (MO) and $f_{i}^{\sigma}$s denote occupation numbers related to ith spin-MO. Alternatively, in terms of the atom-centered basis functions $\{\chi_{\mu}(\mathrm{\mathbf{r}})\}$, one can write,

[TABLE]

with

[TABLE]

where $P_{\mu\nu}$ denotes an element of one-body density matrix $\mathrm{\mathbf{P}}$ while $\mathrm{N_{c}}$ corresponds to the total number of contracted basis functions. In LCAO-MO approximation, the coefficients for expansion of spin-MOs satisfy a set of equations analogous to that in Hartree-Fock theory,

[TABLE]

with the orthonormality condition, $(\mathrm{\mathbf{C}}^{\sigma})^{\dagger}\mathrm{\mathbf{S}}\mathrm{\mathbf{C}}^{\sigma}=\mathrm{\mathbf{I}}$. Here $\mathrm{\mathbf{C}}^{\sigma}$ contains the respective spin-MO coefficients $\{C_{\nu i}^{\sigma}\}$, $\mathrm{\mathbf{S}}$ is the overlap matrix corresponding to elements $S_{\mu\nu}$, $\bm{\epsilon}^{\sigma}$ refers to diagonal matrix of respective spin-MO eigenvalues $\{\epsilon_{i}^{\sigma}\}$, while $F_{\mu\nu}^{\sigma}$ is an element of KS-spin matrix conveniently partitioned as,

[TABLE]

In this equation, $H_{\mu\nu}^{\mathrm{core}}$ contains all one-electron contributions including kinetic energy, nuclear-electron attraction and pseudopotential matrix elements. All one-electron integrals are generated by standard recursion relations (obara86, ) using Cartesian Gaussian-type orbitals as primitives basis functions. We employ the angular-momentum dependent pseudopotential form as proposed by (stevens84, ; stevens92, ), whereas $J_{\mu\nu}$ term signifies contribution from Hartree potential and the last (XC) term arises from non-classical effects.

Now all the relevant quantities like basis functions, electron densities, MOs as well as various two-electron potentials are directly set up on the $3$D CCG simulating a cubic box,

[TABLE]

where $h_{r}$ denotes grid spacing along each directions and $N_{x},N_{y},N_{z}$ signify total number of grid points along $x,y,z$ directions respectively. In case of non-uniform grid, we usually vary $N_{x},N_{y},N_{z}$ independently keeping the value of $h_{r}$ fixed. Thus, electron density $\rho(\mathrm{\mathbf{r}})$ in this grid may be simply written as (“g” implies the discretized grid),

[TABLE]

A major concern in grid-based approach constitutes an accurate estimation of classical electrostatic potential. Here, we invoke the conventional Fourier convolution method (martyna99, ; minary02, ). It has been well documented earlier (roy08, ; roy08a, ; roy10, ; roy11, ; ghosal16, ); so here only the essential details are summarized. In the end, the classical Coulomb potential is calculated from the following,

[TABLE]

The quantities $v_{h}(\mathrm{\mathbf{k}}_{g}),\rho(\mathrm{\mathbf{k}}_{g})$ stand for Fourier integrals of Coulomb interaction kernel and density respectively. The electron density in $k$ space, $\rho(\mathrm{\mathbf{k}}_{g})$ can be obtained easily using discrete Fourier transform of respective real-space values. Thus the real crux of our problem is calculation of Coulomb interaction kernel, which has singularity at origin. In order to overcome this, we exploit an Ewald summation-type approach (chang12, ), expanding the Coulomb kernel into long- and short-range components,

[TABLE]

where erf(x) and erfc(x) denote error function and its complement respectively. Fourier transform of short-range part can be treated analytically whereas the long-range portion needs to be computed directly from FFT of corresponding real-space values. A convergence parameter $\zeta$ is used to adjust the range of $v_{h_{\mathrm{short}}}(\mathrm{\mathbf{r}}_{g})$, such that the error is minimized.

A very crucial step in DFT calculations is choice of an appropriate XC functional. While the exact form remains elusive as yet, highly accurate functionals have been reported including so-called local, “non-local” (gradient and Laplacian-dependent) and hybrid ones. Finally, all the two-electron KS matrix elements are calculated directly through numerical integration on the grid as ($v_{hxc}$ refers to Hartree and XC potential combined),

[TABLE]

The response properties of a many-electron system can be defined by expanding field-dependent dipole moment, calculated from the field-induced charge distribution, as a power series in the external electric field $\mathrm{\mathbf{F}}$, if the field strength remains small, i.e.,

[TABLE]

In this equation, the three terms on right-hand side characterize static dipole moment $\mu_{i}(0)$, dipole polarizability $\alpha_{ij}=\frac{\partial{\mu_{i}}}{\partial F_{j}}$ and first-hyperpolarizability $\beta_{ijk}=\frac{\partial^{2}\mu_{i}}{\partial F_{j}\partial F_{k}}$ (mclean67, ) respectively. Alternatively one may also represent it in terms of field-induced energy; and both the definitions are equivalent according to Hellmann-Feynman theorem (feynman39, ). The components of $\bm{\alpha}$ and $\bm{\beta}$ can be deduced from the following well-known finite-difference formulas (smith78, ) as,

[TABLE]

Moreover, in addition to $\bm{\alpha},\bm{\beta}$ tensors, the experimentally determined so-called average polarizability, defined as,

[TABLE]

can also be calculated for a given species. Now in order to obtain all these above mentioned tensors from dipole moment of the system (expressed as a function of external electric field $\mathrm{\mathbf{F}}$), one needs the perturbed density matrix at different field strengths, which is obtained from the self-consistent solution of Eq. (1). Hence the core part of the Hamiltonian (denoted by a prime) will now be modified by a field-dependent term accordingly as,

[TABLE]

Here $H_{\mu\nu}^{\mathrm{core}}$ refers to unperturbed core Hamiltonian described above, $F_{i}$ denotes ith component of applied field $\mathrm{\mathbf{F}}$ and $\langle\mu|\mathrm{\mathbf{r}}|\nu\rangle$ gives the dipole moment integral corresponding to length vector $\mathrm{\mathbf{r}}$. All two-body matrix elements of KS matrix will remain intact during FF calculations. Finally, dipole moment of a molecule can be expressed as below,

[TABLE]

where $Z_{a}$ and $\mathrm{\mathbf{R}}_{a}$ are nuclear charge and position of atom “a”, respectively.

III Computational aspects

This section provides some of the computational details in our current calculation. Methods like the present one, require the choice of a suitable basis set that can correctly describe the response of electrons to an external perturbation. It is quite well known that it should typically contain diffuse functions to provide accurate results (werner76, ). While there exist several options for full calculations which contain $d,f$ orbitals, the choice is much limited for pseudopotential approximations. In the current work, we have adopted the so-called LFK basis as proposed in (labello05, ), based on a procedure to incorporate diffuse and polarization functions in familiar Sadlej basis set (sadlej92, ). It was been pointed out that response properties (mainly $\bm{\alpha}$) comparable to all-electron results (using Sadlej basis) can be recovered at a lower computational cost using this basis. These are taken from EMSL Basis Set Library (feller96, ).

The pertinent molecular properties are calculated for four different XC functionals: (i) LDA–with the homogeneous electron gas correlation proposed by VWN (parametrization formula V of Ref. (vosko80, )) (ii) BLYP–incorporating the popular Becke (becke88a, ) exchange along with LYP (lee88, ) correlation (iii) PBE (perdew96, ) functional (iv) LBVWN–including LB94 exchange (leeuwen94, ) along with VWN correlation. All XC functionals were adopted from density functional repository program (repository, ) except LDA and LB94.

The self-consistent convergence criteria imposed in this communication is slightly tighter than our earlier work (roy08, ; roy08a, ; roy10, ; roy11, ; ghosal16, ); this is to generate a more accurate perturbed density matrix. Changes in following quantities were checked, viz., (i) orbital energy difference between two successive iterations and (ii) absolute deviation in a density matrix element. They both were required to remain below a certain prescribed threshold set to $10^{-8}$ a.u.; this ensured that the total energy maintained a convergence of at least this much. To accelerate FF calculations, unperturbed (field-free) density matrix was used as trial input. The value of $\zeta$ in Eq. (11) is fixed in such a way that $\zeta\times L=7$, where $L$ is chosen as the smallest side of simulating box. Such a conjecture was put forth in martyna99 and quite successfully implemented in CCG before (ghosal16, ). In order to perform discrete Fourier transform, standard FFTW3 package (fftw05, ) is invoked. The resulting generalized matrix-eigenvalue problem is solved through standard LAPACK routine (anderson99, ) accurately and efficiently. Relevant pseudopotential matrix elements in Gaussian orbitals are imported from GAMESS (schmidt93, ) package. The scaling properties have been discussed earlier roy08a .

IV Results and discussion

At first, it may be worthwhile to discuss the influence of spatial grid on total energy, $\langle E\rangle$ of a representative system with respect to the sparsity of grid (regulated by $N_{x},N_{y},N_{z}$) and grid spacing (determined by $h_{r}$). For this, we consider a closed shell molecule, such as hydrogen chloride (HCl) at the experimental bond length of $1.275$Å along $z$ axis, taken from NIST database (johnson16, ) as test case. The formal convergence and stability of the grid is illustrated in Table I for LDA XC functional at a grid spacing of $h_{r}=0.3$. As in (ghosal16, ), we first vary $N_{z}$, the number of grid points along internuclear axis, keeping the same along $xy$ plane static at certain reasonable value, say $N_{x}=N_{y}=40$. A glance at this table shows that as $N_{z}$ is gradually increased from 40 to 90 with an increment of 10, there is a smooth convergence in energy at around $N_{z}=80$ with a difference in total energy (or more specifically grid accuracy) of about $5\times 10^{-6}$ a.u. between two successive steps. In the beginning, when $N_{z}$ moves from 40 to 50 to 60, one notices quite dramatic improvement in energy; but after that the changes are relatively less until eventually reaching convergence for $N_{z}$ at around 80. Now, applying the same procedure for fixed values of $(N_{x},N_{y})$, say at 50, offers an energy value of $-$15.432989, converging again in the same neighborhood of $N_{z}\approx 80$ with grid accuracy $5\times 10^{-6}$ a.u. This trend is maintained for other $N_{x},N_{y}$ pairs; which is not presented here to save space. It is clear from Set I that for each such $(N_{x}$, $N_{y})$ pair and a particular grid accuracy, energy attains convergence for nearly the same value of $N_{z}$, which in this case happens to be around 80. Then in Set II in right hand side, we vary $N_{x},N_{y}$ along $xy$ plane keeping $N_{z}$ fixed at $80$. It is apparent from Set II that the convergence in energy takes place at $N_{x}=N_{y}=74$ with same grid accuracy of Set I. As a further check on the numerical stability in our solution, several additional calculations were performed in much extended grids (last five entries); as anticipated energy remains stable in all such grids. Besides, it is also verified that different increments in $N_{x},N_{y},N_{z}$ do not bring any significant change in the optimal number of grid points as estimated above. Repeating same steps for other functionals gave only slightly different $(N_{x},N_{y},N_{z})$ triplet for desired energy convergence, which is detailed as below.

Now the effect of grid spacing on energy convergence is analyzed for three different $h_{r}$ ($0.3$, $0.2$ and $0.1$) for all four XC functionals, namely LDA, BLYP, PBE, LBVWN. This is accomplished by following the simple grid optimization strategy as delineated above, maintaining a grid accuracy of $5\times 10^{-6}$ a.u., all throughout. The final converged energies with respect to various spacings for these functionals are offered in Table II. Clearly for a fixed $h_{r}$, the optimum $N_{x},N_{y},N_{z}$ only marginally vary from functional to functional. Also it is evident that to produce similar-quality results in a dense grid (smaller $h_{r}$) typically requires larger numbers of grid points. The above discussion thus suggests that a $h_{r}$ of $0.3$ and an optimal grid of ($74,74,80$) is sufficiently good for all practical purposes in our test molecule.

Next, we move towards the electric response properties mainly $\bm{\alpha}$ and $\bm{\beta}$ tensors. It is now well established that the FF procedure significantly depends on choice of an appropriate field strength and its distribution. Two opposite and counter-intuitive effects govern the overall behavior of response, and they both should be satisfied at same time. First, the field must be sufficiently large so that it can overcome the finite-precision artifacts when one numerically differentiates the dipole moment with respect to electric field. On the other hand, it must also be small enough so that the contributions from remaining higher-order derivatives become negligible. In order to optimize the above mentioned parameters, Taylor, polynomial or rational function-based FF methods (mohammed13, ; mohammed17, ; patel17, ) have been implemented quite successfully in the literature. It is also well documented that the effect of field is rather quite delicate, especially in case of higher-order derivatives, such as $\bm{\beta}$ and $\bm{\gamma}$. The problem is critical, for the numerical stability of the latter derivatives satisfies a rather narrow range of suitable field strengths. So our plan is to vary electric field for a given system to determine $F_{opt}$ such that this can be used in general if possible, irrespective of the XC functional involved. Towards this direction, we have employed a fine field as prescribed in (mohammed17, ); accordingly the calculations are performed at discrete field strengths as given by the following equation,

[TABLE]

with $n$ ranging from 0-160, and $F_{0}=0.0005$ a.u. This yields a maximum field strength greater than $0.01$ a.u. At each $F_{n}$, the required properties are calculated using Eq. (14) for a fixed grid and basis set as mentioned above.

Following Buckingham (buckingham67, ), there are two independent components associated with $\bm{\alpha}\ (\alpha_{xx=yy},\alpha_{zz})$ and $\bm{\beta}\ (\beta_{xxz=yyz},\beta_{zzz})$ respectively, for a heteronuclear diatomic molecule belonging to $C_{\infty v}$ symmetry. The maximum field response towards the electron density is then found along $z$ direction as it is the molecular axis. So we choose only $\alpha_{zz}$ and $\beta_{zzz}$ among these, and show the effect of field strength on these, in Fig. 1 using three functionals (LDA, BLYP, PBE). From panel (a), one notices that, $\alpha_{zz}$ is practically independent of field strength for a broad region ranging from 0.0005 to 0.01 a.u.; this holds true all three functionals. This is in keeping with the previous work of (calaminici98, ) along this direction. Therefore one is free to choose a suitable field strength in above region for calculation of $\bm{\alpha}$ for all these functionals. In this way, a field strength of $10^{-3}$ a.u. is used for $\bm{\alpha}$ for HCl for all the functionals. This technique is quite general, and has been applied to other systems considered in this work as well. Now, a similar kind of analysis was performed for $\beta_{zzz}$ in panel (b), again in case of HCl, for same three XC functionals as above. As expected, generally the changes with respect to field remains rather small in low field strength, and tends to grow in the higher range. A careful analysis reveals that, for all three functionals, $\beta_{zzz}$ values show a maximum deviation of up to 0.02 a.u., for a field strength ranging up to about 0.004 a.u.

In order to analyze this effect in some detail, a field-dependent parameter $\zeta_{\bm{\beta}}$ is introduced which is a pointer to the relative error in $\bm{\beta}_{FF}$ with respect to a standard reference value, the minimum of which corresponds to $F_{opt}$ (mohammed17, ). This is defined as,

[TABLE]

Without any loss of generality, the FF result obtained from standard GAMESS package (schmidt93, ) is chosen as our reference. It may be noted that the primary objective of current methodology is to establish the suitability of this grid for molecular properties calculations. Accordingly, we have performed the GAMESS calculations with default options: 96 radial and 302 angular points for the spatial grid, and 0.001 for the field strength. In this context, it is interesting to note that, a recent study of grid effects (based on atom-centered grid), reported by Castet et al. (castet12, ), suggested a grid of 99 radial and 974 angular points, to be an optimally good solution for such calculations. We have verified that for all three molecules considered here, the default option delivers results which are practically coincident with that from the finer grid. Thus the default grid serves our current purpose, as the main concern here is to validate the present grid with traditional grid. So, taking into account this field-dependent parameter $\zeta$, the field variations for all functionals except LBVWN (for which the reference results are unavailable), are displayed in panel (c) of Fig. 1. One finds minima in case of both LDA and BLYP corresponding to the desired $F_{opt}$, but they both lie within a very narrow range of field strength from $10^{-3}$ to $3\times 10^{-3}$ in a.u. However, PBE result does not show any such minimum; the plot rather gradually increases. Moreover the relative error ($\zeta_{\bm{\beta}}$) in PBE plot always remains above both LDA or BLYP values for the entire field strength considered. Therefore it does not appear straightforward to discern an $F_{opt}$ for $\bm{\beta}$ in a given electronic system, whatever be the functional used. Clearly this occurs due to the narrow range of stability with respect to electric fields. On the light of above facts, calculation of $\bm{\beta}$ for our test molecule is performed with $F_{opt}$ corresponding to the minima in (c) for LDA, BLYP functionals. However, for remaining two functionals, these are chosen in the neighborhood of $1\times 10^{-3}$ a.u.

In order to extend the scope and applicability of current scheme, we now report the non-zero components of FF $\bm{\alpha}$, $\bm{\beta}$ along with $\overline{\alpha}$ and static dipole moment $\mu$, of two selective molecules (HBr, HI) along with the test molecule in Table III. These are provided for all four XC functionals. To put things in perspective, we also quote the reference values (except LBVWN) obtained from GAMESS software (schmidt93, ). Same strategy as in HCl has been followed for the other two–experimental geometries are taken from NIST database (johnson16, ), active grid for each molecule was optimized accordingly, keeping the molecular axis along z direction to be fixed, and then followed up the FF procedure as mentioned earlier. The fact that $\bm{\alpha}$ remains stable for a relatively broad range of field strength, enables us to fix the $F_{opt}$ safely at $10^{-3}$ for all these molecules (for all functionals). The same field optimization protocol as in HCl was applied for all of them to get an idea about $F_{opt}$, for $\bm{\beta}$. The maximum absolute deviation (MAD) in $\mu_{z}$ in HCl is around $5\times 10^{-5}$ a.u. for PBE, whereas the same match perfectly with reference (in all the digits quoted) for LDA and BLYP. The same in case of HBr, HI agree quite nicely with reference–MAD for LDA and PBE being 0.00001 and 0.00002 a.u. respectively, while for BLYP it is 0.0001 (in HI). The $\bm{\alpha},\bm{\beta}$ tensors of our calculations are also equally consistent with the reference data, the respective MAD’s in $\alpha_{xx=yy}$ and $\alpha_{zz}$ being 0.02 and 0.01 respectively (both occur for HI). For $\beta_{xxz=yyz}$ and $\beta_{zzz}$, respective MAD’s turn out to be 0.11 in both occasions (again for HI). We also quote some relevant theoretical results for HCl and HBr in the footnote, along with the methods (such as higher-order perturbation theory, MCSCF, CCSD(T) etc.) and basis set. Similarly, experimental values for $\overline{\alpha}$ are also recorded from two different kinds of experimental techniques (olney97, ; hohm13, ). These values contain only electronic part and neither geometry relaxation in presence of electric field, nor vibrational contribution are considered. It reveals an interesting fact that, all three traditional functionals (LDA, BLYP, PBE) overestimate both experimental results in the range of $3$-$11\%$, $1$-$9\%$ and $3$-$8\%$ respectively; however, LBVWN underestimates by $5$-$9\%$. Similar conclusions have also been drawn regarding the pattern behavior of these functionals for $\overline{\alpha}$ in (cohen99, ; vasiliev10, ), where it has been conjectured that, a significant improvement may be achieved for $\bm{\alpha}$ by combining LB94 potential (in asymptotic region) with LDA (corrected for derivative discontinuity in the bulk region) exchange suitably. This leads a more accurate representation of exchange which approaches the experimental results quite closely (casida98, ; casida00, ) at a lower level of computational cost than standard XC functionals.

In the last part of this study, we investigate the efficacy of CCG in determining the non-zero components of $\mu$ as well as $\bm{\alpha}$ and $\bm{\beta}$ tensors at different internuclear separations ($R$) in Table IV. As an illustration, once again, HCl has been chosen with $R$ ranging from 1.5-3.0 a.u. In general, beyond equilibrium geometry the static correlation becomes dominant; hence, the role of XC functional is of utmost importance. Moreover, the role of basis set is also a major factor in the estimation of $\bm{\alpha}$, $\bm{\beta}$ in such regions, and it will be discussed shortly in next paragraph. The optimal spatial mesh is determined at each $R$ using the automated grid optimization technique of Table I. The same optimal field strengths of Table III were adopted, as these do not affect the qualitative nature of present results with respect to $R$. All relevant quantities are recorded in Table IV at four distinct (1.5, 2, 2,5, 3) $R$. The computed $\mu_{z}$ values are in very excellent agreement with theoretical references, for all XC functionals throughout the whole region; maximum discrepancy ($5\times 10^{-5}$) occurs in case of PBE at $R=2$ a.u. A similar comparison for $\alpha_{zz}$ and $\alpha_{xx=yy}$ reveals that, in the former, reference results are completely reproduced by CCG, again for all functionals for all $R$’s considered, while in the latter case, the two results remain separated by a MAD of 0.01 in few occasions. Similarly, the agreement of $\bm{\beta}$ components with reference is also excellent for all $R$’s with MAD being 0.16 a.u., for the lone case of PBE at $R=1.5$. It may occur due to the fact that $F_{opt}$ is in general, affected by molecular size in case of $\bm{\beta}$, and we have not taken this into account (same $F_{opt}$ is tacitly assumed to be valid for all $R$, which may not hold true). A closer look at this table further reveals that there is a change in sign in $\beta_{xxz=yyz}$ on varying $R$ from 2.5 to 3 a.u., which is quite satisfactorily captured in our results. Lastly to conclude this portion, no attempt was made to do field optimization at each $R$. While this has little or practically no bearing on $\alpha$, the estimation of $F_{opt}$ may have appreciable effect in estimating $\bm{\beta}$ at different distorted geometries.

Next we examine the impact of XC functional and basis set on changes in $\bm{\alpha}$ at distorted geometry (of same HCl molecule) and compare with all-electron results. As far as our knowledge goes, such an analysis of the performance of LFK basis in determining the response properties beyond equilibrium geometry $(R_{eq})$, has not been undertaken before. Figure 2 portrays the calculated $\alpha_{xx=yy}$ and $\alpha_{zz}$ (segments (a), (b) respectively) with $R$, covering a broad region of 1.5-3 a.u., at intervals of $0.1$ a.u. The all-electron results are done using Sadlej basis (sadlej92, ) and standard B3LYP functional through the GAMESS program. All the functionals reproduce the qualitative shape of $\alpha_{xx=yy}$ and $\alpha_{zz}$ very well for the entire range. It is noticed that in both panels, PBE results are the closest to Sadlej-B3LYP results around $R_{eq}$. While in panel (a) all the plots remain well separated, a distinct crossover is noticed in panel (b) as one moves farther beyond $R_{eq}$. On closer inspection, PBE plot in (b) tends to deviate maximum from the all-electron results amongst all functionals. Thus our present pseudopotential FF DFT calculation (with the aid of LFK basis) in CCG can produce comparable results for $\alpha_{xx=yy}$, $\alpha_{zz}$ with the more elaborate full calculations, both around and far from equilibrium.

A few remarks may now be made before passing. In this work our primary objective was to demonstrate that the real-space CCG coupled with FF method (for electric response calculations) could deliver accurate and physically meaningful results within the chemical accuracy, for diatomic molecular systems. This was mostly done through comparisons with standard-grid results. No effort was made to reproduce either accurate theoretical (within DFT as well as wave-function based) or experimental results, which may be considered in future, along with non-linear molecular systems. These would necessarily require the inclusion of more precise $F_{opt}$, especially for higher-order derivatives, more accurate XC functionals having correct short- as well as long-range properties.

V Future and outlook

We have presented a detailed study on the performance of CCG in the context of $\mu,\bm{\alpha},\bm{\beta}$ of diatomic molecules, using first-principles pseudopotential DFT formalism in amalgamation with the FF method. This was achieved through an accurate representation of FF formalism in real-space grid. The viability and feasibility of this approach has been demonstrated by applying it to a set of three diatomic molecules. Four different representations of XC functional was invoked. The effect of spatial grid as well as optimal electric field was analyzed in detail. It is quite gratifying that our results are in excellent agreement with those from standard program using atom-centered grid. In addition, for the first time, the effectiveness of pseudopotential LFK basis set (in CCG) is compared with the all-electron results, far from equilibrium. Application of this approach to more chemical systems would further enhance its success, which may be pursued in future. To conclude, pseudopotential-CCG can offer fairly accurate and reliable results for electric response properties of many-electron systems.

VI Acknowledgement

AG is grateful to UGC for a Senior Research Fellowship (SRF). TM acknowledges IISER Kolkata for a Junior Research fellowship (JRF). Financial support from DST SERB, New Delhi, India (sanction order number EMR/2014/000838) is gratefully acknowledged. We sincerely thank the two anonymous referees for their valuable and constructive suggestion, which have improved the manuscript greatly.