Density-functional theory study of ionic inhomogeneity in metal clusters using SC-ISJM

TL;DR

This study applies a self-consistent inhomogeneous jellium model to predict electronic and geometric properties of various metal clusters, validating results with first-principles calculations and capturing size and structural changes.

Contribution

The paper introduces a generalized inhomogeneous stabilized jellium model with two regions, improving predictions of cluster size and inter-layer distances compared to previous models.

Findings

Model accurately predicts size reductions in clusters.

Model captures size increases in some clusters.

Results agree well with first-principles calculations.

Abstract

In this work we have applied the recently formulated self-compressed inhomogeneous stabilized jellium model [M. Payami and T. Mahmoodi, Can. J. Phys. {\bf 89}, 967 (2011)] to describe the equilibrium electronic and geometric properties of atomic-closed-shell simple metal clusters of Al$_N$ ($N$=13, 19, 43, 55, 79, 87, 135, 141), Na$_N$, and Cs$_N$ ($N$=9, 15, 27, 51, 59, 65, 89, 113). To validate the results, we have also performed first-principles pseudo-potential calculations and used them as our reference. In the model, we have considered two regions consisting of "surface" and "inner" ones, the border separating them is sharp. This generalization makes possible to decouple the relaxations of different parts of the system. The results show that the present model correctly predicts the size reductions seen in most of the clusters. It also predicts increase in size of some clusters, as observed from first-principles results. The other property is the change of inter-layer distances, which in most cases were realized as contractions and in a few cases as expansions, are reproduced in good agreement with the atomic simulation results. For a better description of the properties, it is possible to improve the method of choosing the thicknesses or generalize the model to include more regions than just two.

Tables

Table 1. Table 1: Ground state energies of Al N , Na N , and Cs N clusters resulted from SC-ISJM. Here, Δ ​ r s † = r s † − r s B Δ superscript subscript 𝑟 𝑠 † superscript subscript 𝑟 𝑠 † superscript subscript 𝑟 𝑠 𝐵 \Delta r_{s}^{\dagger}=r_{s}^{\dagger}-r_{s}^{B} . For comparison, the energies of SC-SJM and SJM are also listed. As is seen, in all cases the inequality E SC − ISJM < E SC − SJM < E SJM subscript 𝐸 SC ISJM subscript 𝐸 SC SJM subscript 𝐸 SJM E_{\rm SC-ISJM}<E_{\rm SC-SJM}<E_{\rm SJM} holds.

		SC-ISJM				SC-SJM		SJM
atom	$N$	${\nu }^{†}$	$\Delta r_s{1}^{†}$	$\Delta r_s{2}^{†}$	$-E_{\mathrm{SC}-\mathrm{ISJM}}$	$\Delta r_s^{†}$	-$E_{\mathrm{SC}-\mathrm{SJM}}$	-$E_{\mathrm{SJM}}$
Al	13	1	+0.205	-0.209	27.1684	-0.143	27.1106	26.9684
	19	1	+0.185	-0.101	39.7020	-0.105	39.5816	39.4659
	43	1	+0.206	-0.169	90.0310	-0.079	89.6954	89.5462
	55	1	+0.206	-0.068	114.7170	-0.074	114.6332	114.4648
	79	1	+0.033	-0.163	165.4508	-0.063	164.9392	164.7558
	87	1	+0.014	-0.147	182.2520	-0.064	182.0282	181.8216
	135	1	+0.038	-0.052	282.6374	-0.051	282.3383	282.1326
	141	1	+0.030	-0.052	295.3159	-0.052	295.0944	294.8699
Na	9	1	+0.315	-0.202	1.9610	-0.205	1.9578	1.9524
	15	1	+0.394	-0.288	3.2851	-0.137	3.2705	3.2664
	27	1	-0.067	-0.096	5.9634	-0.114	5.9532	5.9478
	51	1	-0.078	-0.081	11.3871	-0.100	11.3735	11.3652
	59	1	+0.009	-0.268	13.2415	-0.103	13.2301	13.2197
	65	2	-0.174	-0.073	14.5532	-0.105	14.5438	14.5320
	89	1	-0.016	-0.083	20.0237	-0.088	20.0146	20.0033
	113	1	-0.042	-0.064	25.3881	-0.076	25.3646	25.3538
Cs	9	2	+0.500	-0.269	1.4594	-0.206	1.4577	1.4556
	15	2	+0.546	-0.242	2.4430	-0.143	2.4398	2.4381
	27	2	-0.051	-0.182	4.4394	-0.121	4.4341	4.4318
	51	1	-0.116	-0.067	8.4618	-0.106	8.4548	8.4512
	59	2	-0.054	-0.215	9.8225	-0.109	9.8225	9.8181
	65	1	-0.161	-0.059	10.8073	-0.105	10.8011	10.7966
	89	2	-0.035	-0.066	14.8623	-0.093	14.8590	14.8541
	113	2	-0.013	-0.143	18.8588	-0.078	18.8460	18.8415

Equations

$$E_{\mathrm{ISJM}}\left[n;[\{n_{+}^{(\alpha )}\}]\right]$$

$$E_{\mathrm{JM}}\left[n;[\{n_{+}^{(\alpha )}\}]\right]=T_s[n]+E_{xc}[n]+\frac{1}{2}\int d\mathbf{r}\varphi ([n,n_{+}];\mathbf{r})[n(\mathbf{r})-n_{+}(\mathbf{r})]$$

$$\varphi ([n,n_{+}];\mathbf{r})=\int d{\mathbf{r}}^{\prime }\frac{[n({\mathbf{r}}^{\prime })-n_{+}({\mathbf{r}}^{\prime })]}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}$$

$$n_{+}(\mathbf{r})=\sum _{\alpha }{n}_{+}^{(\alpha )}(\mathbf{r});n_{+}^{(\alpha )}(\mathbf{r})\equiv {\bar{n}}_{\alpha }{\Theta }^{(\alpha )}(\mathbf{r}),$$

$$⟨\delta v{⟩}_{\mathrm{WS}}({\bar{n}}_{\alpha })=\frac{1}{{\Omega }^{(\alpha )}}{\int }_{{\Omega }^{(\alpha )}}d\mathbf{r}\delta v(\mathbf{r});{\Omega }^{(\alpha )}=\frac{z}{{\bar{n}}_{\alpha }},$$

$$E_0\equiv E\left[n_0;[\{n_{+}^{(\alpha )†}\}]\right]=\underset{\{n_{+}^{(\alpha )}\}}{\min }E\left[n_0;[\{n_{+}^{(\alpha )}\}]\right]=\underset{\{n_{+}^{(\alpha )}\}}{\min }\left\{\underset{n}{\min }E\left[n;[\{n_{+}^{(\alpha )}\}]\right]\right\},$$

$$E_0$$

$${\bar{n}}^B{V}_{\alpha }^B={\bar{n}}_{\alpha }{V}_{\alpha }={\bar{n}}_{\alpha }^{†}{V}_{\alpha }^{†};\alpha =1,2$$

$${\left(R_2^B\right)}^3=\frac{3zN}{4\pi {\bar{n}}^B},$$

$$R_1^B=R_2^B-t,$$

$${\bar{n}}^B{\left(R_1^B\right)}^3={\bar{n}}_1{\left(R_1\right)}^3,$$

$${\bar{n}}^B\left[{\left(R_2^B\right)}^3-{\left(R_1^B\right)}^3\right]={\bar{n}}_2\left[{\left(R_2\right)}^3-{\left(R_1\right)}^3\right].$$

$$V(r)=\left\{\begin{array}[]{lll}-2\pi\left[\bar{n}_{1}(R_{1}^{2}-\frac{r^{2}}{3})+\bar{n}_{2}(R_{2}^{2}-R_{1}^{2})\right]&,&r\leq R_{1}\\ -\left\{\frac{Nz}{R_{2}}+\frac{4\pi}{3}\left[R_{1}^{3}(\bar{n}_{1}-\bar{n}_{2})\left(\frac{1}{r}-\frac{1}{R_{2}}\right)-\frac{\bar{n}_{2}}{2}\left(r^{2}-R_{2}^{2}\right)\right]\right\}&,&R_{1}\leq r\leq R_{2}\\ -\frac{Nz}{r}&,&r>R_{2}.\end{array}\right.$$

$${\Psi }_i(\mathbf{r})=\psi (r)Y_{l,m}(\theta ,\varphi ),$$

$${\hat{H}}^{(l)}{\psi }_n^{(l)}(r)={\epsilon }_n^{(l)}{\psi }_n^{(l)}(r),n=1,2,\dots ,$$

$${\hat{H}}^{(l)}\equiv -\frac{1}{2}\frac{d^2}{d{r}^2}-\frac{1}{r}\frac{d}{dr}+V_{eff}([n,n_{+}];r)+\frac{l(l+1)}{2r^2}$$

$$V_{eff}([n,n_{+}];r)$$

$$n(r)=2\times \frac{1}{4\pi }\sum _{l,n}(2l+1)f_{l,n}{\left[{\psi }_n^{(l)}(r)\right]}^2$$

$$f_{l,n}\equiv f({\epsilon }_n^{(l)},\mu )=\frac{1}{\exp \{({\epsilon }_n^{(l)}-\mu )/\sigma \}+1}.$$

$$Nz=2\sum _{l,n}(2l+1)f_{l,n}.$$

$$F=E_{\mathrm{ISJM}}[n,\{n_{+}^{(\alpha )}\}]+2\sigma \sum _{l,n}(2l+1)\left[f_{l,n}\log f_{l,n}+(1-f_{l,n})\log (1-f_{l,n})\right]$$

$$V(r)=\left\{\begin{array}[]{lll}0&;&r<R\\ \infty&;&r\geq R,\end{array}\right.$$

$$W_s^{(l)}(r)={\left\{R^3{\left[j_{l+1}({\beta }_{l,s})\right]}^2/2\right\}}^{-1/2}{j}_l(\frac{{\beta }_{l,s}}{R}r),$$

$$j_l({\beta }_{l,s})=0,l=0,1,2,\dots ;s=1,2,3,\dots ,$$

$${\int }_{r=0}^{R}dr{r}^2{W}_s^{(l)}(r)W_{s^{\prime }}^{(l)}(r)={\delta }_{s,s^{\prime }},$$

$$E_{l,s}=\frac{{\beta }_{l,s}^2}{2R^2}.$$

$${\psi }^{(l)}(r)=\sum _s{a}_s^{(l)}{W}_s^{(l)}(r)$$

$${\left[{\hat{H}}^{(l)}\right]}_{s^{\prime },s}$$

Full text

I Introduction

Atomic clusters have attracted a great deal of attention in the past three decades because of their unique properties different from those of single atoms and bulk materials.de Heer et al. (1987); Jena et al. (1987, 1992); Haberland (1994); Ellis (1995); Jena and Behera (1996); Martin (1996); Sugano and Koizumi (1998); Jellinek (1999); Ekardt (1999); Meiwes-Broer (2000); Reinhard and Suraud (2004); Alonso (2005); Connerade and Solov’yov (2008); Philip et al. (2012); Stambula et al. (2014); Ditmire (2016); Zhang et al. (2016); Verkhovtsev et al. (2017) The physical properties of clusters depend on the atomic content and the geometry as well as the number of constituent atoms $N$. Nowadays, due to the significant progress in the experimental methods, it is possible to synthesize clusters with some desired properties. Among the various types of clusters, metal clusters (MC) are of much interest, firstly because of their biological and technological applications,Lang and Bernhardt (2012); Wang et al. (2013); Yin et al. (2012); Yin and Bernstein (2012); Yudanov et al. (2012); Collins et al. (2004); Shang et al. (2011); Lu and Chen (2012); Shimojo et al. (2010); Chen et al. (2013); Bonačić-Kouteckỳ et al. (2012) and secondly because the theoretical description of metallic systems are relatively simple.

In fact, with the first-time production of free alkali clusters,Knight et al. (1984) it was revealed that the quantized motion of the delocalized valence electrons in the average potential created by the ions is dominantly responsible for the shell effects observed in the abundance spectra. That is, the properties of alkali and other simple metal clusters are not much sensitive to the detailed structure of the ions.de Heer (1993) This fact led to the use of simple model, called jellium model (JM), in which the discrete ionic cores are replaced with a uniform positive background.Ekardt (1984a, b); Brack (1993)

Historically, development of the JM aimed to describe physical properties of metals within a small computational effort. It was extensively shown that this model was successful in describing some properties of simple s-p bonded metals.Ashcroft and Langreth (1967); Lang and Kohn (1970, 1971); Perdew and Monnier (1976, 1980); Perdew (1982); Kiejna and Wojciechowski (1996) However, rapid increasing of computational power along with the development of accurate quantum-chemical ab initio methods did not discourage researchers from finding simple models for accurate description of the materials’ properties. In contrast to the use of quantum-chemical ab initio methods which are appropriate for the accurate description of only small clusters (containing at most a few hundred atoms), using the JM enables one to go far beyond this limitation and describe the properties of very large simple MCs qualitatively and semi-quantitatively within an acceptable time period. In the first-time JM description of metal clusters, the positive ionic background was replaced by a jellium sphere of uniform charge in the electrostatic field of which the interacting valence electrons were moving. In subsequent works, spheroidal and ellipsoidal shapes for the jellium were considered which were useful in describing the odd-even effects seen experimentally. In one of the works, to describe the observed hollow clusters, the JM was applied using spherical shells of uniform charge.Polozkov et al. (2009)

Although the simple JM was successful in describing some effects observed in metal clusters, it was unable to predict accurate values for the energetics. It was also shown that simple JM, because of its mechanical instability, would lead to some inaccurate results for systems of very high or very low electron densities.Lang and Kohn (1970); Ashcroft and Langreth (1967)

Introducing the stabilized jellium model (SJM) was the first successful attemptPerdew et al. (1990) to overcome the above-mentioned drawbacks of the JM and get some accurate results for the bulk properties. In the SJM, the observed bulk density parameter, $r_{s}^{B}$, was used as input parameter for the background jellium density.

In applying the SJM for finite systems, because of the ionic relaxation near the surface, the $r_{s}^{B}$ would not bring the system to mechanical stability. In that case, the appropriate parameter, $r_{s}^{\dagger}$, which is close to $r_{s}^{B}$ for sufficiently large system, minimizes the total energy and brings the system to mechanical stability. The method was called self-compressed SJM (SC-SJM),Perdew et al. (1993); Payami (2001, 2004, 2006); Mahmoodi and Payami (2009) which could only predict an overall expansion or contraction for the whole system, and was not capable to explain the ionic relaxation in the surface region.

To explain the properties near the surface of a system, one should take a special care for the surface region, and to this end it is natural to assume that the whole system may be divided into two ”surface” and ”bulk” regions with two different jellium densities (the boundaries between regions are taken to be sharp), $\bar{n}_{1}$ and $\bar{n}_{2}$, as independent variables for the inner bulk and surface regions, respectively. This generalized model which was called the self-compressed inhomogeneous stabilized jellium model (SC-ISJM) has recently been formulated by the authors and successfully applied to describe the surface effects in simple metal thin films.Payami and Mahmoodi (2011)

In this work, we have applied the SC-ISJM to determine the surface- and bulk-region relaxations for closed-atomic-shell $N$-atom clusters of AlN ($N$=13, 19, 43, 55, 79, 87, 135, 141) as well as NaN and CsN ($N$=9, 15, 27, 51, 59, 65, 89, 113). Our results show that for most of them, the surface region contracts to higher electron and ionic densities, in good agreement with the first-principles method results. The inner regions for some of them expand to lower electron and ionic densities, and for others it contracts. That is, some clusters in their equilibrium states tend to have a relatively hollow shapes than being uniformly distributed over the space. Also, our calculations for Na${}_{65},$Cs51 and Cs65 clusters show that, in contrast to others, they undergo surface expansions. The results for Na65 and Cs65 are also in good agreement with those obtained from first-principles pseudo-potential method.Giannozzi et al. (2009) Comparison of the two sets of results reveals that the present model, in which the ionic constraints are somewhat released, give lower energies for the ground states and is more realistic than previously established ones.

The organization of this paper is as follows. In section II, the SC-ISJM is casted in the form suitable for clusters. Section III is dedicated to the calculation details, and the results are discussed in section IV. Finally, in section V, we conclude this work.

II Formulation of SC-ISJM for metal clusters

The total energy of a system in the SC-ISJM is given by Payami and Mahmoodi (2011)

[TABLE]

where,

[TABLE]

is the simple-JM energy, in which $T_{s}$ and $E_{xc}$ being the non-interacting kinetic and exchange-correlation energies, respectively; and the last term in the right hand side, the electrostatic energy, is given by

[TABLE]

with the positive charge distribution

[TABLE]

in which $\Theta^{(\alpha)}({{\bf r}})$ takes the value of unity inside the region $\alpha$, and zero outside. $\bar{n}_{\alpha}$ is the constant jellium density of region $\alpha$. All equations throughout the paper are expressed in Hartree atomic units.

In equation  (1), the averaged difference potential is given by

[TABLE]

where $\delta v$ is the difference potential between the pseudo-potentialAshcroft (1966) of ions and the electrostatic potential of the jellium background. $\bar{w}_{R}$ and $\varepsilon_{\rm M}$ are the average of the repulsive part of the pseudo-potential and the Madelung energy, respectively.

To determine the equilibrium state of the system with specified regions, the following minimization rule is applied:

[TABLE]

in which the inner minimization after the last equality sign gives the ground-state energy of the system (with $n_{0}$ as the ground-state electron density) constrained to the external parameters $\{n_{+}^{(\alpha)}\}$, and the outer minimization determines the equilibrium distribution of positive charges, $\{n_{+}^{(\alpha)\dagger}\}$, in the system.

For metal clusters, $\alpha$ takes the two values of 1 and 2 for the inner and surface regions, respectively. Assuming spherical geometry (which is appropriate for closed-atomic-shell clusters), the inner region is a sphere with radius $R_{1}$ and uniform density $\bar{n}_{1}$, while the surface region is a shell of thickness $t$ which is confined between two concentric spheres of radii $R_{1}$ and $R_{2}$ and having uniform density $\bar{n}_{2}$. The situation is schematically shown in Figure 1.

We may consider that the isolated cluster is cut out from an infinite bulk system. Before taking the cluster out of the bulk, the two above-mentioned regions have the same uniform bulk densities, $\bar{n}_{1}=\bar{n}_{2}=\bar{n}^{B}$, and the radius of the outer sphere is given by $R_{2}^{B}=(zN)^{1/3}r_{s}^{B}$ with $N$ being the number of constituent atoms of the cluster, $z$ being the valence of the atom, and $r_{s}^{B}=(3/4\pi\bar{n}^{B})^{1/3}$. The surface region is specified by the parameter $t=(R_{2}^{B}-R_{1}^{B})$ which could in principle take any continuous real number that should not be less than some critical value which corresponds to states with large repulsive forces between the ions. The starting structure of a closed-atomic-shell cluster is constructed from an infinite lattice by choosing a central atom and adding spherical shells of atoms successively. The first, second, etc., shells consist of first, second, etc., nearest-neighbor atoms, respectively, as shown in figure 2.

The result for aluminum is AlN clusters with $N$=13, 19, 43, 55, 79, 87, 135, 141 atoms, whereas for sodium and cesium are NaN and CsN clusters with $N$=9, 15, 27, 51, 59, 65, 89, 113 atoms. In figure 3, we have shown the typical closed-atomic-shell clusters Al87, Na59, and Cs65 resulted from the above procedure.

In applying the SC-ISJM for clusters, the minimization rule, Eq. (6), now reduces to

[TABLE]

in which $t^{\dagger}$, $\bar{n}_{1}^{\dagger}$, $\bar{n}_{2}^{\dagger}$ are the equilibrium surface thickness, inner- and surface-region jellium densities, respectively. The values of $r_{s}^{B}$, used here for Al, Na, and Cs are 2.07, 3.99, and 5.67 Bohrs.

III Calculation methods

The equilibrium states of closed-atomic-shell clusters are calculated using the SC-ISJM as well as the first-principles pseudo-potential method, which are explined as follows.

III.1 SC-ISJM

To find the equilibrium state of the cluster through Eq. (7), that is, minimization over a three-parameter space, we first take a trial surface thickness, $t$ (through which the two regions are specified), and calculate the ground-state energy for the given external parameters $t$, $\bar{n}_{1}^{B}$ and $\bar{n}_{2}^{B}$, using the density functional theoryHohenberg and Kohn (1964) (DFT) and the self-consistent solution of the Kohn-Sham (KS) equationsKohn and Sham (1965) within the local-density approximationKohn and Sham (1965) (LDA). The calculations then are repeated for different points in the {$\bar{n}_{1}$, $\bar{n}_{2}$} parameter space (keeping the number of atoms in each region as fixed) until a local minimum is achieved.Powell (1964) The global minimum-energy state is then selected from the calculated local-minimum states for different $t$ values, and thereby the equilibrium value, $t^{\dagger}$, is singled out. For $t$ we arbitrarily use the values obtained from $t=\nu d_{100}$ in which $\nu$ takes positive integer values and $d_{100}$ is the distance between two adjacent (100) planes, being 3.82, 4.05, and 5.72 for Al, Na, and Cs, respectively.

The constraint, that keeps constant the number of atoms in each region during the variations of jellium densities, is given by

[TABLE]

in which $\bar{n}_{\alpha}^{\dagger}$ and $V_{\alpha}^{\dagger}$ are the density and volume of region $\alpha$ in the equilibrium state, respectively; while $\bar{n}_{\alpha}$ and $V_{\alpha}$ are the corresponding quantities in an arbitrary non-equilibrium state. $V_{\alpha}^{B}$ corresponds to the volume of region $\alpha$ in its bulk un-relaxed state (before cutting out the cluster from the infinite bulk material). The constraints are easily expressed in terms of the densities of the two regions and the radii before and after relaxation as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Given the number of constituent atoms, $N$, and valence $z$, the outer bulk radius of the cluster, $R_{2}^{B}$, is determined from equation (9). Choosing a proper value for $t$, the value of the bulk inner radius, $R_{1}^{B}$, is determined from equation (10). The set of equations (9)-(12) completely specify the set of external parameters {$t,\bar{n}_{1},\bar{n}_{2}$}, from which the electrostatic potential energy of an electron for an arbitrary state is given by:

[TABLE]

In the SC-ISJM calculations, we assume a full rotational symmetry for the clusters, and therefore, the physical quantities will depend on the single radial variable $r$. Under this constraint, by taking the KS orbitals as a product of a radial function and spherical harmonics

[TABLE]

the KS equations reduce to one-dimensional eigenvalue equations:

[TABLE]

where,

[TABLE]

in which

[TABLE]

In this study, we have considered the closed-atomic-shell clusters. This does not imply that the resultant cluster is an electronic closed-shell one. On the other hand, if at the starting point ($\bar{n}_{1}=\bar{n}_{2}=\bar{n}^{B}$) we choose the number of constituent atoms in such a way that results in an electronic closed-shell state, then as soon as $\bar{n}_{\alpha}\neq\bar{n}^{B}$, the hierarchy of electronic levels differs from that in the starting configuration, and therefore, the relaxed cluster may become an electronic open-shell one. Moreover, in some cases where $\varepsilon_{\rm HOMO}\approx\varepsilon_{\rm LUMO}$, achieving the self-consistent point in the iterative solution of the KS equations becomes problematic. The workaround we adopted was to add a tiny “smearing” factor, $\sigma$, for the levels. That is, we consider the occupation of all states by the Fermi-Dirac weight. The resulting spherical averaged electron density then reads as

[TABLE]

in which $f_{l,n}$ is the Fermi-Dirac function,

[TABLE]

The chemical potential, $\mu$, which is close to the Fermi energy $\varepsilon_{\rm F}$ for a tiny $\sigma$ value, is determined from the solution of the equation

[TABLE]

This procedure is equivalent to treating the electrons in a heat bath of temperature $\sigma/k_{\rm B}$ ($k_{\rm B}$ is Boltzmann constant), and then the variational quantity would be the free energy given by

[TABLE]

in which the second term in the right hand side is the contribution of the entropy of the electronic energy levels. The Fermi-Dirac weight should also be included properly in the “band energy” (sum of KS eigenvalues) which is used in the calculation of the total energy.

To solve the KS equations for a cluster, we expand the KS orbitals in equation (15), in terms of basis functions, $W_{s}^{(l)}(r)$, which are the normalized radial eigenfunctions of an infinite spherical potential well,

[TABLE]

and are given by spherical Bessel’s functions as

[TABLE]

where, $\beta_{l,s}$ is the $s$th root of the spherical Bessel’s function of order $l$, that satisfies the equation

[TABLE]

and the radius of the well, $R$, is so chosen that $R\approx R_{2}^{B}+2\lambda_{\rm F}^{B}$ is satisfied ($\lambda_{\rm F}^{B}$ is the Fermi wavelength of the valence electrons). The basis functions satisfy the orthonormality condition

[TABLE]

and the eigenvalues of the infinite well are given by

[TABLE]

Inserting the expansion

[TABLE]

into equation (15), we get the matrix elements as

[TABLE]

from which we solve the secular equation and obtain the KS eigenvalues and eigenfunctions.

III.2 First-principles method

For the First-principles calculations, we have employed super-cell method with Martyna-Tuckerman correctionMartyna and Tuckerman (1999) for isolated systems as implemented in the Quantum-ESPRESSO code package.Giannozzi et al. (2009) The size of the cubic super-cell for each cluster is taken to be enough large suitable for an isolated cluster, and all calculations are done at the $\Gamma$ point. For the Al, Na, and Cs atoms, we have used the pseudo-potentials “Al.pz-vbc.UPF”, “Na.pz-n-vbc.UPF”, and “Cs.pz-nc.UPF”, respectively.Giannozzi et al. (2009) The calculations are based on the above-mentioned norm-conserving pseudo-potentials with a cutoff of 40 Ry for the plane-wave basis set and using smearing with 0.02 Ry of broadening. The geometries are fully relaxed until the total force on the cluster is less than 10*-5* Ry/au. Because of the symmetry inherent in the atomic positions of the closed-atomic-shell clusters, the relaxations take place in the radial direction of the atomic shells, and therefore, the changes in the inter-shell spacings can be compared with the SC-ISJM results.

IV Results and discussion

In this section, we present the calculation results of first-principles method and of the SC-ISJM method, and then discuss how the SC-ISJM results are consistent with those of ab initio method.

IV.1 First-principles results

In order to determine the finite-size effects on the clusters, we have first calculated the optimized lattice constants for bulk Al, Na and Cs crystals with fcc, bcc, and bcc structures, and obtained the values of 7.48, 7.79, and 11.75 a.u. for the equilibrium lattice constants, respectively. In the bulk calculations, the plane-wave basis set cutoffs for the wave functions and the density were chosen as 40 and 160 Ry, respectively. Also, for Brillouin-zone integrations, we used a mesh of size $20\times 20\times 20$ in the k space. We have then constructed our un-relaxed clusters out of these bulk-optimized structures as closed-atomic-shell ones containing successive shells of first-nearest-neighbors, second-nearest-neighbors, and so on, from the central atom. These structures are highly symmetric ones which may correspond not to a bound state of global-minimum energy, but this approximation makes them suitable for comparison with spherical SC-ISJM results.

The radial positions and numbers of atoms enclosed by shells for the un-relaxed clusters are shown by dashed blue lines; and the optimized atomic configurations which preserve the $n^{th}$-nearest neighbor structures but with some new radial distances, are represented by solid red lines in Figs. 4-6. Looking at the radial distance of the outermost shell (which we call the first shell) for relaxed systems, we see that all clusters, but Cs15 and Cs65, undergo size reductions. For the mentioned two exceptions the sizes are increased.

The results show that the second shells of all clusters; but Al43, Na27, Na51, Cs27, and Cs51; were also contracted. As to the third shell, all were contracted except for Al135, Na59, Cs51, Cs59, and Cs113 which were expanded or remained more or less unchanged. The behaviors of the innermost shells can also be inspected from the figures, but since we are mostly interested in the surface behaviors, we ignore those details and just mention their overall effects.

The relative distances between the first and second layers were reduced for all except Al87, Al135, Na15, Na65, Na113, Cs15, Cs65, Cs113, and for these it is increased or remained unchanged even though the mean radii might be changed. The relative distance between the second and third layers show reductions for all except Al43, Al141, Na27, Na51, Na89, Cs27, Cs51, Cs89. In these exceptions it is increased or remained more or less unchanged. Finally, the relative distances between the first and third layers were reduced for all except Na65 and Cs65. For these two exceptions, they were increased.

IV.2 SC-ISJM results

Using the SC-ISJM we have determined, through Eq. (7), the equilibrium surface region thickness, $t^{\dagger}=\nu^{\dagger}d_{100}$, and the corresponding equilibrium jellium densities, $\bar{n}_{1}^{\dagger}$, $\bar{n}_{2}^{\dagger}$, for AlN ($N$=13, 19, 43, 55, 79, 87, 135, 141), NaN and CsN ($N$=9, 15, 27, 51, 59, 65, 89, 113) clusters. The results are summarized in Table 1.

For comparison, the Table 1 also include the SC-SJM and SJM results. The values $\nu^{\dagger}=1,2$ imply that the surface region has approximately $t^{\dagger}=1\times d_{100}$ or $t^{\dagger}=2\times d_{100}$ thicknesses. The negativity of the values for $\Delta r_{s}2^{\dagger}$ implies that SC-ISJM predicts an increase for the ionic density of surface region. However, $\Delta r_{s}1^{\dagger}$ values are positive for all AlN, and mostly negative for NaN and CsN clusters. In cases that $\Delta r_{s}1^{\dagger}$ are positive, it implies that the inner-region atoms of the corresponding cluster tend to occupy a larger volume, which in turn indicates the expansion of inner region. However, this does not mean the overall expansion of the cluster because the sign of relaxation for the surface region must also be taken into account. These results show the advantages of SC-ISJM over SC-SJM in which all parts of jellium undergo contractions (see column 7 of Table 1).

In the SC-ISJM, the equilibrium sizes of all clusters, except for Al55, are smaller than their un-relaxed values. For Al55, the results show an expansion of the jellium. We notice that even though the surface region in this case is contracted ($\Delta r_{s}2^{\dagger}$=-0.068), the expansion of the inner region ($\Delta r_{s}1^{\dagger}$=+0.206) is dominated and resulted in an overall expansion.

For Al clusters, the SC-ISJM predicts expansions for inner regions, because $\Delta r_{s}1^{\dagger}>0$ for all of them. In NaN clusters, the inner regions of all, except for $N=9,15,59$, undergo contractions. For CsN clusters, the inner region of only Cs9 and Cs15 undergo expansions and the others show contractions.

In Figs. 7-9, we have plotted the equilibrium SC-ISJM electron and jellium densities of the clusters, and compared with the spherical averaged electron densities of first-principles calculation results. The unit for density is that of bulk jellium, $n^{\rm B}$, and the distances are expressed in terms of Fermi wavelength for bulk jellium, $\lambda_{\rm F}^{\rm B}$. The thick dashed green lines correspond to SC-ISJM electron densities; the thin dashed black lines and solid red lines correspond to jellium and first-principles electron densities, respectively. Jellium densities of less than unity correspond to expansion with respect to un-relaxed bulk jellium, whereas for larger than unity the system is contracted.

In Fig. 7 which corresponds to Al clusters, we observe that the inner regions are expanded for all clusters, and this expansion is significant for $N$=13, 19, 43, 55. On the other hand, these expansions are accompanied by the contractions of surface regions.

In Na clusters, as Fig. 8 shows, the expansion of inner region is realized for $N$=9, 15, 27 whereas for others it is contracted or remained unchanged. On the other hand, these changes are accompanied by contractions of surface regions. Careful inspection of Fig. 8 corresponding to Na65, we notice that even though both surface and inner regions are contracted with respect to the bulk, but the contraction of the surface region is less than that of the inner region. This behavior is in agreement with the relative surface-expansion seen in Fig. 5 corresponding to Na65.

Finally, in Fig. 9, we have plotted the SC-ISJM results for CsN clusters. For these systems in most cases we have $\nu^{\dagger}=2$, and for both regions we have jellium contractions, except for $N$=9, 15. Similar to Na65, here we have Cs51 and Cs65 which show relative expansions with respect to inner regions. Here, the Cs65 behavior is in agreement with the first-principles result in Fig. 6 which shows relative expansion of the surface layer.

As to electron densities, the spherical-averaged values of the first-principles calculations are quite similar to the SC-ISJM ones, both in extensions and in the form of Friedel oscillations.Friedel (1958) The discrepancy between SC-ISJM and first-principles electron densities near $r=0$ originates from the $r^{2}$ weight factor which is used in spherical average.

IV.3 Discussion

As was mentioned in Section I, the JM was developed to describe some electronic properties of metals within a low-cost computational effort. The SJM, which was a mature version of the JM, was able to describe the properties in a broader areas than the JM. However, to describe the geometric properties of finite metals, it was necessary to release the constraint of uniform positive background density, $\bar{n}^{\rm B}$, and let it take a new value, $\bar{n}^{\dagger}$, which stabilizes the system. The latter procedure which gave rise to the SC-SJM, was successful in giving the accurate sizes of simple metal clusters. The next development was to generalize the model in such a way to decouple the relaxations of surface and inner part of the system, which was called the SC-ISJM.

The SC-ISJM calculation results for AlN, NaN, and CsN clusters show that all, but Al55, undergo size reduction in good agreement with the atomic simulation results. However, in the simulations, the Cs15 and Cs65 were predicted to increase their sizes. This phenomenon which is shared in both simulation and the SC-ISJM, shows the advantages of the model and some miss-predictions can be attributed to the rough method of choosing the thickness as $t=\nu d_{100}$ or to the fact that assuming “only two regions” for the system is possibly not sufficient and one needs to straightforwardly generalize it to three or more regions.

The first-principles results showed that the second shells of most clusters were also contracted. This is also in agreement with the SC-ISJM results because, for some clusters the minimum-energy state occur at $\nu^{\dagger}$=2 and the sign of $\Delta r_{s}2^{\dagger}$ is negative. If in the application of SC-ISJM we would have considered more than two regions, we could then be able to reproduce the relaxations of the third or higher layers as in the simulations.

Finally, in the first-principles results, the relative distances between the first and second atomic shells were increased for Na15, Na65, Cs15, Cs65, and Cs113. This relative expansion was also reproduced in the SC-ISJM results for Na65, Cs51, and Cs65, which shows the good agreement of these results with those of atomic simulations.

The good agreement of the SC-ISJM results with those of first-principles method reveals that the model is more realistic than its predecessors and could even be improved.

V Conclusions

It is desirable to explain the properties of materials in simple ways by realistic models, because on the one hand it demands much less computational facilities and efforts, and on the other hand it is feasible by a personal computer system. In this respect, the JM and its variants were developed to describe the electronic structure properties of metals. The earlier versions of the JM were utilized to describe the electronic properties of mostly simple metal systems. The newer versions such as SC-SJM and SC-ISJM (subject of present work) can be used to describe some geometric properties either. As was discussed, the SC-SJM was able to predict only the overall size changes, but the SC-ISJM has the ability to reproduce the ionic inhomogeneity of a system. The SC-ISJM showed that it correctly predicts the size reduction of most clusters of AlN, NaN, and CsN, as was predicted by ab initio calculations. On the other hand, we showed that SC-ISJM can also predict the size enlargements in some clusters, as was seen from first-principles results. Finally, the expansions of inter-layer distances which was predicted by the atomic simulations, were correctly reproduced by the SC-ISJM. To overcome some discrepancies observed in the results of the model with those of first-principles method, one should improve the method of choosing surface thicknesses or straightforwardly generalize the model to include more than two regions.

Acknowledgements.

This work is part of research program in Theoretical and Computational Physics Group, NSTRI, AEOI.