An ab initio Langevin dynamics approach is developed based on stochastic density functional theory (sDFT) within a new embedded saturated 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–dissipation relation. The overall approach scales linearly with the 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 Ne = 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.

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.1–11 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 method. 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 increase 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 theory1,11 and circumvents the need for SCF iterations by the 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,12–15 which 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) methods16–19 which scale linearly (or even sublinearly) with the system size and do not rely on the sparsity of the density matrix. sDFT is a general approach to the 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,20 have been developed for the tight-binding electronic structure,21–23 for molecular electronics24 and for multi-exciton generation in nanocrsytals.25 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 nanocrystals (NCs), despite the statistical fluctuations in the force estimates. For this, we invoke the Langevin equation following the work of Attaccalite and Sorella,26 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 Ne = 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. 17. 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. 17, but also to covalently bonded systems, such as silicon NCs. The method is described in detail in  Appendix A.

Kohn-Sham density functional theory27,28 maps a system of Ne interacting electrons in an external electron-nucleus potential veNr=e4π𝜖0NZNerRN, where RN (N=1,2,) are the nuclei positions and ZNe are their charge (e is the electron charge), onto a system of non-interacting electrons (the KS system), having the same ground-state density nr. This mapping is performed by solving the KS equations27,28

h^KSϕnr=𝜀nϕnr,
(1)

where the KS Hamiltonian is

h^KS=22me2+vKSr,
(2)

and the KS potential vKSr is the sum of the external electron-nuclear potential veNr, the density-dependent Hartree potential vHr=e24π𝜖0nrrrdr, and the exchange-correlation potential vxcr,

vKSr=veNr+vHr+vxcr.
(3)

In the KS system, the density is expressed in terms of the normalized single electron KS eigenstates ϕnr and eigenvalues 𝜀n,

nr=2n𝜃μ𝜀nϕnr2,
(4)

where 𝜃x is the Heaviside function and μ is the chemical potential chosen so that 2n𝜃μ𝜀n=Ne. Equations (1)–(4) must be solved self-consistently since h^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 ONe3.

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

nr=2Tr𝜃μh^KSδrr^.
(5)

The problem now shifts into calculating self-consistently the trace in Eq. (5) [since h^KS depends on nr] rather than solving the KS equations by brute-force diagonalization. When the trace is performed using the KS eigenstates, the computational cost remains ONe3 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.20 

Using the stochastic trace formula, the density can be estimated as a symmetrized stochastic trace formula, given by16 

nIr=χ𝜃βμh^KS2δrr^,  ×𝜃βμh^KS|χχ,
(6)

where χ denotes an average over I stochastic orbitals |χ, defined as

rχ=h32eiφr
(7)

for each grid point r, the parameter h (not to be confused with the KS Hamiltonian h^KS operator) is the grid spacing, and φr are statistically independent random variables in the range [0,2π] (eiφreiφrφ=δrr). The density nr is, strictly speaking, given by the limit nr=limInIr and we approximate it with a finite I. The Heaviside function in Eq. (6) is smoothed by the function 𝜃β𝜀12erfcβ𝜀, where β is a large constant satisfying βEg1, where Eg is the KS-DFT fundamental gap. Throughout this paper we set the value of β to 100Eh1. The action of 𝜃βμh^KS on |χ is evaluated by a Chebyshev expansion in powers of the sparse KS Hamiltonian, h^KS.29 The length of the Chebyshev series is determined by the value of μEminΔE and βΔE, where ΔE=EmaxEmin2 and Emin/max are the minimal and maximal eigenvalues of 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 ONe and for certain properties even to a sub-linear scaling.16 Linear-scaling complexity is achieved due to the following facts: (1) the application of a Hamiltonian to a stochastic orbital, h^KS|χ, 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 the 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,12,30 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 (α=1,,N), which can be evaluated through the Hellmann-Feynman theorem,31,32

fα=veNrRαnIrd3r.
(8)

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.

FIG. 1.

The x-component of the atomic force statistics for the 71 atoms of Si35H36 calculated by sDFT (black dashed line for I = 16) and efsDFT (solid lines, using passivated fragments of size smaller or equal to Si5, depending on the way surface atoms are treated). For each atom α=1,,N and number of stochastic orbitals I, we present the standard deviation (STD) σfx=fαxfαx2I (left) and the mean-absolute-deviation (MAD) from the deterministic DFT (dDFT) value, fαxfαxdetI. These statistics were calculated using 60 independent efsDFT/sDFT runs. Atoms are ordered by their distance from the origin, the first 35 atoms are Si atoms followed by 36 H atoms.

FIG. 1.

The x-component of the atomic force statistics for the 71 atoms of Si35H36 calculated by sDFT (black dashed line for I = 16) and efsDFT (solid lines, using passivated fragments of size smaller or equal to Si5, depending on the way surface atoms are treated). For each atom α=1,,N and number of stochastic orbitals I, we present the standard deviation (STD) σfx=fαxfαx2I (left) and the mean-absolute-deviation (MAD) from the deterministic DFT (dDFT) value, fαxfαxdetI. These statistics were calculated using 60 independent efsDFT/sDFT runs. Atoms are ordered by their distance from the origin, the first 35 atoms are Si atoms followed by 36 H atoms.

Close modal

These sDFT forces can be expressed as

fα=fαdet+fαfluc+fαbias,
(9)

where fαdet is the deterministic (generally unknown) force, fαfluc is the pure fluctuating term, and fαbias is the bias expected to be proportional to 1I in leading order. The choice of I should be large enough to reduce fαbias to negligible values, and the only source of error in the procedure is then the statistical fluctuations proportional to 1I with vanishing mean fαfluc=0.

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. 17. 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,

nr=nFr+Δnr,
(10)

where nFr=f=1Fnfr is the density generated by the individual fragments obtained from a deterministic KS-DFT (dDFT) calculation for each fragment and Δnr=nIrnFIr is a correction term evaluated using stochastic orbitals. Here, nIr is given by Eq. (6) and nFIr=f=1FnfIr is a sum over a stochastic estimate of the fragments density. In the limit I, 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, nFr, provides a reasonable approximation for the full density nr, the correction term, Δnr, 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:

nr=nIr+f=1FnfrnfIr,
(11)

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 (MADs) with respect the deterministic DFT of the atomic forces in an Si35H36 NC using hydrogen passivated Si5 fragments. The results are shown in Fig. 1.33,43–45 The STDs and MADs decrease as 1I, indicating that the bias in the force estimation is negligible. The standard deviations in the sDFT forces are larger by a factor of 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).

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 trajectory34–37 is a sequence of configurations p,qm=ptm,qtm at discrete “times” tm=mΔt, where Δt is the time step, and qq1,,qN and pp1,,pN are the Cartesian coordinates and conjugate momenta, respectively, for the N atoms. The trajectory is a solution of the Langevin equation (LE) of motion,38 

μαq̈α=fαqγαpα+ηα,
(12)

where μα is the mass of the atom α, γα is its friction constant, and fα=fαdet+fαfluc is the total efsDFT force acting on it, including deterministic and fluctuating parts [see Eq. (9)]. The bias is assumed negligible, so that fα=fαdet. In Eq. (12)ηα 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 obeys

ηαt=fαfluct=0

and

ηαt+fαfluctηαt+fαfluct=ηαtηαt+fαfluctfαfluct=I3×3σα2δααδtt,
(13)

where α,α=1,N are atom indices, designates average over the atomic force distribution, I3×3 is the 3×3 unit matrix, and σα is the atomic force STD of atom α, which is taken to satisfy the fluctuation-dissipation relation,

σα2=2μαγαkBT.
(14)

We use the Verlet-like algorithm39 for numerically integrating the LE of motion at a fixed temperature T and a predefined time step Δ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 ηαm is sampled from a Gaussian distribution such that the discretized version of Eq. (13) holds: ηαm+fαflucηαn+fαflucΔt=Iσα2δααδmn,

qαm+1=qαm+bαΔtμα1pαm+12bαΔt2μα1fαm+ηαm+1,pαm+1=aαpαm+12Δtaαfαm+fαm+1+2bαηαm+1,
(15)

where aα=bα112γαΔt and bα1=1+12γαΔ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 ηα of the force differently from the force fα=fαdet+fαfluc that result from the efsDFT calculation containing deterministic and fluctuating components that cannot be separated.

In Fig. 2, we plot for Si35H36 the running average of the transient temperature, Tm, calculated from the kinetic energy

TKm=23NkBαpαm22μα
(16)

and from the virial estimator,

TVm=13NkBαfαm+ηαmqαmqα.
(17)

In the above, qα is the time average of the coordinate of atom α. 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 = 300 K. This non-equilibrium initial configuration relaxes towards equilibrium.

FIG. 2.

The Verlet-like39 temperature in Si35H36 evaluated with efsDFT using transient kinetic energy (dotted line), running average kinetic energy for Si (red), and the running average virial (black) using I = 30 stochastic orbitals, a time step of Δt=1.2fs, and the friction coefficients γSi=γH=0.04fs1. In the canonical distribution, the average kinetic energy T is equal to 3N2δT, where δT is the fluctuation in the kinetic energy.40 The standard deviation of the transient shown as a dotted line is δT=29 K and multiplied by 3N2, where N = 71 is the number of atoms in the system, gives 299.3 K, which is close to the designated temperature. The average kinetic energy of Si and H is 315 K and 285 K, respectively.

FIG. 2.

The Verlet-like39 temperature in Si35H36 evaluated with efsDFT using transient kinetic energy (dotted line), running average kinetic energy for Si (red), and the running average virial (black) using I = 30 stochastic orbitals, a time step of Δt=1.2fs, and the friction coefficients γSi=γH=0.04fs1. In the canonical distribution, the average kinetic energy T is equal to 3N2δT, where δT is the fluctuation in the kinetic energy.40 The standard deviation of the transient shown as a dotted line is δT=29 K and multiplied by 3N2, where N = 71 is the number of atoms in the system, gives 299.3 K, which is close to the designated temperature. The average kinetic energy of Si and H is 315 K and 285 K, respectively.

Close modal

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 ±5% deviations in the temperatures. These deviations are slightly inconsistent with equipartition theorem and may result from several factors, such as the finite time step of the Langevin propagator,39 incomplete SCF convergence,40 and insufficiently accurate estimate of the amount of white noise ηαm in Eq. (13) required to fulfill the fluctuation-dissipation relation.

The effect of γSi on the configurational relaxation and on the velocity autocorrelation decay is illustrated in Fig. 3 for Si36H35. In order to decrease the number of unknown parameters, we set the values of γα and σα 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 σ, we introduced white noise ηαi=σαi2fαi,fluc2 for each degree of freedom [see Eq. (13)], where fαi,fluc2 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 ffluc2 is not sensitive to the particular configuration used.

FIG. 3.

The normalized velocity autocorrelation function (top panel) and the mean nearest-neighbor Si–Si distance (bottom panel) in Si35H36 as a function of time for a LD trajectory at T = 300 K with time step Δt=1.2 fs calculated using a dDFT based LD for different values of γ=γSi=γH. The dashed curve corresponds to efsDFT based LD calculation with γ=0.04fs1 and I = 30 stochastic orbitals. The simulation started intentionally from an inflated configuration in order to measure the relaxation time.

FIG. 3.

The normalized velocity autocorrelation function (top panel) and the mean nearest-neighbor Si–Si distance (bottom panel) in Si35H36 as a function of time for a LD trajectory at T = 300 K with time step Δt=1.2 fs calculated using a dDFT based LD for different values of γ=γSi=γH. The dashed curve corresponds to efsDFT based LD calculation with γ=0.04fs1 and I = 30 stochastic orbitals. The simulation started intentionally from an inflated configuration in order to measure the relaxation time.

Close modal

As expected, the configurational relaxation time increases with increasing values of γSi with the opposite trend for the decay time of the velocity autocorrelation function. Based on the results of Si36H35 presented in Fig. 3, we conclude that a friction coefficient of 0.04fs1 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, γSi=0.04fs1 for Silicon and γH=0.12fs1 for the lighter H atoms. These values were used for the larger systems described in Sec. III (see Table I). 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 γSi=0.04fs1 (dotted red line) proving that the relaxation times are similar to those of the dDFT based LD calculation with the same value of γSi.

TABLE I.

Value of various parameters for the LD based on efsDFT calculations: The friction coefficients γ, number of stochastic orbitals I, time-step Δt, and the wall time per single SCF iteration titer.

γfs1
NCT(K)HSiIΔtfs1titermin
Si35H36 30 0.12 0.04 120 1.2 
 300 0.04 0.04 30 1.2 
Si147H100 30 0.12 0.04 120 1.2 
 300 0.12 0.04 30 1.2 
Si705H300 30 0.12 0.04 120 1.2 10 
 300 0.12 0.04 92 1.2 10 
γfs1
NCT(K)HSiIΔtfs1titermin
Si35H36 30 0.12 0.04 120 1.2 
 300 0.04 0.04 30 1.2 
Si147H100 30 0.12 0.04 120 1.2 
 300 0.12 0.04 30 1.2 
Si705H300 30 0.12 0.04 120 1.2 10 
 300 0.12 0.04 92 1.2 10 

Validation of the structure obtained using efsDFT based LD is demonstrated using the pair distribution function gr.34 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. Figure 4 shows a close agreement between the dDFT and efsDFT based LD estimates of gr of the Si35H36 NC at two temperatures. The inset focuses on the first peak in gr, comparing the efsDFT to dDFT at T = 30 and 300K.

FIG. 4.

The Si–Si pair distribution function gr for Si35H36 calculated using dDFT (dotted lines) and efsDFT based LD (solid lines, see Table I for parameters) at T=30K (blue curves) and 300K (red curves). Inset: Details of the first (nearest neighbor) peak. The dotted lines are for different friction coefficients γ=γSi=γH in the range 0.020.5fs1.

FIG. 4.

The Si–Si pair distribution function gr for Si35H36 calculated using dDFT (dotted lines) and efsDFT based LD (solid lines, see Table I for parameters) at T=30K (blue curves) and 300K (red curves). Inset: Details of the first (nearest neighbor) peak. The dotted lines are for different friction coefficients γ=γSi=γH in the range 0.020.5fs1.

Close modal

In Secs. I and II, we presented the methods and assessed the accuracy and validity of the efsDFT based LD. Here we apply the method to study the structural properties of larger NCs exceeding Ne = 3000 electrons. The Si–Si pair distribution functions gr at two temperatures T (30 and 300K) are displayed in the upper panel of Fig. 5 for Si147H100 and Si705H300. Temperature broadens the peaks by a factor of 2–3 without significantly changing the peak position.

FIG. 5.

The Si–Si pair distribution function gr for Si147H100 (left) and Si705H300 (right) calculated using efsDFT based LD. Upper panels: gr for T=30K (blue curves) and T=300K (red curves). Lower panels: The first peak of gr at 30K shown for several times. The calculation parameters are given in Table I.

FIG. 5.

The Si–Si pair distribution function gr for Si147H100 (left) and Si705H300 (right) calculated using efsDFT based LD. Upper panels: gr for T=30K (blue curves) and T=300K (red curves). Lower panels: The first peak of gr at 30K shown for several times. The calculation parameters are given in Table I.

Close modal

In the lower panel of Fig. 5, we plot the transient and relaxed gr at 30K for the two systems, focusing on the first, nearest neighbor peak. As described also for Si35H36, 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 was 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 30K, the relaxation times are 180 and 650 fs for Si147H100 and Si705H300, respectively. For 300K, 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 Si147H100 (spherical shells A-B) and Si705H300 (shells A-E); see Table II for the definition and properties of the shells. In Si705H300 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, ω, of the breathing mode of the larger NC. In the limit of an overdamped motion (as is the case here since γ2ω2), the relaxation is dominated by two time scales proportional to γ1 and (ω2γ)1. The former leads to a fast relaxation, while the latter is slower and depends on the value of ω2. The ratio of the breathing mode frequency for the two particles is ωL2ωS22.8 (L/S for large/small) assuming that the breathing mode frequency scales linearly with the NC diameter.41 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 Si147H100 than in Si705H300. This results from the difference in the bond distance of atoms in the outer shell, while the inner shell atoms have similar bond distances.

FIG. 6.

Si–Si nearest neighbor distance averaged over atoms in shells A-E (see Table II for definition), at 30K for Si147H100 (left) and Si705H300 (right) as a function of time. The calculation parameters are given in Table I.

FIG. 6.

Si–Si nearest neighbor distance averaged over atoms in shells A-E (see Table II for definition), at 30K for Si147H100 (left) and Si705H300 (right) as a function of time. The calculation parameters are given in Table I.

Close modal
TABLE II.

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 NSi, and the number of nearest neighbor (NN) Si–Si pairs NNN.

ShellRinRoutNSiNNN
A 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 
ShellRinRoutNSiNNN
A 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 

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. Choosing the friction coefficients and the number of stochastic orbitals according to the criteria presented in Sec. II D, thermalization is reached within 100 time steps for these materials, since the method is trivially parallelizable, larger computer resources would 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 the structural properties of silicon NCs of 3nm diameter containing more than 3000 electrons, larger systems still, of unprecedented size, are now coming within our grasp due to the linear-scaling highly parallelizable features of sDFT.

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 No. DE-AC02-05CH11232. D.N. acknowledges support from the NSF Grant No. DMR/BSF-1611382. R.B. acknowledges support from the US-Israel Binational Science foundation under the BSF-NSF program, Grant No. 2015687.

Here, we provide the technical details for the embedded fragment method described in Sec. II. The method corrects the stochastic estimate ÂI for the expectation value  of a one-body operator Â, using calculations performed on F separate fragments [see Eq. (11)],

Â=ÂI+f=1FΔAfI,
(A1)

where the stochastic correction due to fragment f is

ΔAfI=A^fA^fI,
(A2)

and the deterministic and the stochastic estimates, A^f and A^fI, 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.17 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 an 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 𝜀nf and eigenfunctions ψnfr. Further, occupation numbers pnf2=12erfcβ𝜀nfμf are introduced for determining the saturated fragment density nsfr=npnf2ψnfr2. The fragment density nfr=cfr2nsfr is “carved out” of nsfr using a carving function cfr2. Thus

nfr=cfr2npnf2ψnfr2,
(A3)

where, inspired by Hirshfeld partitioning,42 the carving function is defined as

cf(r)=afna0rasfna0r,

where na0r is the spherical density of neutral atom a. The temperature parameter β in the definition of the population pnf is chosen to be the same value as that of the sDFT calculation, while the chemical potential μf of each fragment is determined by the condition of neutrality of the fragment,

nfrdr=afna0rdr.
(A4)

Defining non-orthogonal functions ψ̃nfr=cfrpnfψnfr, the fragment density of Eq. (A3) becomes nfr=2nψ̃nfr2, so the chemical potential is determined from the condition

2nψ̃nfψ̃nf=afna0rdr.
(A5)

After determining μf and in order to construct the reduced density matrix (RDM), we orthogonalize the functions ψ̃nfr by diagonalizing the overlap matrix Snnf=ψ̃nfψ̃nf, obtaining the unitary matrix Uf of eigenvectors and the eigenvalues snf>0 (so that UfTSfUf=diags1f,s2f). The orthogonal wavefunctions are ϕmfr=nψ̃nfrUnmf and the norm is ϕmfϕmf=smf. Using the new wave functions, the unsaturated fragment density is given by

nfr=2mϕmfr2

and the RDM by

𝜃^f=2mϕmfϕmf.

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

A^ftr𝜃^fÂ=tr𝜃^fÂ𝜃^f,

where

𝜃^f=2msmf12ϕmfϕmf.
(A6)

By choosing the fragment grid-points to be a subset of the full system grid, each stochastic orbital χi (i=1,,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. (A2) as

A^fI=1Iiχi𝜃̂fÂ𝜃̂fχif,

where the subscript f on the left denotes integration over the fragment grid. The difference ΔAfI=A^fA^fI in Eq. (A2) can now be written in a unified form as

ΔAfI=2mmΔmmfIϕmfÂϕmff,
(A7)

where

ΔmmfIδmm1Iiχiϕmffϕmfχifsmfsmf.
(A8)

Hence, by calculating the matrix ΔmmfI all types of expectation value corrections can be obtained from Eq. (A7).

The efficacy of embedded fragments in sDFT force calculations is achieved through a reduction of the STD σfx of a force component. The STD σfx is proportional to 1I, where I is the number of stochastic orbitals and the proportionality constant, denoted by σ1fx=Iσfx, 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 Si35H36 and Si705H300 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 Si705H300) to 2.3 (for Si35H36). Using larger fragments reduces the STD by an additional factor of 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 slightly higher (5%) STD.

FIG. 7.

Efficacy of fragments on the inherent sDFT STD σ1fx of the x-component of the force on each atom of Si35H36 (left) and Si705H300 (right) NCs. The inherent STD σ1 is the actual STD σ times I. Si atoms are shown first followed by H atoms, where atoms are ordered by distance from the NC center. Calculations are done on the configuration cut out from the bulk silicon, where H atoms were placed near the surface for saturating the dangling bonds.

FIG. 7.

Efficacy of fragments on the inherent sDFT STD σ1fx of the x-component of the force on each atom of Si35H36 (left) and Si705H300 (right) NCs. The inherent STD σ1 is the actual STD σ times I. Si atoms are shown first followed by H atoms, where atoms are ordered by distance from the NC center. Calculations are done on the configuration cut out from the bulk silicon, where H atoms were placed near the surface for saturating the dangling bonds.

Close modal

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 22=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 103. 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.

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

nr=2χP^δrr^P^χ,
(B1)

where P^=𝜃μ is the Chebyshev expansion of the projection operator, depending on β and μ, on the occupied space of h^KS, and the energy is

E=2χP^T^P^χ+veNr;Rnrdr+EHXCn,=2χP^T^+veNr;RP^χ+EHXCn,
(B2)

where EHXCn is the Hartree-exchange-correlation energy functional, depending only on the electronic density nr. Under variation in position of nuclei R,

δE=2χP^T^+vr^,RδP^χ+2χδP^T^+vr^,RP^χ+χP^δvr^,RP^χ+vHXCrδnrdr,
(B3)

which using

δnr=2χδP^δrr^P^χ+2χP^δrr^δP^χ
(B4)

can be written as

δE=2χh^KSP^δP^+δP^P^h^KSχ+χP^δvr^,RP^χ,=2χh^KSP^δP^+δP^P^h^KSχ+nrδvr;Rd3r
(B5)

The average of the second term on the right leads to the Hellmann-Feynman force Eq. (8), while the first term can be shown to vanish when a full sampling is made on χ and when β, owing to the fact that within these limits PδPP=0.

1.
R.
Car
and
M.
Parrinello
,
Phys. Rev. Lett.
55
,
2471
(
1985
).
2.
R. N.
Barnett
,
U.
Landman
,
A.
Nitzan
, and
G.
Rajagopal
,
J. Chem. Phys.
94
,
608
(
1991
).
3.
M. C.
Payne
,
M. P.
Teter
,
D. C.
Allan
,
T. A.
Arias
, and
J. D.
Joannopoulos
,
Rev. Mod. Phys.
64
,
1045
(
1992
).
4.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
5.
M.
Tuckerman
,
K.
Laasonen
,
M.
Sprik
, and
M.
Parrinello
,
J. Chem. Phys.
103
,
150
(
1995
).
6.
D.
Marx
and
J.
Hutter
, “
Ab initio molecular dynamics: Theory and implementation
,” in
Modern Methods and Algorithms of Quantum Chemistry, Proceedings
, NIC Series Vol. 3, edited by
J.
Grotendorst
(
John von Neumann Institute for Computing
,
Julich
,
2000
), p.
329
.
7.
H. B.
Schlegel
,
J. M.
Millam
,
S. S.
Iyengar
,
G. A.
Voth
,
A. D.
Daniels
,
G. E.
Scuseria
, and
M. J.
Frisch
,
J. Chem. Phys.
114
,
9758
(
2001
).
8.
J. M.
Herbert
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
7
,
3269
(
2005
).
9.
J.-Y.
Raty
,
F.
Gygi
, and
G.
Galli
,
Phys. Rev. Lett.
95
,
096103
(
2005
).
10.
M.
Bockstedte
,
A.
Kley
,
J.
Neugebauer
, and
M.
Scheffler
,
Comput. Phys. Commun.
107
,
187
(
1997
).
11.
M.
Cawkwell
,
A. M.
Niklasson
, and
D. M.
Dattelbaum
,
J. Chem. Phys.
142
,
064512
(
2015
).
12.
S.
Goedecker
,
Rev. Mod. Phys.
71
,
1085
(
1999
).
13.
M. J.
Gillan
,
D. R.
Bowler
,
A. S.
Torralba
, and
T.
Miyazaki
,
Comput. Phys. Commun.
177
,
14
(
2007
).
14.
R.
Baer
and
M.
Head-Gordon
,
Phys. Rev. Lett.
79
,
3962
(
1997
).
15.
J. M.
Soler
,
E.
Artacho
,
J. D.
Gale
,
A.
Garcia
,
J.
Junquera
,
P.
Ordejon
, and
D.
Sanchez-Portal
,
J. Phys.: Condens. Matter
14
,
2745
(
2002
).
16.
R.
Baer
,
D.
Neuhauser
, and
E.
Rabani
,
Phys. Rev. Lett.
111
,
106402
(
2013
).
17.
D.
Neuhauser
,
R.
Baer
, and
E.
Rabani
,
J. Chem. Phys.
141
,
041102
(
2014
).
18.
Y.
Gao
,
D.
Neuhauser
,
R.
Baer
, and
E.
Rabani
,
J. Chem. Phys.
142
,
034106
(
2015
).
19.
D.
Neuhauser
,
E.
Rabani
,
Y.
Cytter
, and
R.
Baer
,
J. Phys. Chem. A
120
,
3071
(
2015
).
20.
M. F.
Hutchinson
,
Commun. Stat. Simul. Comput.
19
,
433
(
1990
).
21.
D. A.
Drabold
and
O. F.
Sankey
,
Phys. Rev. Lett.
70
,
3631
(
1993
).
22.
L. W.
Wang
,
Phys. Rev. B
49
,
10154
(
1994
).
23.
H.
Röder
,
R.
Silver
,
D.
Drabold
, and
J. J.
Dong
,
Phys. Rev. B
55
,
15382
(
1997
).
24.
R.
Baer
,
T.
Seideman
,
S.
Ilani
, and
D.
Neuhauser
,
J. Chem. Phys.
120
,
3387
(
2004
).
25.
R.
Baer
and
E.
Rabani
,
Nano Lett.
12
,
2123
(
2012
).
26.
C.
Attaccalite
and
S.
Sorella
,
Phys. Rev. Lett.
100
,
114501
(
2008
).
27.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
28.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
29.
H.
Tal-Ezer
and
R.
Kosloff
,
J. Chem. Phys.
81
,
3967
(
1984
).
30.
R.
Baer
and
M.
Head-Gordon
,
J. Chem. Phys.
107
,
10003
(
1997
).
31.
J.
Hellman
,
Einfuhrung in die Quantenchemie
(
Springer
,
Deuticke, Leipzig
,
1937
).
32.
R. P.
Feynman
,
Phys. Rev.
56
,
340
(
1939
).
33.

All calculations in this work use real-space grids of spacing Δx=0.5a0 and Troullier-Martins norm-conserving pseudopotentials43 within the Kleinman-Bylander approximation.44 Fast Fourier transforms were used for applying the kinetic energy operator and for determining the Hartree potentials, and the method of Ref. 45 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).

34.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
Oxford
,
1987
), p.
xix, 385
.
35.
N. G.
Van Kampen
,
Stochastic Processes in Physics and Chemistry
(
Elsevier
,
1992
), Vol. 1.
36.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
USA
,
2001
).
37.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
, 2nd ed. (
Academic Press
,
San Diego
,
2002
), p.
xxii, 638
.
38.
P.
Langevin
,
C. R. Acad. Sci. Paris
146
,
530
(
1908
).
39.
N.
Grønbech-Jensen
and
O.
Farago
,
Mol. Phys.
111
,
983
(
2013
).
40.
E.
Martínez
,
M. J.
Cawkwell
,
A. F.
Voter
, and
A. M.
Niklasson
,
J. Chem. Phys.
142
,
154120
(
2015
).
41.
E.
Ghavanloo
,
S.
Fazelzadeh
,
T.
Murmu
, and
S.
Adhikari
,
Phys. E (Amsterdam, Neth.)
66
,
228
(
2015
).
42.
F.
Hirshfeld
,
Theor. Chim. Acta
44
,
129
(
1977
).
43.
N.
Troullier
and
J. L.
Martins
,
Phys. Rev. B
43
,
1993
(
1991
).
44.
L.
Kleinman
and
D.
Bylander
,
Phys. Rev. Lett.
48
,
1425
(
1982
).
45.
G. J.
Martyna
and
M. E.
Tuckerman
,
J. Chem. Phys.
110
,
2810
(
1999
).