Static polarizability and hyperpolarizability in atoms and molecules   through a Cartesian-grid DFT

Tanmay Mandal; Abhisek Ghosal; Amlan K. Roy

arXiv:1904.11169·physics.chem-ph·April 26, 2019

Static polarizability and hyperpolarizability in atoms and molecules through a Cartesian-grid DFT

Tanmay Mandal, Abhisek Ghosal, Amlan K. Roy

PDF

Open Access

TL;DR

This paper presents a real-space Cartesian grid DFT approach to calculate static electric response properties of atoms and molecules, demonstrating high accuracy and feasibility across various properties and functionals.

Contribution

It introduces a systematic Cartesian-grid DFT method for static polarizability and hyperpolarizability, validated against literature and experimental data.

Findings

01

Excellent agreement with atom-centered-grid results

02

Effective use of non-uniform grid-optimization technique

03

Viable alternative for response property calculations

Abstract

Static electric response properties of atoms and molecules are reported within the real-space Cartesian grid implementation of pseudopotential Kohn-Sham (KS) density functional theory (DFT). A detailed systematic investigation is made for a representative set of atoms and molecules, through a number of properties like total ground-state electronic energies, permanent dipole moment ( $μ$ ), static average dipole polarizability ( $\overline{α}$ ). This is further extended to first-hyperpolarizability ( $β$ ) in molecules. It employs a recently developed non-uniform grid-optimization technique, with a suitably chosen fixed initial applied field. A simple variant of the finite-field method, using a rational function fit to the dipole moment with respect to electric field, is adopted. We make use of Labello-Ferreira-Kurtz (LFK) basis set, which has performed quite well in…

Tables4

Table 1. Table 1: Average static polarizability, α ¯ ¯ 𝛼 \overline{\alpha} for atoms in Set CCG-A, using FF KS method, for four XC functionals, along with the respective MSE and MAE’s. All results in a.u. See text for details.

	Average static polarizability ( $\bar{α}$ )
Atom	LDA	BLYP	PBE	LBVWN	Literature^a
Be( ${}^{1}S$ )	44.49	43.43	43.10	40.81	37.79
B( ${}^{2}P$ )	22.24	21.88	21.11	18.81	20.45
C( ${}^{3}P$ )	13.68	13.55	13.50	10.68	11.88
N( ${}^{4}S$ )	8.31	8.29	8.29	6.83	7.42
O( ${}^{3}S$ )	5.62	5.48	5.47	4.21	5.41
Mg( ${}^{1}S$ )	76.91	75.02	74.23	70.51	71.53
Al( ${}^{2}P$ )	48.26	48.37	47.03	40.50	45.89
Si( ${}^{3}P$ )	37.50	37.89	36.23	35.13	36.31
P( ${}^{4}S$ )	28.68	28.27	27.91	24.17	24.50
S( ${}^{3}P$ )	21.75	20.91	20.54	18.28	19.57
Cl( ${}^{2}P$ )	16.25	16.51	15.73	13.84	14.71
Ar( ${}^{1}S$ )	12.14	12.02	11.88	10.43	11.07
Ca( ${}^{1}S$ )	165.32	160.14	159.41	150.78	153.86
Kr( ${}^{1}S$ )	18.21	18.14	17.86	15.59	16.77
Xe( ${}^{1}S$ )	28.67	28.42	28.04	25.02	27.29
MSE	$-$ 2.91	$-$ 2.26	$-$ 1.74	1.25	–
MAE	2.91	2.26	1.74	1.65	–

Table 2. Table 2: Ground-state energy and permanent static dipole moment (absolute) for all the molecules in CCG-M set, for four XC functionals. All results are in a.u. See text for details.

Molecule	$- ⟨ E ⟩$				$𝝁$
	LDA	BLYP	PBE	LBVWN	LDA	BLYP	PBE	LBVWN	Expt.^†
F₂	47.91383	48.13149	48.12858	47.91426	–	–	–	–	–
Cl₂	29.65556	29.70811	29.40441	29.65489	–	–	–	–	–
Br₂	26.63157	26.66731	26.73636	26.63055	–	–	–	–	–
I₂	22.81646	22.82037	22.90226	22.81509	–	–	–	–	–
HF	24.64387	24.76034	24.76552	24.64197	0.70315	0.68988	0.69307	0.75623	0.71604
CO	21.49163	21.56884	21.59654	21.48968	0.09824	0.07539	0.09929	0.05770	0.04406
HCl	15.43342	15.47383	15.50513	15.43286	0.43825	0.42337	0.43420	0.45357	0.42490
HBr	13.90305	13.93459	13.96864	13.90218	0.31610	0.29881	0.31207	0.30760	0.32536
HI	11.98212	11.99865	12.03800	11.98063	0.18224	0.16100	0.17944	0.11947	0.17625
H₂O	17.06428	17.14538	17.16016	17.06321	0.71610	0.69956	0.70607	0.76583	0.7278
H₂S	11.25031	11.28181	11.31490	11.24899	0.44115	0.41985	0.43639	0.39553	0.38162
H₂Se	10.47641	10.50418	10.53809	10.47477	0.28488	0.25717	0.28047	0.19486	0.24668
NH₃	11.61575	11.67003	11.69468	11.61508	0.57940	0.57091	0.57498	0.59063	0.57834
PH₃	8.23604	8.26939	8.29887	8.23367	0.23137	0.19833	0.22442	0.11764	0.22818
P₄	26.10747	26.10747	26.20419	26.10558	–	–	–	–	–
PCl₃	51.15601	51.20379	51.34330	51.15280	0.17708	0.21547	0.18162	0.35481	0.30687
SiH₃Cl	20.49531	20.55085	20.60549	20.49152	0.50313	0.50656	0.49827	0.58014	0.51539
SiH₄	6.16406	6.21318	6.23146	6.16040	–	–	–	–	–
CH₃Cl	22.23814	22.28972	22.35264	22.23688	0.73076	0.73269	0.72914	0.71637	0.73571
CH₂Cl₂	36.51145	36.56834	36.66719	36.50986	0.63290	0.65394	0.63676	0.62398	0.62948
CH₃Br	20.70879	20.75325	20.81866	20.70720	0.71377	0.72486	0.71875	0.63353	0.71210
CHCl₃	50.78407	50.84088	50.97698	50.78216	0.42682	0.43710	0.42855	0.42664	0.39736
CH₄	7.96412	7.96412	8.03747	7.96306	–	–	–	–	–
CCl₄	65.05251	65.10947	65.28327	65.05034	–	–	–	–	–
C₂H₂	12.31925	12.35341	12.39988	12.31819	–	–	–	–	–
C₂H₄	13.55522	13.55522	13.55522	13.55386	–	–	–	–	–
Si₂H₆	11.19320	11.25222	11.29862	11.18691	–	–	–	–	–
C₃H₈	21.58015	21.64097	21.73504	21.57769	0.04065	0.03925	0.03844	0.03102	0.03304
C₄H ${}^{‡}_{6}$	25.98404	26.03255	26.14425	25.98181	–	–	–	–	–

Table 3. Table 3: Average static polarizability, α ¯ ¯ 𝛼 \overline{\alpha} for CCG-M Set, using FF KS method for different XC functionals, along with MSE and MAE. All results are given in a.u. See text for details.

Molecule	$\bar{α}$				Experiment^†
	LDA	BLYP	PBE	LBVWN	(zero frequency)
F₂	8.96	9.04	8.89	7.56	8.34
Cl₂	32.03	31.76	32.20	28.42	30.37
Br₂	46.60	45.89	45.60	41.65	43.7
I₂	72.68	71.54	71.07	66.10	70.36
CO	13.66	13.57	13.43	11.68	13.04
HF	6.24	6.25	6.15	4.88	5.09
HCl	18.79	18.55	18.37	16.07	17.40
HBr	25.75	25.19	25.09	22.22	23.78
HI	37.97	37.09	36.96	33.35	35.30
H₂O	10.52	10.41	10.26	8.91	9.52
H₂S	26.45	26.01	25.76	22.65	24.68
H₂Se	31.74	31.29	30.85	27.00	29.68
NH₃	15.43	15.31	15.08	12.59	14.61
PH₃	32.49	32.69	29.44	26.29	30.90
P₄	95.75	94.90	93.85	88.44	91.70
PCl₃	74.13	73.34	72.86	67.03	70.30
SiH₃Cl	44.93	43.81	43.86	39.93	35.8
SiH₄	34.07	32.80	33.01	29.91	31.94
CH₃Cl	31.95	31.72	31.28	27.50	30.00
CH₂Cl₂	47.48	46.94	46.55	42.46	44.89
CH₃Br	38.70	38.33	37.87	34.35	36.76
CHCl₃	60.12	59.52	58.98	52.02	56.22
CH₄	17.98	17.69	17.44	15.70	17.24
CCl₄	74.41	73.66	73.06	67.44	69.23
C₂H₂	24.46	24.78	24.02	20.86	22.67
C₂H₄	29.21	29.34	28.60	25.53	27.72
Si₂H₆	65.76	63.32	63.81	59.02	63.53
C₃H₈	43.94	43.23	42.69	39.27	42.12
C₄H ${}^{‡}_{6}$	59.34	59.64	58.27	53.07	54.64
MAE	2.59	1.95	1.61	2.33	–
MSE	$-$ 2.59	$-$ 1.93	$-$ 1.51	2.18	–

Table 4. Table 4: The components of first hyperpolarizability for some selected molecules from CCG-M set, for four XC functionals. All results are in a.u. See text for details.

Molecule	LDA	BLYP	PBE	LBVWN	LDA	BLYP	PBE	LBVWN	LDA	BLYP	PBE	LBVWN
	$β_{x x z}$				$β_{y y z}$				$β_{z z z}$
H₂S^a	$-$ 12.41	$-$ 14.30	$-$ 12.21	$-$ 4.93	6.96	5.80	6.07	8.78	$-$ 25.07	$-$ 27.54	$-$ 24.92	$-$ 7.31
H₂Se	$-$ 23.65	$-$ 32.56	$-$ 27.08	$-$ 8.95	52.99	46.07	48.75	40.85	9.24	0.58	6.58	15.67
PH₃^b	6.07	4.56	6.38	5.19	6.07	4.56	6.39	5.19	20.73	5.63	14.79	6.12
SiH₃Cl	23.01	18.19	19.86	14.93	23.01	18.18	19.86	14.93	58.69	53.77	55.75	46.66
CH₃Cl	$-$ 6.64	$-$ 8.82	$-$ 8.44	0.34	$-$ 6.65	$-$ 8.81	$-$ 8.43	0.33	21.24	18.87	17.19	19.55
CH₂Cl₂	14.45	14.62	14.66	2.13	$-$ 26.39	$-$ 26.68	$-$ 25.51	$-$ 19.88	23.99	24.12	24.48	31.15
CHCl₃^c	$-$ 19.42	$-$ 18.65	$-$ 17.82	$-$ 11.75	$-$ 18.49	$-$ 17.80	$-$ 16.92	$-$ 11.42	15.31	17.06	15.94	2.17
H₂O^d	18.73	14.76	16.67	8.60	16.19	16.26	15.85	11.40	31.78	26.75	28.75	19.75
NH₃^e	$-$ 13.97	$-$ 15.65	$-$ 14.15	$-$ 9.73	$-$ 13.97	$-$ 15.65	$-$ 14.15	$-$ 9.73	$-$ 55.43	$-$ 51.10	$-$ 51.20	$-$ 26.37
	$β_{x x y}$				$β_{y z z}$				$β_{y y y}$
C₃H₈	1.04	2.99	1.69	1.82	$-$ 25.01	$-$ 25.15	$-$ 23.92	$-$ 14.13	$-$ 28.59	$-$ 26.24	$-$ 26.30	$-$ 11.06
	$β_{x y y}$				$β_{x z z}$				$β_{x x x}$
CH₃Br	42.62	45.80	4.84	21.94	42.59	45.86	44.13	21.94	17.24	18.72	21.43	5.40
	$β_{x x z}$ = $β_{y y z}$				$β_{z z z}$
CO^f	8.55	7.88	8.42	3.44	31.83	30.29	29.80	21.48
HF^g	$-$ 3.59	$-$ 3.24	$-$ 3.22	$-$ 1.51	$-$ 14.09	$-$ 14.04	$-$ 13.52	$-$ 9.24
HCl^h	8.27	6.30	7.19	3.78	20.77	19.60	18.91	15.19
HBrⁱ	7.28	4.16	5.79	2.21	23.13	20.90	20.52	15.61
HI	$-$ 3.32	1.19	$-$ 1.80	2.39	$-$ 16.48	$-$ 12.64	$-$ 13.22	$-$ 9.22

Equations27

μ_{i} (F) = μ_{i} (0) + j = 1 \sum 3 α_{ij} F_{j} + \frac{1}{2} j, k = 1 \sum 3 β_{ij k} F_{j} F_{k} + \dots .

μ_{i} (F) = μ_{i} (0) + j = 1 \sum 3 α_{ij} F_{j} + \frac{1}{2} j, k = 1 \sum 3 β_{ij k} F_{j} F_{k} + \dots .

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

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

v_{i o n}^{p} (r) = R_{a} \sum v_{i o n, a}^{p} (r - R_{a}) .

v_{i o n}^{p} (r) = R_{a} \sum v_{i o n, a}^{p} (r - R_{a}) .

n \sum F_{mn}^{σ} C_{ni}^{σ} = ϵ_{i}^{σ} n \sum σ S_{mn} C_{ni}^{σ},

n \sum F_{mn}^{σ} C_{ni}^{σ} = ϵ_{i}^{σ} n \sum σ S_{mn} C_{ni}^{σ},

F_{mn}^{σ} = H_{mn}^{core} + J_{mn} + F_{mn}^{x c^{σ}} .

F_{mn}^{σ} = H_{mn}^{core} + J_{mn} + F_{mn}^{x c^{σ}} .

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},

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},

ρ (r_{g}) = m, n \sum P_{mn} χ_{m} (r_{g}) χ_{n} (r_{g}) .

ρ (r_{g}) = m, n \sum P_{mn} χ_{m} (r_{g}) χ_{n} (r_{g}) .

⟨ χ_{m} (r_{g}) ∣ v_{h x c} (r_{g}) ∣ χ_{n} (r_{g})⟩ = h_{x} h_{y} h_{z} g \sum χ_{m} (r_{g}) v_{h x c} (r_{g}) χ_{n} (r_{g}) .

⟨ χ_{m} (r_{g}) ∣ v_{h x c} (r_{g}) ∣ χ_{n} (r_{g})⟩ = h_{x} h_{y} h_{z} g \sum χ_{m} (r_{g}) v_{h x c} (r_{g}) χ_{n} (r_{g}) .

μ (F) = \frac{a + b F + c F ^{2} + d F ^{3} \dots}{1 + B F + C F ^{2} + D F ^{3} \dots} .

μ (F) = \frac{a + b F + c F ^{2} + d F ^{3} \dots}{1 + B F + C F ^{2} + D F ^{3} \dots} .

F_{n} = x^{n} F_{0}, x = 2^{\frac{p}{100}}

F_{n} = x^{n} F_{0}, x = 2^{\frac{p}{100}}

μ (F) = \frac{a + b F + c F ^{2} + d F ^{3}}{1 + B F + C F ^{2}},

μ (F) = \frac{a + b F + c F ^{2} + d F ^{3}}{1 + B F + C F ^{2}},

- F_{0} μ (F_{0}) F_{0} μ (- F_{0}) - x F_{0} μ (x F_{0}) x F_{0} μ (- x F_{0}) - x^{2} F_{0}^{2} μ (x^{2} F_{0}^{2}) x^{2} F_{0}^{2} μ (- x^{2} F_{0}^{2}) - F_{0}^{2} μ (F_{0}) - F_{0}^{2} μ (- F_{0}) - x^{2} F_{0}^{2} μ (x F_{0}) - x^{2} F_{0}^{2} μ (- x F_{0}) - x^{4} F_{0}^{4} μ (x^{2} F_{0}^{2}) - x^{4} F_{0}^{4} μ (- x^{2} F_{0}^{2}) F_{0} - F_{0} x F_{0} - x F_{0} x^{2} F_{0}^{2} - x^{2} F_{0}^{2} F_{0}^{2} F_{0}^{2} x^{2} F_{0}^{2} x^{2} F_{0}^{2} x^{4} F_{0}^{4} x^{4} F_{0}^{4} F_{0}^{3} - F_{0}^{3} x^{3} F_{0}^{3} - x^{3} F_{0}^{3} x^{6} F_{0}^{6} - x^{6} F_{0}^{6} B C b c d = μ (F_{0}) - μ (0) μ (- F_{0}) - μ (0) μ (x F_{0}) - μ (0) μ (- x F_{0}) - μ (0) μ (x^{2} F_{0}^{2}) - μ (0) μ (- x^{2} F_{0}^{2}) - μ (0)

- F_{0} μ (F_{0}) F_{0} μ (- F_{0}) - x F_{0} μ (x F_{0}) x F_{0} μ (- x F_{0}) - x^{2} F_{0}^{2} μ (x^{2} F_{0}^{2}) x^{2} F_{0}^{2} μ (- x^{2} F_{0}^{2}) - F_{0}^{2} μ (F_{0}) - F_{0}^{2} μ (- F_{0}) - x^{2} F_{0}^{2} μ (x F_{0}) - x^{2} F_{0}^{2} μ (- x F_{0}) - x^{4} F_{0}^{4} μ (x^{2} F_{0}^{2}) - x^{4} F_{0}^{4} μ (- x^{2} F_{0}^{2}) F_{0} - F_{0} x F_{0} - x F_{0} x^{2} F_{0}^{2} - x^{2} F_{0}^{2} F_{0}^{2} F_{0}^{2} x^{2} F_{0}^{2} x^{2} F_{0}^{2} x^{4} F_{0}^{4} x^{4} F_{0}^{4} F_{0}^{3} - F_{0}^{3} x^{3} F_{0}^{3} - x^{3} F_{0}^{3} x^{6} F_{0}^{6} - x^{6} F_{0}^{6} B C b c d = μ (F_{0}) - μ (0) μ (- F_{0}) - μ (0) μ (x F_{0}) - μ (0) μ (- x F_{0}) - μ (0) μ (x^{2} F_{0}^{2}) - μ (0) μ (- x^{2} F_{0}^{2}) - μ (0)

α = μ^{'} (0) = b - a B,

α = μ^{'} (0) = b - a B,

β = μ^{''} (0) = 2 c - 2 b B - 2 a C - 2 a B^{2} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsAdvanced Chemical Physics Studies · Nonlinear Optical Materials Research · Spectroscopy and Quantum Chemical Studies

Full text

Static polarizability and hyperpolarizability in atoms and molecules through a Cartesian-grid DFT

Tanmay Mandal, Abhisek Ghosal

Amlan K. Roy

Email: [email protected], [email protected].

Department of Chemical Sciences

Indian Institute of Science Education and Research (IISER) Kolkata, Mohanpur-741246, Nadia, West Bengal, India

Abstract

Static electric response properties of atoms and molecules are reported within the real-space Cartesian grid implementation of pseudopotential Kohn-Sham (KS) density functional theory (DFT). A detailed systematic investigation is made for a representative set of atoms and molecules, through a number of properties like total ground-state electronic energies, permanent dipole moment ( $\bm{\mu}$ ), static average dipole polarizability ( $\overline{\alpha}$ ). This is further extended to first-hyperpolarizability ( $\bm{\beta}$ ) in molecules. It employs a recently developed non-uniform grid-optimization technique, with a suitably chosen fixed initial applied field. A simple variant of the finite-field method, using a rational function fit to the dipole moment with respect to electric field, is adopted. We make use of Labello-Ferreira-Kurtz (LFK) basis set, which has performed quite well in these scenarios. To assess the efficacy and feasibility, four XC functionals such as LDA, BLYP, PBE and LBVWN are chosen. Present results are compared with available literature (both theoretical and experimental) values, whenever possible. In all instances, they show excellent agreement with the respective atom-centered-grid results, very widely used in many quantum chemistry programs. This demonstrates a viable alternative towards accurate prediction of response properties of many-electron systems in Cartesian coordinate grid.

Keywords: Density functional theory, polarizability, hyperpolarizability, exchange-correlation functional, Cartesian grid, basis set.

I Introduction

In past several decades, we have witnessed remarkable stride in ab initio quantum mechanical methods for accurate computation of various properties of many-electron systems, ranging from atoms, molecules to complex extended systems. To a large extent, this has been possible due to rapid advances in mathematical formalism, smart numerical algorithm, coupled with the developments in computer architecture. Amongst different theoretical approaches, Kohn-Sham density functional theory (KS-DFT) (hohenberg64, ; kohn65, ) has proved to be a highly effective, versatile tool for structural and electronic characterization of such systems. Within this framework, a considerable number of elegant and computationally efficient schemes are available to tackle these complex systems. They cover a broad range of approximations regarding the exchange-correlation (XC) functionals, keeping in mind the chemical accuracy as far as practicable. Due to its favorable computational cost, KS-DFT calculations are now routinely used in physical, chemical, biological and materials sciences, making them the preferred workhorse of electronic structure calculations for “real-world” materials (guliamov07, ). The progress in the subject has been reviewed in many elegant references (becke14, ; jones15, ).

The electric-dipole polarizability i.e., response of electric dipole to an external electric field, is a very important property of an atom, molecule or cluster (maroulis06, ; champagne10, ). Its accurate, efficient description plays an important role in understanding electro-magnetic field-matter interaction and inter-particle collision phenomena such as, second-order Stark effect, Rayleigh and Raman scattering, electron-detachment process (ghazaly05, ) etc. On the other hand, synthesis and design of novel non-linear optical materials and molecular assemblies having large non-linear optical coefficients, are particularly interesting from a technological point of view (kummel06, ). Theoretical advances and dramatic enhancements in computational power in recent years have made it possible to provide guiding principles in this direction.

A number of different theoretical approaches have been proposed in literature to compute electric response properties in linear and non-linear regime within the KS-DFT rubric. Some of the notable ones are: coupled-perturbed Kohn-Sham (CPKS) method (fournier90, ; colwell93, ), linear-response time dependent DFT (jansik05, ; helgaker12, ), perturbative sum-over states expression over all dipole-allowed electronic transitions (orr71, ; bishop94, ), numerical method using Sternheimer approach (talman12, ), auxiliary DFT (moreno08, ; espindola12, ), non-iterative CPKS (sophy08, ) as well as the fully numerical finite field (FF) method (bishop87, ; maroulis88, ). From a computational viewpoint, FF method is the least expensive as there is no requirement of analytical derivatives or excited-state information, and also easy to implement as only one-body Hamiltonian is perturbed by the applied field (kurtz90, ). That is why FF method is usually considered at a first glance to justify the suitability and viability of a newly proposed quantum chemistry method for computation of electric response (bulat05, ; wouters12, ; wergifosse14, ; wouters16, ). The latter in a system is governed by its static dipole polarizability $\left(\alpha_{ij}=\frac{\partial{\mu_{i}}}{\partial F_{j}}\right)$ as well as first $\left(\beta_{ijk}=\frac{\partial^{2}\mu_{i}}{\partial F_{j}\partial F_{k}}\right)$ and second $\left(\gamma_{ijkl}=\frac{\partial^{3}\mu_{i}}{\partial F_{j}\partial F_{k}\partial F_{l}}\right)$ hyperpolarizability, where $i,j$ denote the Cartesian directions $x,y,z$ . Pertinent coefficients are then obtained from a Taylor expansion of the dipole moment,

[TABLE]

where $\mu_{i}(\mathrm{\mathbf{F}})=\sum_{mn}P_{mn}(\mathrm{\mathbf{F}})\langle m|r_{i}|n\rangle+\sum_{\mathrm{\mathbf{R}}_{nuclei}}Z_{\mathrm{\mathbf{R}}}R_{i}$ , denotes molecular dipole moment and $P_{mn}(\mathrm{\mathbf{F}})$ the $(m,n)$ th element of perturbed density matrix $\mathrm{\mathbf{P}}$ in presence of $-\mathrm{\mathbf{r}}\cdot\mathrm{\mathbf{F}}$ perturbation, in dipole approximation. Alternatively one may also represent it in terms of field-induced energy, and both definitions are equivalent according to Hellmann-Feynman theorem. Now a vital factor in FF method is its numerical accuracy i.e, the system’s dipole moment is calculated in presence of applied electric field, and the appropriate derivatives are approximated by finite differences. Hence the strength of applied field should be chosen carefully such that (i) it is sufficiently large to overcome finite-precision artifacts for a meaningful evaluation of necessary finite differences, especially in the non-linear regime for hyperpolarizabilities (ii) it must also be small enough so that the contributions from remaining higher-order derivatives become negligible for one particular coefficient. It is well documented that this effect is quite pronounced, in particular for higher-order derivatives (bishop85, ), and to optimize it, some novel techniques, viz., polynomial or rational-function-based FF methods (with or without extrapolation), were proposed in the literature (mohammed17, ; patel17, ).

A fundamental problem in a given DFT calculation lies in the correct approximation of hitherto unknown and difficult XC functional, regarding its proper asymptotic behavior at long-range as long as spurious self-interaction correction (perdew81, ) which produces considerable error in these properties. In general, approximate local or semi-local functionals pose rather serious limitations regarding these issues. On the other hand, hybrid functionals (such as B3LYP or CAM-B3LYP) and their asymptotic-corrected variants were reported quite successfully in literature. Usually the obtained results are found to be as accurate as those from correlated wave-function based approaches (castet12a, ). This work employs following four XC functionals, namely, (i) local density approximation (LDA) using Vosko-Wilk-Nusair (VWN) (vosko80, ) correlation, (ii) BLYP (Becke becke88 exchange with Lee-Yang-Parr (LYP) lee88 correlation) (iii) Perdew-Burke-Ernzerhof (PBE) (perdew96, ) (iv) asymptotically corrected van Leeuwen-Baerends (LB94) (leeuwen94, ) exchange in conjunction with VWN correlation.

Recently KS-DFT was implemented in Cartesian coordinate grid (CCG) roy08 ; roy08a ; roy11 ; ghosal16 , offering accurate results on atoms, molecules within an LCAO-MO ansatz. For a decent number of systems, several electronic properties including energy, $\mu$ , as well as its derivatives $\alpha,\beta$ were correctly presented. The electric responses were provided for three diatomic hydrides (HCl, HBr, HI) by including a coupling term (interaction of applied electric field and matter) in the respective one-body KS Hamiltonian (ghosal18, ). Through an accurate representation of FF formalism in real-space grid, corresponding results were compared and contrasted with the widely used, popular atom-centered grid ones (two were virtually indistinguishable) within a pseudopotential approximation. A detailed analysis of grid convergence and its impact on response was considered by taking four XC functionals, In order to broaden its scope of applicability, here we extend this to a larger set of atoms, molecules, which is the primary motivation of this communication. A secondary objective is to engage and ascertain the suitability of a recently proposed rational function approximation patel17 for their estimation. Note that in previous work, this was done through a Taylor FF approach.

It is well acknowledged that the success of such calculations crucially depends on choice of an appropriate good-quality basis set. In other words, it should accurately describe the response of electrons to an external perturbation; an essential requisite being that it preferably includes polarization and diffused functions miadokova97 . A good candidate that satisfies these criteria is the Labello-Ferreira-Kurtz (LFK) basis labello05 utilized in our previous work (ghosal18, ), which has been found to be quite successful in both non-relativistic and relativistic domain. Current results are compared with available references (both theoretical and experimental), whenever possible, to gauge its performance. The article is organized as follows. Section II gives principal aspects of the method along with required computational details. Results are discussed in Sec. III, while the future and prospect is commented in Sec. IV.

II Methodology and computational aspects

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 ionic pseudopotential, expressed as in the following,

[TABLE]

Here $v^{p}_{ion,a}$ represents ion-core pseudopotential associated with atom A, located at $\mathrm{\mathbf{R}}_{a}$ . In presence of applied electric field, $v_{ext}(\mathrm{\mathbf{r}})$ contains contributions from both nuclear and field effects. The classical Coulomb (Hartree) potential, $v_{h}[\rho(\mathrm{\mathbf{r}})]$ describes electrostatic interaction amongst valence electrons whereas $v_{xc}[\rho(\mathrm{\mathbf{r}})]$ signifies the non-classical XC potential, and $\{\psi^{\sigma}_{i},\sigma=\uparrow\mathrm{or}\downarrow\mathrm{spin}\}$ corresponds to a set of $N$ occupied orthonormal spin molecular orbitals (MO). In LCAO-MO approximation, the coefficients for expansion of 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 respective MO coefficients $\{C_{ni}^{\sigma}\}$ for a given MO, $\psi_{i}^{\sigma}(\mathrm{\mathbf{r}})$ , whereas $\mathrm{\mathbf{S}}$ signifies the usual overlap matrix corresponding to elements $S_{mn}$ and $\bm{\epsilon}^{\sigma}$ refers to diagonal matrix of respective MO eigenvalues, $\{\epsilon_{i}^{\sigma}\}$ . The KS matrix has elements $F_{mn}^{\sigma}$ partitioned as,

[TABLE]

In this equation, $H_{mn}^{\mathrm{core}}$ contains all one-electron contributions and in presence of applied electric field this could be further separated as, $H_{mn}^{\mathrm{core}}=H_{mn}^{(0)}+F_{i}\langle m|r_{i}|n\rangle$ , where $H_{mn}^{(0)}$ includes kinetic energy, nuclear-electron attraction and pseudopotential matrix elements, whereas the latter term signifies external electric field perturbation within dipole approximation. Furthermore, $F_{i}$ ( $i=x,y,z$ ) corresponds to ith component of field $\mathrm{\mathbf{F}}$ , while $\langle m|r_{i}|n\rangle$ provides the dipole moment integral corresponding to length vector $\mathrm{\mathbf{r}}$ . Last two terms $J_{mn},F_{mn}^{xc\sigma}$ refer to usual contributions from classical Hartree and XC potentials respectively.

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 a correct, reliable estimation of two-electron KS matrix elements. In this work, they are calculated directly through straightforward numerical integration on grid as ( $v_{hxc}$ refers to Coulomb and XC potential combined),

[TABLE]

The detailed construction of various potentials in grid has been well documented (roy08, ; roy08a, ; roy11, ; ghosal16, ).

As already mentioned, a critical problem of finite difference method for optical property calculation arises from the delicate nature of field–they satisfy a rather narrow range of field strengths. In order to alleviate this problem, we follow a recently published procedure (patel17, ), where the energy expression is fitted in terms of a rational function with respect to electric-field coefficients. An analogous fitting strategy for the induced dipole moment in terms of electric field is in invoked here, i.e.,

[TABLE]

where a,b,c,d, $\cdots$ and B,C,D, $\cdots$ are fitting coefficients. This is a generalized form of the Taylor expansion, as by setting all the denominator coefficients to zero, the latter expansion is recovered. This recipe provides a less sensitive (and hence consequently more effective) dependence on the electric field, as it increases latter’s range. In FF technique, dipole moment of a given species needs to be computed at different field strengths, and its proper selection begins with a convenient choice of initial field strength ( $F_{0}$ ), around which the field is distributed. This is achieved here via the following proposal patel17 ,

[TABLE]

Here we choose their recommended value of $p$ (namely 50), corresponding to a geometry progression–this was reached on the basis of an elaborate analysis of $\bm{\alpha}$ and $\bm{\gamma}$ for a set of 121 and $\bm{\beta}$ for 91 molecules. It is well established that the optimal value of $F_{0}$ plays a decisive role to ensure accuracy in these non-linear properties. On the basis of our recent work ghosal18 , where this aspect was discussed at length for $\bm{\beta}$ , it is found to be appropriate to select an initial value of $10^{-2.5}$ , which we exploit here. Next the optimal (numerator and denominator degrees) form of rational function are taken from (patel17, ) as,

[TABLE]

corresponding to four and three terms in the numerator and denominator respectively. Now putting the value of $\bm{\mu}$ at $\mathrm{\mathbf{F}}=0$ in Eq. (11) gives rise to the trivial relation $\bm{\mu}(0)=a$ . The remaining unknown coefficients are determined using different $F_{n}$ values in this equation. For five unknown coefficients the minimum equations that can be constructed is six, because both $+F_{n}$ and $-F_{n}$ are used, with each $F_{n}$ producing two equations. By putting $\bm{\mu}(0)$ for $a$ , the substituted equations may be recast in the matrix form $A\mathrm{\mathbf{x}}=\mathrm{\mathbf{b}}$ as,

[TABLE]

Because both ( $+$ )ve and ( $-$ )ve fields are used, solution of this matrix equation is overdetermined, this can be solved by means of either least-squares method or by disregarding one of the equations. We adapt the second procedure, where the arbitrary elimination of any one of six equations was found to be valid for all the systems concerned. Finally, our desired properties can be determined by taking appropriate derivatives of Eq. (11) at F=0 as,

[TABLE]

A few practical aspects regarding the current compilation is in order. All calculations are performed using norm-conserving pseudopotential at the experimental geometries taken from NIST database (johnson16, ). The needed one-electron (excepting the pseudopotential) integrals were generated by standard recursion relations (obara86, ) using Cartesian Gaussian-type orbitals as primitive basis functions, whereas the latter matrix elements in Gaussian orbitals are imported from GAMESS (schmidt93, ) package. The influence of spatial grid on energy components has been detailed in a previous article (ghosal18, ) with respect to sparsity of grid (governed by $N_{x},N_{y},N_{z}$ ) and grid spacing (determined by $h_{r}$ ); and hence not repeated here. A simple grid optimization strategy has been maintained, which guarantees a grid accuracy of at least $5\times 10^{-6}$ a.u., throughout, at a fixed $h_{r}=0.3$ . It was observed that the optimal non-uniform grid marginally varies from functional to functional. As done earlier, here also we employ the so-called Labello-Ferreira-Kurtz (LFK) basis set advocated in (labello05, ) based on the procedure that incorporates diffuse and polarization functions in the well-known Sadlej (sadlej92, ) basis–these are adopted from EMSL library (feller96, ). For efficient calculation of Hartree potential in real space, the standard discrete Fourier transform package FFTW3 (fftw05, ) has been used.

The properties are explored for four representative XC functionals, viz., LDA, BLYP, PBE, LBVWN. The middle two functionals were imported from density functional repository program (repository, ). The self-consistent convergence criteria imposed in this communication is slightly tighter than our earlier work (roy08, ; roy08a, ; roy11, ; ghosal16, ); this is to generate a more accurate perturbed density matrix. Convergence of following quantities was followed, viz., (i) orbital energy difference between two successive iterations (ii) absolute deviation in a density matrix element. They both were required to remain below a certain prescribed threshold set to $1\times 10^{-8}$ a.u.; this ensured that total energy maintained a convergence of at least this much. To accelerate the calculations, unperturbed (field-free) density matrix was used as trial input. The resulting generalized matrix-eigenvalue problem as well as Eq. (12) was solved through standard LAPACK (anderson99, ) routines. Scaling properties have been discussed elsewhere roy08 . The present calculations have been performed using an in-house computer program originally initiated by Roy in the references roy08 ; roy08a ; roy11 , and subsequently extended ghosal16 ; ghosal18 by the two co-authors of this communication.

III Results and Discussion

We begin this section by categorizing two sets: first one for atoms containing 15 of them, labeled CCG-A and the second set having 29 selected molecules, represented by CCG-M. Both sets contain open and close-shell electronic systems, with maximum number of total valence electrons being 32. The self-consistent field convergence was thoroughly checked with respect to all relevant parameters (like grid and field optimization) as delineated in previous section and elaborated in an earlier publication ghosal18 , both in absence/presence of the electric field. Since the converged properties thus obtained reproduce standard GAMESS results very well for all XC functionals available therein, these reference values are omitted henceforth. This matching agreement is expected and in keeping with the findings of ghosal18 , where this was demonstrated for three diatomic hydrides. A systematic investigation is carried out for these two sets through the following quantities, namely, non-relativistic ground-state total energy, $\bm{\mu}$ , non-zero components of $\bm{\alpha}$ , $\bm{\beta}$ , along with the average or mean polarizability, $\overline{\alpha}$ .

At first, the LDA, BLYP, PBE and LBVWN results for $\overline{\alpha}$ are provided for Set CCG-A in Table I. As atoms contain inversion of symmetry, $\bm{\beta}$ vanishes, and the lowest non-vanishing non-linear coefficient is $\gamma$ (which is not considered in this work). For the sake of comparison, available theoretical values from NIST database (johnson16, ) are quoted in last column. It reveals an interesting fact that, all three traditional functionals (LDA, BLYP, PBE) overestimate the references in the ranges of $3$ - $11\%$ , $1$ - $9\%$ and $3$ - $8\%$ respectively; however LBVWN offers underestimation in all cases (with exception for Be) by 5-9%. Moreover we also provide the respective mean signed error (MSE) and mean absolute error (MAE) with respect to experimental result, for all the XC functionals, at the bottom of this table. It is observed that the first three functionals systematically overestimate the experimental values, while LBVWN shows underestimation except the lone case of Be. Furthermore the errors also diminish as one goes from columns 2-4, corresponding to LDA, BLYP and PBE, having same magnitudes for MSE and MAE. For LBVWN, the errors are lowest amongst all the four functionals; however these two deviations differ in magnitude considerably. Similar conclusions have also been drawn regarding the pattern behavior of these functionals for $\overline{\alpha}$ in (vasiliev10, ), where it was conjectured that, a significant improvement may be accomplished for $\bm{\alpha}$ by combining LB94 potential (in asymptotic region) with LDA (corrected for derivative discontinuity in the bulk region) exchange suitably. This leads to a more proper representation of exchange which approaches the experimental results quite closely (casida98, ) at a lower level of computational cost than standard XC functionals.

Next in Table II, we present ground-state energy and $\bm{\mu}$ (absolute value) for all 29 molecules in Set CCG-M. This covers diatomic to hexa-atomic systems containing both close and open shells–the equilibrium experimental geometries are taken from NIST computational chemistry database johnson16 . Note that these values pertain to only the electronic part; neither geometry relaxation in presence of electric field nor vibrational contribution is considered. In all occasions, total energies using BLYP and PBE XC functional are observed to be consistently lower than LDA and LBVWN. Keeping in mind the performance of results for various XC functionals in our earlier work roy08 ; roy08a ; roy11 ; ghosal16 ; ghosal18 , one can safely conclude that current total electronic energies for this data set are equally accurate and trustworthy. Several theoretical and experimental results are available in literature. In order to compare, we report only its non-zero components along with some selected experimental results. The computed value of zero components of $\bm{\mu}$ for non-polar molecules have been correctly reproduced in this calculation; and henceforth not mentioned. In polar molecules also they show reasonably good agreement with experimental results. Leaving aside the two lone cases (CO plus PCl3 for LDA, BLYP, PBE, and PH3 for LBVWN), the overall maximum absolute deviation hovers around $13\%$ (for LDA, PBE) and $10\%$ (for BLYP, LBVWN) respectively from experimental counterparts.

Next we move towards the average polarizability, $\overline{\alpha}$ of molecules in Set CCG-M. Table III reports these at equilibrium geometries of NIST computational chemistry database johnson16 , for same four XC functionals of previous tables. To put things in perspective, we also quote corresponding experimental values of zero frequency (containing only electronic part) in column six, from the compilation of (hohm13, ). It uncovers that, all three traditional functionals (LDA, BLYP, PBE) overestimate experimental references in the range of $3$ - $11\%$ , $1$ - $9\%$ and $3$ - $8\%$ respectively; however, LBVWN separates out from them by underestimating within $5$ - $9\%$ (with exception of Be). We have also performed the MAE and MSE analysis, as given in the bottom. It follows similar kind of behavior regarding first three functionals (LDA, BLYP and PBE), in harmony with those of CCG-A set in Table I. However, it is observed that the results using LBVWN shows much more error compared to BLYP and PBE in contrast to Table I.

In the last part, Table IV reports the non-zero components of $\bm{\beta}$ (using the $T$ convention) for some representative molecules (CO, HCl, HBr, HI, H2S, H2Se, PH3, SiH3Cl, CH3Cl, CH2Cl2, CHCl3, C3H8, CH3Br, H2O, NH3, HF) of CCG-M set (for same four XC functionals) to extend the scope of applicability. The results for HCl, HI and HBr using celebrated FF method are already mentioned in our previous article (ghosal16, ) and these remain nearly same in the present case too. Same field strength, as mentioned in Sec. II has been employed irrespective of the molecular system and XC functional. It is clear that for a particular molecule the components differ significantly from functional to functional–sometimes even the signs vary. One such candidate is HI, where $\beta_{xxz},\beta_{yyz}$ signs for LDA, PBE functionals are opposite from those of BLYP, PBE. In literature these properties have been pursued through a wide range experiments and theoretical calculations. As a matter of comparison, some selected high-level all-electron calculations (such as CCSD, CAS, CCSD(T)) in large basis sets (Sadlej, taug-cc-pVTZ, qaug-sadlej, NLO-II) are quoted, along with some selected experiments. Of course, as expected, for obvious reasons present results differ from these extended and elaborate results rather considerably. However this is aside the main focus of this work; here we wanted to extend the domain of CCG in the context of response properties. As evident from the foregoing discussion (along with the previous tables), this scheme very satisfactorily fulfills that objective, for a decent number of atoms, molecules. However, as well known, for better comparison with experiment and other sophisticated theoretical calculations, more suitable basis set as well as XC functionals having correct asymptotic potential at long range, would be necessary, which may be pursued in our future works.

IV Future and outlook

We have analyzed the performance of CCG in the context of $\bm{\mu},\bm{\alpha}$ for a set of 44 species (sets CCG-A and CCG-M) using first-principles pseudopotential DFT formalism combined with an optimized FF procedure. Additionally, $\bm{\beta}$ values were provided for ten molecules. This was realized through a recently prescribed rational function approximation of dipole moment in real-space grid within the LFK basis set. The suitability and viability of this simple yet effective scheme has been demonstrated for four different XC functionals. Comparison with existing theoretical and experimental results have been made, wherever possible. In all occasions, our present results show excellent agreement with the corresponding values from familiar atom-centered grid. For further accuracy, more precise basis set as well as XC functionals need to be invoked. In some DFT works yanai04 ; limacher09 , significant improvements in these quantities have been reported through the use of popular CAMB3LYP functional; we would like to consider this in a forthcoming communication. Moreover, it will also be interesting to study the dynamical hyper-polarizability within the rubric of time-dependent DFT. To conclude, pseudopotential-CCG can offer fairly accurate and reliable results for electric response in atoms and molecules.

V Acknowledgement

AG thanks UGC for a Senior Research Fellowship. TM very much appreciates a Junior Research fellowship from IISER Kolkata. Financial support from DST SERB, New Delhi, India (sanction order number EMR/2014/000838) is sincerely acknowledged. Constructive comments from anonymous referee have helped in improving the quality of this manuscript.

Bibliography69

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) P. Hohenberg and W. Kohn. Phys. Rev. , 1964, 136 , B 864.
2(2) W. Kohn and L. J. Sham. Phys. Rev. , 1965, 140 , A 1133.
3(3) O. Guliamov and L. Kronik. J. Phys. Chem. A , 2007, 111 , 2028.
4(4) A. D. Becke. J. Chem. Phys. , 2014, 140 , 18A 301.
5(5) R. O. Jones. Rev. Mod. Phys. , 2015, 87 , 897.
6(6) G. Maroulis (Ed.), Atoms, Molecules and Clusters in Electric Fields: Theoretical Approaches to the Calculation of Electric Polarizability; Imperial College Press: London, 2006.
7(7) B. Champagne. Specialist Periodical Reports: Chemical Modelling, Applications and Theory; M. Springborg (Ed.) Vols. 6, 7 . Royal Society of Chemistry, London, 2009, 2010.
8(8) M. O. A. EI Ghazaly, A. Svendsen, H. Bluhme, S. B. Nielsen and L. H. Andersen. Chem. Phys.Lett. , 2005, 405 , 278.