Equilibrium configurations of large nanostructures using the embedded saturated-fragments stochastic density functional theory

TL;DR

This paper introduces an  ab initio Langevin dynamics method based on stochastic density functional theory with an embedded saturated fragment formalism, enabling scalable simulations of large covalently bonded nanostructures like silicon nanocrystals.

Contribution

It develops a linear-scaling, stochastic DFT-based Langevin dynamics approach applicable to large covalent systems using an embedded fragment formalism.

Findings

Successfully simulated silicon nanocrystals up to 3 nm in diameter.

Generated canonical configurations across a range of temperatures.

Analyzed surface structure and geometry reconstruction of nanocrystals.

Abstract

An \emph{ab initio} Langevin dynamics approach is developed based on stochastic density functional theory (sDFT) within a new \emph{embedded saturated } \emph{fragment }formalism, applicable to covalently bonded systems. The forces on the nuclei generated by sDFT contain a random component natural to Langevin dynamics and its standard deviation is used to estimate the friction term on each atom by satisfying the fluctuation\textendash dissipation relation. The overall approach scales linearly with system size even if the density matrix is not local and is thus applicable to ordered as well as disordered extended systems. We implement the approach for a series of silicon nanocrystals (NCs) of varying size with a diameter of up to $3$nm corresponding to $N_{e}=3000$ electrons and generate a set of configurations that are distributed canonically at a fixed temperature, ranging from cryogenic to room temperature. We also analyze the structure properties of the NCs and discuss the reconstruction of the surface geometry.

Tables

Table 1. Table 1: Value of various parameters for the LD based on efsDFT calculations: The friction coefficients γ 𝛾 \gamma , number of stochastic orbitals I 𝐼 I , time-step Δ ​ t Δ 𝑡 \Delta t and the wall time per single SCF iteration t iter subscript 𝑡 iter t_{\text{iter}} .

NC	T(K)	$\gamma \left({\text{fs}}^{-1}\right)$		$I$	$\Delta t\left({\text{fs}}^{-1}\right)$	$t_{\text{iter}}$$\left(\text{min}\right)$
		H	Si
${\text{Si}}_{35}{\text{H}}_{36}$	30	0.12	0.04	120	1.2	1
	300	0.04	0.04	30	1.2	1
${\text{Si}}_{147}{\text{H}}_{100}$	30	0.12	0.04	120	1.2	2
	300	0.12	0.04	30	1.2	2
${\text{Si}}_{705}{\text{H}}_{300}$	30	0.12	0.04	120	1.2	10
	300	0.12	0.04	92	1.2	10

Table 2. Table 2: The shells of the silicon NCs used for analyzing the bond length relaxation in Fig.  6 : Their inner and outer radii (in Å), the number of Si atoms N Si subscript 𝑁 Si N_{\text{Si}} , and the number of nearest neighbor (NN) Si-Si pairs N NN subscript 𝑁 NN N_{\text{NN}} .

Shell	$R_{\text{in}}$	$R_{\text{out}}$	$N_{\text{Si}}$	$N_{\text{NN}}$
$A$	0	5.5	35	52
$B$	5.5	8.5	113	158
$C$	9.0	11.6	153	163
$D$	11.6	13.6	200	189
$E$	13.6	15.1	205	168

Equations

$${\hat{h}}_{\text{KS}}{\varphi }_n\left(\mathbf{r}\right)={\epsilon }_n{\varphi }_n\left(\mathbf{r}\right),$$

$${\hat{h}}_{\text{KS}}=-\frac{{\hslash }^2}{2m_e}{\nabla }^2+v_{KS}\left(\mathbf{r}\right),$$

$$v_{KS}\left(\mathbf{r}\right)=v_{eN}\left(\mathbf{r}\right)+v_H\left(\mathbf{r}\right)+v_{xc}\left(\mathbf{r}\right).$$

$$n\left(\mathbf{r}\right)=2\sum _n\theta \left(\mu -{\epsilon }_n\right){\left|{\varphi }_n\left(\mathbf{r}\right)\right|}^2,$$

$$n\left(\mathbf{r}\right)=2\text{Tr}\left[\theta \left(\mu -{\hat{h}}_{\text{KS}}\right)\delta \left(\mathbf{r}-\hat{\mathbf{r}}\right)\right].$$

$$n^I\left(\mathbf{r}\right)$$

$${\times \sqrt{{\theta }_{\beta }\left(\mu -{\hat{h}}_{\mathrm{KS}}\right)}|\chi ⟩⟩}_{\chi }$$

$$⟨\mathbf{r}|\chi ⟩=h^{-3/2}{e}^{i{\phi }_{\mathbf{r}}}$$

$${\mathbf{f}}_{\alpha }=-\int \frac{\partial v_{eN}\left(\mathbf{r}\right)}{\partial {\mathbf{R}}_{\alpha }}{n}^I\left(\mathbf{r}\right)d^{3}r.$$

$${\mathbf{f}}_{\alpha }={\mathbf{f}}_{\alpha }^{\text{det}}+{\mathbf{f}}_{\alpha }^{\text{fluc}}+{\mathbf{f}}_{\alpha }^{\text{bias}}$$

$$n\left(\mathbf{r}\right)=n_F\left(\mathbf{r}\right)+\Delta n\left(\mathbf{r}\right)$$

$$n\left(\mathbf{r}\right)=n^I\left(\mathbf{r}\right)+\sum _{f=1}^F\left(n_f\left(\mathbf{r}\right)-n_f^I\left(\mathbf{r}\right)\right),$$

$${\mu }_{\alpha }{\ddot{\mathbf{q}}}_{\alpha }={\mathbf{f}}_{\alpha }\left(\mathbf{q}\right)-{\gamma }_{\alpha }{\mathbf{p}}_{\alpha }+{\mathbf{\eta }}_{\alpha }.$$

$$⟨\mathbf{\eta }{}_{\alpha }\left(t\right)⟩$$

$$⟨\left(\mathbf{\eta }{}_{\alpha }\left(t\right)+{\mathbf{f}}_{\alpha }^{\text{fluc}}\left(t\right)\right)\otimes \left(\mathbf{\eta }{}_{{\alpha }^{\prime }}\left(t^{\prime }\right)+{\mathbf{f}}_{{\alpha }^{\prime }}^{\text{fluc}}\left(t^{\prime }\right)\right)⟩$$

$$⟨\mathbf{\eta }{}_{\alpha }\left(t\right)\otimes \mathbf{\eta }{}_{{\alpha }^{\prime }}\left(t^{\prime }\right)⟩+⟨{\mathbf{f}}_{\alpha }^{\text{fluc}}\left(t\right)\otimes {\mathbf{f}}_{{\alpha }^{\prime }}^{\text{fluc}}\left(t^{\prime }\right)⟩$$

$$\text{I3×3}{\sigma }_{\alpha }^2{\delta }_{\alpha {\alpha }^{\prime }}\delta \left(t-t^{\prime }\right),$$

$${\sigma }_{\alpha }^2=2{\mu }_{\alpha }{\gamma }_{\alpha }{k}_{\mathrm{B}}T.$$

$${\mathbf{q}}_{\alpha }^{m+1}$$

$${\mathbf{p}}_{\alpha }^{m+1}$$

$$T_K^m=\frac{2}{3N{k}_{\mathrm{B}}}\sum _{\alpha }{\left({\mathbf{p}}_{\alpha }^m\right)}^2/2{\mu }_{\alpha }$$

$$T_V^m=-\frac{1}{3N{k}_{\mathrm{B}}}\sum _{\alpha }\left({\mathbf{f}}_{\alpha }^m+{\mathbf{\eta }}_{\alpha }^m\right)\cdot \left({\mathbf{q}}_{\alpha }^m-⟨{\mathbf{q}}_{\alpha }⟩\right).$$

$$⟨\hat{A}⟩={⟨\hat{A}⟩}^I+\sum _{f=1}^F\Delta A_f^I,$$

$$\Delta A_f^I=⟨{\hat{A}}_f⟩-{⟨{\hat{A}}_f⟩}^I,$$

$$n_f\left(\mathbf{r}\right)=c_f{\left(\mathbf{r}\right)}^2\sum _n{\left(p_n^f\right)}^2{\psi }_n^f{\left(\mathbf{r}\right)}^2,$$

$$c_f(\mathbf{r})=\sqrt{\frac{{\sum }_{a\in f}{n}_a^{\left(0\right)}\left(\mathbf{r}\right)}{{\sum }_{a\in sf}{n}_a^{\left(0\right)}\left(\mathbf{r}\right)}},$$

$$\int n_f\left(\mathbf{r}\right)d\mathbf{r}=\int \sum _{a\in f}{n}_a^{\left(0\right)}\left(\mathbf{r}\right)d\mathbf{r}.$$

$$2\sum _n⟨{\tilde{\psi }}_n^f|{\tilde{\psi }}_n^f⟩=\int \sum _{a\in f}{n}_a^{\left(0\right)}\left(\mathbf{r}\right)d\mathbf{r}.$$

$$n_f\left(\mathbf{r}\right)=2\sum _m{\varphi }_m^f{\left(\mathbf{r}\right)}^2$$

$${\hat{\theta }}_f=2\sum _m|{\varphi }_m^f⟩⟨{\varphi }_m^f|.$$

$$⟨{\hat{A}}_f⟩\equiv tr\left[{\hat{\theta }}_f\hat{A}\right]=tr\left[\sqrt{{\hat{\theta }}_f}\hat{A}\sqrt{{\hat{\theta }}_f}\right],$$

$$\sqrt{{\hat{\theta }}_f}=\sqrt{2}\sum _m{\left(s_m^f\right)}^{-1/2}|{\varphi }_m^f⟩⟨{\varphi }_m^f|.$$

Full text

I Introduction

*Ab initio *molecular dynamics based on density functional theory (DFT) is becoming an important tool for studying the plethora of structural and dynamical processes in a broad range of systems in material science, chemistry, biology and physics.Car and Parrinello (1985); Barnett et al. (1991); Payne et al. (1992); Kresse and Hafner (1993); Tuckerman et al. (1995); Marx and Hutter (2000); Schlegel et al. (2001); Herbert and Head-Gordon (2005); Raty, Gygi, and Galli (2005); Bockstedte et al. (1997); Cawkwell, Niklasson, and Dattelbaum (2015) The application of this approach to very large systems is still limited by the computational scaling of the electronic structure portion of the calculation, regardless of whether one uses a Lagrangian-based or Born-Oppenheimer-based methods. This is because of the cubic scaling involved in solving the Kohn-Sham equations coupled with the need to iterate to self-consistency or to propagate the Kohn-Sham (KS) orbitals, as both of these options further increases the computational times by an order of magnitude.

Significant advances in these respects have been made along two major directions. One primary direction is based on a Lagrangian formulation of density functional theory Car and Parrinello (1985); Cawkwell, Niklasson, and Dattelbaum (2015) and circumvents the need for SCF iterations by propagation of the KS orbitals. This venue does not eliminate the cubic scaling and is therefore limited to relatively small systems. Another approach is based on linear-scaling techniques Goedecker (1999); Gillan et al. (2007); Baer and Head-Gordon (1997a); Soler et al. (2002), that reduces the algorithmic complexity by finding the density matrix directly, relying on its asymptotic sparseness in real-space. However, sparsity sets in only for very large systems, limiting the applicability sparse-matrix methods, especially in 3D.

In a recent set of papers we have introduced the stochastic DFT (sDFT) methods Baer, Neuhauser, and Rabani (2013); Neuhauser, Baer, and Rabani (2014); Gao et al. (2015); Neuhauser et al. (2015) which scales linearly (or even sublinearly) with the system size and does not rely on the sparsity of the density matrix. sDFT is a general approach to electronic structure based on a stochastic process and is applicable to extended ordered as well as disordered materials. Some of the techniques we use, based on the stochastic trace formula.Hutchinson (1990), have been developed for tight-binding electronic structure Drabold and Sankey (1993); Wang (1994); Röder et al. (1997), for molecular electronics Baer et al. (2004) and for multi-exciton generation in nanocrsytals Baer and Rabani (2012). The success of sDFT in reducing the scaling comes at a price of introducing a stochastic error in all its predictions, including forces, and that precludes application to *ab initio *molecular dynamics.

In this paper we show that sDFT can be used to study equilibrium structural properties of large NCs, despite the statistical fluctuations in the force estimates. For this, we invoke the Langevin equation following the work of Attaccalite and Sorella,Attaccalite and Sorella (2008) and generate a sequence of configurations distributed according to the canonical ensemble. These configurations can be used in a variety of applications for studying the structural, electronic and optical properties of NCs. Here we demonstrate their use for studying the structural properties of silicon nanocrystals (NCs) with a diameter of up to $3$ nm, and $N_{e}=3000$ electrons.

The paper includes development of the embedded saturated fragments method which allows reducing the statistical errors in sDFT. This new method is inspired by, but more general than, the embedded fragments method developed in Ref. Neuhauser, Baer, and Rabani, 2014. It uses small saturated fragments of the system, and carves out the relevant part of the density to be embedded in the system. Hence, it is applicable not only to clusters of molecules, like Ref. Neuhauser, Baer, and Rabani, 2014, but also to covalently-bonded systems such as silicon NCs. The method is described in detail in Appendix A.

II Methods

II.1 Stochastic DFT

Kohn-Sham density functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965) maps a system of $N_{e}$ interacting electrons in an external electron-nucleus potential $v_{eN}\left(\mathbf{r}\right)=-\frac{e}{4\pi\epsilon_{0}}\sum_{N}\frac{Z_{N}e}{\left|\mathbf{r}-\mathbf{R}_{N}\right|}$, where $\mathbf{R}_{N}$ ($N=1,2,...$) are the nuclei positions and $Z_{N}e$ are their charge ($e$ is the electron charge), onto a system of *non-interacting *electrons (the KS system), having the same ground-state density $n\left(\boldsymbol{r}\right)$. This mapping is performed by solving the KS equations Hohenberg and Kohn (1964); Kohn and Sham (1965)

[TABLE]

where the KS Hamiltonian is:

[TABLE]

and the KS potential $v_{KS}\left(\mathbf{r}\right)$ is the sum of the external electron-nuclear potential $v_{eN}\left(\mathbf{r}\right)$, the density-dependent Hartree potential $v_{H}\left(\mathbf{r}\right)=\frac{e^{2}}{4\pi\epsilon_{0}}\int\frac{n\left(\mathbf{r}^{\prime}\right)}{\left|\mathbf{r-r^{\prime}}\right|}d\mathbf{r}^{\prime}$, and the exchange-correlation potential $v_{xc}\left(\mathbf{r}\right)$:

[TABLE]

In the KS system, the density is expressed in terms of the normalized single electron KS eigenstates $\phi_{n}\left(\boldsymbol{r}\right)$ and eigenvalues $\varepsilon_{n}$

[TABLE]

where $\theta\left(x\right)$ is the Heaviside function and $\mu$ is the chemical potential chosen so that $2\sum_{n}\theta\left(\mu-\varepsilon_{n}\right)=N_{e}$. Eqs. (1)-(4) must be solved self-consistently, since $\hat{h}_{\text{KS}}$ depends on the density. While the entire scheme is a significant simplification over the original many-electron problem, it remains a challenge for large systems since the computational effort scales as $O\left(N_{e}^{3}\right)$.

An important step towards reducing the computational scaling of KS-DFT was recently proposed by Baer*,* Neuhauser, and Rabani (BNR),Baer, Neuhauser, and Rabani (2013) where the density of Eq. (4) was expressed as a trace over the projected density operator:Baer, Neuhauser, and Rabani (2013)

[TABLE]

The problem now shifts into calculating self-consistently the trace in Eq. (5) (since $\hat{h}_{\text{KS}}$ depends on $n\left(\boldsymbol{r}\right)$) rather than solving the KS equations by brute-force diagonalization. When the trace is performed using the KS eigenstates, the computational cost remains $O\left(N_{e}^{3}\right)$ similar to the traditional approach. However, since the trace is invariant to the basis, alternative schemes that potentially lead to improved scaling can be used. One such scheme is based on the concept of a stochastic trace formula, which reduces the scaling of the trace operation by introducing a controlled statistical error.Hutchinson (1990)

Using the stochastic trace formula, the density can be estimated as a symmetrized stochastic trace formula, given by:Baer, Neuhauser, and Rabani (2013)

[TABLE]

where $\left\langle\cdots\right\rangle_{\chi}$ denotes an average over $I$ stochastic orbitals $\left|\chi\right\rangle$, defined as:

[TABLE]

for each grid point $\boldsymbol{r}$, the parameter $h$ (not to be confused with the KS Hamiltonian $\hat{h}_{KS}$ operator) is the grid spacing, and $\varphi_{\boldsymbol{r}}$ are statistically independent random variables in the range $[0,2\pi]$ ($\left\langle e^{i\varphi_{\boldsymbol{r}}}e^{-i\varphi_{\boldsymbol{r}^{\prime}}}\right\rangle_{\varphi}=\delta_{\boldsymbol{r}\,\boldsymbol{r}^{\prime}}$). The density $n\left(\boldsymbol{r}\right)$ is, strictly speaking, given by the limit $n\left(\boldsymbol{r}\right)=\lim_{I\rightarrow\infty}n^{I}\left(\boldsymbol{r}\right)$ and we approximate it with a finite $I$. The Heaviside function in Eq. (6) is smoothed by the function $\theta_{\beta}\left(\varepsilon\right)\equiv\frac{1}{2}\mbox{erfc}\left[\beta\varepsilon\right]$, where $\beta$ is a large constant satisfying $\beta E_{g}\gg 1$, where $E_{g}$ is the KS-DFT fundamental gap. Throughout this paper we set the value of $\beta$ to $100E_{h}^{-1}$. The action of $\sqrt{\theta_{\beta}\left(\mu-\hat{h}_{{\rm KS}}\right)}$ on $|\chi\rangle$ is evaluated by a Chebyshev expansion in powers of the sparse KS Hamiltonian, $\hat{h}_{{\rm KS}}$.Tal-Ezer and Kosloff (1984) The length of the Chebyshev series is determined by the value of $\left(\mu-E_{min}\right)/\Delta E$ and $\beta\Delta E$ where $\Delta E=\left(E_{max}-E_{min}\right)/2$ and $E_{min/max}$ are the minimal and maximal eigenvalues of $\hat{h}_{KS}$. Under the conditions of the present systems, the length of the series is ~3000 terms.

The stochastic trace evaluation (Eq. 6) reduces the computational scaling of KS-DFT to $O\left(N_{e}\right)$ and for certain properties even to a sub-linear scaling.Baer, Neuhauser, and Rabani (2013) Linear-scaling complexity is achieved due to the following facts:

1. the application of a Hamiltonian to a stochastic orbital, $\hat{h}_{KS}\left|\chi\right\rangle$ requires a linear scaling effort (irrespective of the structure of the orbital); 2) The length of the Chebyshev series is independent (or at most weakly dependent) of system size and 3) Only a system-size independent number of stochastic orbitals are required. This type of assumptions is different from the linear-scaling approaches depending on density matrix sparsity,Baer and Head-Gordon (1997b); Goedecker (1999) which assume that locality of orbitals is not completely destroyed by the repeated operation of the Hamiltonian.

A converged self-consistent solution of Eq. (6) provides an estimate of the electron density and in addition can be used to generate other quantities, such as the density of states (DOS), the total energy per electron, and the forces acting on the nuclei. All estimates contain a statistical error that can be controlled by increasing the number of stochastic orbitals ($I$) used to evaluate the trace in Eq. (6). Of particular relevance to this work are the Cartesian forces exerted by the electrons on $N$ nuclei ($\alpha=1,\dots,N$), which can be evaluated through the Hellmann-Feynman theorem:Hellman (1937); Feynman (1939)

[TABLE]

It should be stressed that for finite sampling, these forces are only approximately commensurate with the stochastic estimate of the energy (which is not used in the sampling procedure at all), as discussed in Appendix B. The stochastic estimate of the Hellmann-Feynman forces is an excellent estimator of the deterministic forces, as can be seen in Fig. 1, where the mean absolute deviation is dominated by the fluctuations and not by additional bias terms.

These sDFT forces can be expressed as:

[TABLE]

where $\boldsymbol{f}_{\alpha}^{\text{det}}$ is the deterministic (generally unknown) force, $\boldsymbol{f}_{\alpha}^{\text{fluc}}$ is the pure fluctuating term, and $\boldsymbol{f}_{\alpha}^{\text{bias}}$ is the bias expected to be proportional to $\frac{1}{I}$ in leading order. The choice of $I$ should be large enough to reduce $\boldsymbol{f}_{\alpha}^{\text{bias}}$ to negligible values and the only source of error in the procedure is then the statistical fluctuations proportional to $\frac{1}{\sqrt{I}}$ with vanishing mean ($\left\langle\boldsymbol{f}_{\alpha}^{\text{fluc}}\right\rangle=0$).

II.2 Embedded saturated fragments sDFT

The reduction of the scaling in sDFT is achieved by replacing the deterministic, numerically exact, trace evaluation with a stochastic sampling of the density. In return, this leads to statistical errors in the computed observables. To reduce the size of the statistical fluctuations, an embedded saturated fragments method is introduced inspired by (but different from) the method of Ref. Neuhauser, Baer, and Rabani, 2014 . The latter approach was suitable mainly for systems composed of proximate but chemically separated molecules (like clusters of water molecules, for examples). The present method is applicable for fragmenting covalently bonded systems, like silicon NCs.

In this approach, the system is divided into $F$ small fragments that are possibly overlapping. The division to fragments is flexible, and any desired physically motivated fragmentation can be used. The density is then a sum of the fragment density and a small correction term:

[TABLE]

where $n_{F}\left(\boldsymbol{r}\right)=\sum_{f=1}^{F}n_{f}\left(\boldsymbol{r}\right)$ is the density generated by the individual fragments obtained from a deterministic KS-DFT calculation for each fragment and $\Delta n\left(\boldsymbol{r}\right)=\left(n^{I}\left(\boldsymbol{r}\right)-n_{F}^{I}\left(\boldsymbol{r}\right)\right)$ is a correction term evaluated using stochastic orbitals. Here, $n^{I}\left(\boldsymbol{r}\right)$ is given by Eq. (6) and $n_{F}^{I}\left(\boldsymbol{r}\right)=\sum_{f=1}^{F}n_{f}^{I}\left(\boldsymbol{r}\right)$ is a sum over a *stochastic *estimate of the fragments density. In the limit $I\rightarrow\infty$, Eqs. (6) and (10) are identical and equal to the deterministic density. For finite values of $I$, the size of the statistical fluctuations of the two approaches are quite different. Since the deterministic fragmented density, $n_{F}\left(\boldsymbol{r}\right)$, provides a reasonable approximation for the full density $n\left(\boldsymbol{r}\right)$, the correction term, $\Delta n\left(\boldsymbol{r}\right)$, which is evaluated stochastically, is rather small, leading to a reduced variance in the relevant observables (forces, DOS, total energy per electron, etc.) compared to the direct stochastic approach of Eq. (6). An equivalent viewpoint is that the fragmentation is a device for reducing the variance in the stochastic evaluation of the density. This is evident by rewriting Eq. (10) in the following form

[TABLE]

and the implementation of this form is described in Appendix (A).

To assess the accuracy of the embedded saturated fragmented sDFT (efsDFT), we calculated the standard deviations (STDs) and mean absolute deviations with respect the deterministic DFT (MADs) of the atomic forces in a $\text{Si}{}_{35}\text{H}_{36}$ NC using hydrogen passivated $\text{Si}{}_{5}$ fragments. The results are shown in Fig. 1.111All calculations in this work use real-space grids of spacing $\Delta x=0.5a_{0}$, Troullier-Martins norm-conserving pseudopotentials Troullier1991 within the Kleinman-Bylander approximation.Kleinman1982 Fast Fourier Transforms were used for applying the kinetic energy operator and for determining the Hartree potentials and the method of Ref. Martyna1999 was used for treating the long range Coulomb interactions in a finite simulation cell with periodic boundary conditions. DFT calculations were performed under the local density approximation (LDA). The STDs and MADs decrease as $1/\sqrt{I}$, indicating that the bias in the force estimation is negligible. The standard deviations in the sDFT forces are larger by a factor of $\approx 3$ compared to those of efsDFT. This implies that the required number of stochastic orbitals in efsDFT is nearly an order of magnitude smaller than in sDFT for similar STDs. The STDs can be further reduced by using larger fragments as discussed below (cf., Fig. 7).

II.3 Langevin dynamics based on

efsDFT

The standard approach to generate canonically distributed configurations using ab initio techniques is based on molecular dynamics, which requires as input accurate force estimates for each atomic degree of freedom. Since the forces generated by efsDFT contain a stochastic component, we use Langevin dynamics (LD) instead of molecular dynamics to sample configurations according to the Boltzmann distribution. A LD trajectory Allen and Tildesley (1987); Van Kampen (1992); Zwanzig (2001); Frenkel and Smit (2002) is a sequence of configurations $\left(\boldsymbol{p},\boldsymbol{q}\right)^{m}=\left(\boldsymbol{p}\left(t_{m}\right),\boldsymbol{q}\left(t_{m}\right)\right)$ at discrete “times” $t_{m}=m\Delta t$, where $\Delta t$ is the time step, and $\boldsymbol{q}\equiv\left(\boldsymbol{q}_{1},\dots,\boldsymbol{q}_{N}\right)$ and $\boldsymbol{p}\equiv\left(\boldsymbol{p}_{1},\dots,\boldsymbol{p}_{N}\right)$ are the Cartesian coordinates and conjugate momenta, respectively, for the $N$ atoms, The trajectory is a solution of the Langevin equation (LE) of motion:Langevin (1908)

[TABLE]

where $\mu_{\alpha}$ is the mass of the atom $\alpha$, $\gamma_{\alpha}$ is its friction constant, and $\boldsymbol{f}_{\alpha}=\boldsymbol{f}_{\alpha}^{\text{det}}+\boldsymbol{f}_{\alpha}^{\text{fluc}}$ is the total efsDFT force acting on it, including deterministic and fluctuating parts (see Eq. 9). The bias is assumed negligible, so that $\left\langle\boldsymbol{f}_{\alpha}\right\rangle=\boldsymbol{f}_{\alpha}^{\text{det}}$ . In Eq. (12) $\boldsymbol{\eta}_{\alpha}$ is an additional uncorrelated white-noise force introduced so as to satisfy the fluctuation-dissipation (FD) relation. We require that the total random fluctuation on each atom obey:

[TABLE]

and

[TABLE]

where $\alpha,\alpha^{\prime}=1,\dots N$ are atom indices, $\left\langle\cdots\right\rangle$ designates average over the atomic force distribution, $\boldsymbol{\text{I}}_{3\times 3}$ is the $3\times 3$ unit matrix, and $\sigma_{\alpha}$ is the atomic force STD of atom $\alpha$, which is taken to satisfy the fluctuation-dissipation relation:

[TABLE]

We use the Verlet-like algorithm Grønbech-Jensen and Farago (2013) for numerically integrating the LE of motion at a fixed temperature $T$ and a predefined time-step $\Delta t$. The positions and momenta in time step $m+1$ depend on the positions and momenta in time step $m$ as well as on the forces in time step $m$ and the additional white noise $\boldsymbol{\eta}_{\alpha}^{m}$ is sampled from a Gaussian distribution such that the discretized version of Eq. (13) holds: \left\langle\left(\boldsymbol{\eta}{}_{\alpha}^{m}+\boldsymbol{f}_{\alpha}^{\text{fluc}}\right)\otimes\left(\boldsymbol{\eta}{}_{\alpha^{\prime}}^{n}+\boldsymbol{f}_{\alpha^{\prime}}^{\text{fluc}}\right)\right\rangle\Delta t=\text{\boldsymbol{\text{I}}}\sigma_{\alpha}^{2}\delta_{\alpha\alpha^{\prime}}\delta_{mn}:

[TABLE]

where $a_{\alpha}=b_{\alpha}\left(1-\frac{1}{2}\gamma_{\alpha}\Delta t\right)$ and $b_{\alpha}^{-1}=1+\frac{1}{2}\gamma_{\alpha}\Delta t$. The algorithm allows for stable and accurate solutions of the LE with time step comparable to that used in molecular dynamics simulations for similar systems. It treats the additional white noise component $\boldsymbol{\eta}_{\alpha}$ of the force differently from the force $\boldsymbol{f}_{\alpha}=\boldsymbol{f}_{\alpha}^{\text{det}}+\boldsymbol{f}_{\alpha}^{\text{fluc}}$ that result from the efsDFT calculation containing deterministic and fluctuating components that cannot be separated.

In Fig. 2 we plot for $\text{Si}_{35}\text{H}_{36}$ the running average of the transient temperature, $T^{m}$, calculated from the kinetic energy

[TABLE]

and from the virial estimator,

[TABLE]

In the above, $\left\langle\boldsymbol{q}_{\alpha}\right\rangle$ is the time average of the coordinate of atom $\alpha$. The initial positions of the Si atoms were taken from the bulk values. All surface Si atoms with more than two dangling bonds were removed and the remaining surface Si atoms were passivated using one or two H atoms placed in a tetrahedral geometry at the Si-H distance of $1.47$Å. The momenta were sampled from a Boltzmann distribution at $T=300K$. This non-equilibrium initial configuration relaxes towards equilibrium.

The agreement in Fig. 2 between the two temperature estimators is consistent with a proper sampling of the canonical distribution of both positions and velocities. The small discrepancies at the longest averaging time are due to the large fluctuations of the transient temperature, particularly when using the virial estimator. We have also calculated the fluctuations in the kinetic energy and found good agreement with the corresponding analytical value (see caption of the figure). The two atomic species have a slightly $\pm 5\%$ deviations in the temperatures. These may result from several factors such as the finite timestep of the Langevin propagator Grønbech-Jensen and Farago (2013), incomplete SCF convergence,Martínez et al. (2015) insufficiently accurate estimate of the amount of white noise $\eta_{\alpha}^{m}$ in Eq. (13) required to fulfill the fluctuation-dissipation relation.

II.4 Determining the optimal friction

The effect of $\gamma_{\text{Si}}$ on the configurational relaxation and on the velocity autocorrelation decay is illustrated in Fig. 3 for $\text{Si}_{36}\text{H}_{35}$. In order to decrease the number of unknown parameters we set the values of $\gamma_{\alpha}$ and $\sigma_{\alpha}$ to be identical for all atoms of the same type (i.e. Si or H in the systems studied here). To achieve such a uniform value of $\sigma$ we introduced white noise $\eta_{\alpha}^{i}=\sqrt{\left(\sigma_{\alpha}^{i}\right)^{2}-\left\langle\left(f_{\alpha}^{i,\text{fluc}}\right)^{2}\right\rangle}$ for each degree of freedom (see Eq. 13), where $\left\langle\left(f_{\alpha}^{i,\text{fluc}}\right)^{2}\right\rangle$ is estimated by a separate set of runs on the initial NC configuration using several independent sets of stochastic orbitals. Note, that we have tested that the magnitude of the sDFT force fluctuation $\left\langle\left(f^{\text{fluc}}\right)^{2}\right\rangle$ is not sensitive to the particular configuration used.

As expected, the configurational relaxation time increases with increasing values of $\gamma_{\text{Si }}$ with the opposite trend for the decay time of the velocity autocorrelation function. Based on the results of $\text{Si}_{36}\text{H}_{35}$ presented in Fig. 3 we conclude that a friction coefficient of $0.04\text{fs}^{-1}$ is sufficiently small for this system, with respect to minimizing both velocity and pair distance autocorrelation times. Although lower values of the friction coefficients could decrease the correlation time further, they would require reducing the statistical noise, which would be expensive to achieve using sDFT. Thus we chose the friction constants, $\gamma_{Si}=0.04\text{fs}^{-1}$ for Silicon and $\gamma_{H}=0.12\text{fs}^{-1}$ for the lighter H atoms. These values were used for the larger systems described in the next section (see Table 1). Note that the results shown in Fig. 3, which were generated using LD under dDFT, could have been equally well generated under efsDFT. This is shown explicitly for $\gamma_{\text{Si}}=0.04\,\text{fs}^{-1}$ (dotted red line) proving that the relaxation times are similar to those of the dDFT based LD calculation with the same value of $\gamma_{\text{Si}}$.

II.5 Validation of LD within efsDFT

Validation of the structure obtained using efsDFT based LD is demonstrated using the pair distribution function $g\left(r\right)$.Allen and Tildesley (1987) For finite size NCs the average number of neighbors at a distance $r$ is expected to be smaller than the bulk value due to surface atoms with a smaller number of neighbors. Fig. 4 shows a close agreement between the dDFT and efsDFT based LD estimates of $g\left(r\right)$ of the $\text{Si}_{35}\text{H}_{36}$ NC at two temperatures. The inset focuses on the first peak in $g\left(r\right)$, comparing the efsDFT to dDFT at $T=30$ and $300\text{K}$ .

III Results

In the previous sections we presented the methods and assessed the accuracy and validity of the efsDFT based LD. Here we apply the method to study structural properties of larger NCs exceeding $N_{e}=3000$ electrons. The Si-Si pair distribution functions $g\left(r\right)$ at two temperatures $T$ ($30$ and $300\text{K}$) are displayed in the upper panel of Fig. 5 for $\text{Si}_{147}\text{H}_{100}$ and $\text{Si}_{705}\text{H}_{300}$. Temperature broadens the peaks by a factor of $2-3$ without significantly changing the peak position.

In the lower panel of Fig. 5 we plot the transient and relaxed $g\left(r\right)$ at $30\text{K}$ for the two systems, focusing on the first, nearest neighbor peak. As described also for $\text{Si}_{35}\text{H}_{36}$, the initial positions of the Si atoms for both systems were taken from the experimental bulk values and all surface Si atoms with more than two dangling bonds were removed. The remaining surface Si atoms were then passivated using one or two H atoms placed in a tetrahedral position at the Si-H distance of $1.47$Å. The initially sharp peak broadens and shifts to longer Si-Si bond lengths as the system relaxes towards thermal equilibrium. For $30\text{K}$, the relaxation times are $180$ and $650\,$fs for $\text{Si}_{147}\text{H}_{100}$ and $\text{Si}_{705}\text{H}_{300}$, respectively. For $300\text{K}$ they are $180$ and $250\,$fs respectively.

The relaxation transient is studied in greater detail in Fig.6, where the average nearest-neighbor bond lengths are shown for $\text{Si}_{147}\text{H}_{100}$ (spherical shells A-B) and $\text{Si}_{705}\text{H}_{300}$ (shells A-E); see Table 2 for the definition and properties of the shells. In $\text{Si}_{705}\text{H}_{300}$ the deep layer shells (A-D) relax slower than those near the surface showing that relaxation progresses from the surface inwards. The difference between the relaxation times of the two systems is correlated with the smaller frequency, $\omega$, of the breathing mode of the larger NC. In the limit of an over damped motion (as is the case here since $\gamma^{2}\gg\omega^{2}$), the relaxation is dominated by two timescales proportional to $\gamma^{-1}$ and $\left(\omega^{2}/\gamma\right)^{-1}$. The former leads to a fast relaxation while the latter is slower and depends on the value of $\omega^{-2}$. The ratio of the breathing mode frequency for the two particles is $\frac{\omega_{\text{L}}^{2}}{\omega_{\text{S}}^{2}}\approx 2.8$ (L/S for large/small) assuming that the breathing mode frequency scales linearly with the NC diameter.Ghavanloo et al. (2015) This is similar to the ratio of the relaxation times ($650/180=3.6$) for the lower temperature. At the higher temperature, one needs to consider anharmonic effects which are more pronounced in the large NC with lower acoustic phonons. Another noticeable feature in Fig. 6 is that the Si-Si bonds seem slightly shorter in $\text{Si}_{147}\text{H}_{100}$ than in Si_{705}$$H_{300}. This results from the difference in the bond distance of atoms in the outer shell, while the inner shell atoms have similar bond distances.

IV Conclusions

In this paper we developed an ab initio Langevin dynamics approach based on a new embedded saturated fragment stochastic DFT method. We showed how the noisy forces resulting from the efsDFT calculation are used to generate a set of configurations that are distributed canonically at cryogenic and room temperatures. By proper choice of the friction coefficients and the number of stochastic orbitals, thermalization is reached within $\approx 100$ time steps for these materials, since the method is trivially parallelizable, larger computer resources we allow to easily reduce the friction coefficients thus greatly improving the sampling efficiency. While the methods presented here have already allowed impressive achievements, such as determining structural properties of silicon NCs of $3\text{nm}$ diameter containing more than $3000$ electrons, larger systems still, of unprecedented size, are now coming within our grasp due to the fact that linear-scaling highly parallelizable features of sDFT.

Acknowledgements.

E.R. acknowledges support from the Physical Chemistry of Inorganic Nanostructures Program, KC3103, Office of Basic Energy Sciences of the United States Department of Energy under Contract DE-AC02-05CH11232. D.N. acknowledge support by the NSF Grant DMR/BSF-1611382. R.B. acknowledges the US-Israel Binational Science foundation support under the BSF-NSF program, Grant 2015687.

Appendix A The embedded saturated fragments

approach

Here, we provide the technical details for the embedded fragment method described in Section II. The method corrects the stochastic estimate $\left\langle\hat{A}\right\rangle^{I}$ for the expectation value $\left\langle\hat{A}\right\rangle$ of a one-body operator $\hat{A}$, using calculations performed on $F$ separate fragments (see Eq. (11)):

[TABLE]

where the stochastic correction due to fragment $f$ is

[TABLE]

and the deterministic and the stochastic estimates, $\left\langle\hat{A}_{f}\right\rangle$ and $\left\langle\hat{A}_{f}\right\rangle^{I}$ are calculated directly on the fragment itself.

Previous implementations of embedded fragment sDFT were applied to systems of many weakly interacting molecules where the selection of fragments or clusters of such molecules was natural.Neuhauser, Baer, and Rabani (2014) We now describe a new method for defining and carrying calculations with fragments which can break up covalently bonded systems such as silicon NCs. The large system is divided into $F$ small fragments composed of one or more bonded atoms each. The surface dangling bonds of the fragment are passivated using a H atom placed in $1.46$ Å from the Si atom, in the direction of the neighboring atom which is not included in the fragment. This forms a saturated fragment. For a saturated fragment $f$, the deterministic KS-DFT method is applied to determine the KS eigenvalues $\varepsilon_{n}^{f}$ and eigenfunctions $\psi_{n}^{f}\left(\boldsymbol{r}\right)$. Further, occupation numbers $\left(p_{n}^{f}\right)^{2}=\frac{1}{2}\text{erfc}\left(\beta\left(\varepsilon_{n}^{f}-\mu_{f}\right)\right)$ are introduced for determining the saturated fragment density $n_{sf}\left(\boldsymbol{r}\right)=\sum_{n}\left(p_{n}^{f}\right)^{2}\psi_{n}^{f}\left(\boldsymbol{r}\right)^{2}$. The fragment density $n_{f}\left(\boldsymbol{r}\right)=c_{f}\left(\boldsymbol{r}\right)^{2}n_{sf}\left(\boldsymbol{r}\right)$ is “carved out” of $n_{sf}\left(\boldsymbol{r}\right)$ using a carving function $c_{f}\left(\boldsymbol{r}\right)^{2}$. Thus:

[TABLE]

where, inspired by Hirshfeld partitioning, Hirshfeld (1977) the carving function is defined as:

[TABLE]

where $n_{a}^{\left(0\right)}\left(\boldsymbol{r}\right)$ is the spherical density of neutral atom $a$. The temperature parameter $\beta$ in the definition of the population $p_{n}^{f}$ is chosen be the same value as that of the sDFT calculation, while the chemical potential $\mu_{f}$ of each fragment is determined by the condition of neutrality of the fragment:

[TABLE]

Defining non-orthogonal functions $\tilde{\psi}_{n}^{f}\left(\boldsymbol{r}\right)=c_{f}\left(\boldsymbol{r}\right)p_{n}^{f}\psi_{n}^{f}\left(\boldsymbol{r}\right)$, the fragment density of Eq. (20) becomes $n_{f}\left(\boldsymbol{r}\right)=2\sum_{n}\tilde{\psi}_{n}^{f}\left(\boldsymbol{r}\right)^{2}$, so the chemical potential is determined from the condition:

[TABLE]

After determining $\mu_{f}$ and in order to construct the reduced density matrix (RDM), we orthogonalize the functions $\tilde{\psi}_{n}^{f}\left(\boldsymbol{r}\right)$ by diagonalizing the overlap matrix $S_{nn^{\prime}}^{f}=\left\langle\tilde{\psi}_{n}^{f}\left|\tilde{\psi}_{n^{\prime}}^{f}\right.\right\rangle$, obtaining the unitary matrix $U_{f}$ of eigenvectors and the eigenvalues $s_{n}^{f}>0$ (so that $U_{f}^{T}S_{f}U_{f}=diag\left[s_{1}^{f},s_{2}^{f}\dots\right]$). The orthogonal wavefunctions are: $\phi_{m}^{f}\left(r\right)=\sum_{n}\tilde{\psi}_{n}^{f}\left(r\right)U_{nm}^{f}$ and the norm is $\left\langle\phi_{m}^{f}\left|\phi_{m}^{f}\right.\right\rangle=s_{m}^{f}$. Using the new wave functions, the unsaturated fragment density is given by:

[TABLE]

and the RDM by

[TABLE]

Using the RDM we express the unsaturated fragment expectation value appearing in Eq. (19) as:

[TABLE]

where

[TABLE]

By choosing the fragment grid-points to be a subset of the full system grid, each stochastic orbital $\chi_{i}$ ($i=1,\dots,I$) of the full system appears as a stochastic orbital on the fragment grid and can be used to perform the stochastic estimate appearing in Eq. (19) as:

[TABLE]

where the subscript $f$ on the left denotes integration over the fragment grid. The difference $\Delta A_{f}^{I}=\left\langle\hat{A}_{f}\right\rangle-\left\langle\hat{A}_{f}\right\rangle^{I}$ in Eq. (19) can now be written in a unified form as:

[TABLE]

where:

[TABLE]

Hence, by calculating the matrix $\Delta_{mm^{\prime}}^{f\,I}$ all types of expectation value corrections can be obtained from Eq. (24).

The efficacy of embedded fragments in sDFT force calculations is achieved through a reduction of the STD $\sigma\left(f_{x}\right)$ of a force component. The STD $\sigma\left(f_{x}\right)$ is proportional to $1/\sqrt{I}$, where $I$ is the number of stochastic orbitals and the proportionality constant, denoted $\sigma_{1}\left(f_{x}\right)=\sqrt{I}\sigma\left(f_{x}\right)$, is called the *inherent *STD. This quantity depends on the NC characteristics but not on the number of stochastic orbitals. In Fig. 7 we plot the inherent STD on each atom for $\text{Si}_{35}\text{H}_{36}$ and $\text{Si}_{705}\text{H}_{300}$ as a function of fragment size. Even the use of the smallest fragments reduces the inherent force STD by a significant factor, $1.6$ (for $\text{Si}_{705}\text{H}_{300}$) to $2.3$ (for $\text{Si}_{35}\text{H}_{36}$). Using larger fragments reduces the STD by an additional factor of $\approx 1.5$, with increasing effect for larger systems, since the electron density in the larger fragments is similar to that of the full system. It is interesting to see that for the forces there is no noticeable sublinear scaling: the inherent STD for both systems is similar, with the larger system having a slightly ($\approx 5\%$) STD.

In summary, the embedded fragment sDFT method serves as a way to expedite the sDFT calculation by a judicious choice of fragment size and composition. As the fragment size grows, the numerical effort invested in sDFT decreases (due to reduction of STD) while in dDFT it increases. For example, consider Fig. 7 where we showed that increasing the fragment size by a factor of $10-20$ reduces the STD by a factor of $2$ and therefore the sDFT CPU time by a factor of $\approx 2^{2}=4$. On the other hand since the fragments are ten-fold larger, the amount of dDFT work on them increases (cubically) by a factor of more than$\sim 10^{3}$. Clearly then, the optimal fragment size is system dependent. Embedded fragments have the additional benefit of providing an initial density for the SCF calculation, significantly reducing the number of SCF cycles.

Appendix B Stochastic estimates of the forces

and energy perturbations

In the stochastic method, the electronic density is (see Eq. (6)):

[TABLE]

where $\hat{P}=\sqrt{\theta_{\mu}}$ is the Chebyshev expansion of the projection operator, depending on $\beta$ and $\mu$, on the occupied space of $\hat{h}_{KS}$, and the energy is

[TABLE]

where $E_{HXC}\left[n\right]$ is the Hartree-exchange-correlation energy functional, depending only on the electronic density $n\left(\boldsymbol{r}\right)$. Under variation in position of nuclei $R$:

[TABLE]

which using

[TABLE]

can be written as:

[TABLE]

The average of the second term on the right leads to the work of the Hellmann-Feynman force Eq. (8) while the first term can be shown to vanish when a full sampling is made on $\chi$ and when $\beta\to\infty$for then $P\delta PP=0$.