Global density of states method

The recursion method just described is an optimal method for generating a density of states, and hence the band energy, from the moments of the Hamiltonian. For molecular dynamics simulations, however, we also need the atomic forces. These have a contribution from the band energy, which is given by the derivative of the band energy with respect to atomic positions. Combining equations 6, 20 and 24, we get:

$$U_{band}=-{1\over{\pi}}\sum_{i\alpha}\lim_{\eta\to 0} {\rm Im}\int{\rm d}EG_{i\alpha,i\alpha}(E+{\rm i}\eta)Ef(x).$$ (31)

The contribution to the force on atom k from the band energy (at constant chemical potential) is then:

 

$$-{{\partial U_{band}}\over{\partial\vec r_k}}$$ F^(k)_band =

$$-{1\over{\pi}}\sum_{i\alpha}\lim_{\eta\to 0}{\rm Im}\int{\rm d}E {{\partial G_{i\alpha,i\alpha}(E+{\rm i}\eta)}\over{\partial\vec r_k}}Ef(x).$$ = (32)

Thus we see that the force depends on the derivative of the Green's function. Using the recursion expression for G₀₀(Z) (equation 27) with $\mid U_0\rangle = \mid i\alpha \rangle$, and the fact that the recursion coefficients depend on the moments (equation 30) we can use the chain rule for partial differentiation to give:

 

$${{\partial G_{00}(Z)}\over{\partial\vec r_k}} \sum_n\sum_{p=1}^{2n+1} {{\partial G_{00}(Z)}\over{\partial a_n}}... ...ial\mu^{(p)}_{i\alpha}}} {{\partial\mu^{(p)}_{i\alpha}}\over{\partial\vec r_k}}$$ =

$$\sum_n\sum_{p=1}^{2n} {{\partial G_{00}(Z)}\over{\partial b_n}} {... ...al\mu^{(p)}_{i\alpha}}} {{\partial\mu^{(p)}_{i\alpha}}\over{\partial\vec r_k}}.$$ + (33)

The derivatives of the Green's function with respect to the recursion coefficients can be evaluated straighforwardly ( ${{\partial G_{00}(Z)}/{\partial a_n}}=G_{0n}(Z)G_{n0}(Z)$ and ${{\partial G_{00}(Z)}/{\partial b_n}}=2G_{0(n-1)}(Z)G_{n0}(Z)$), as can the derivatives of the recursion coefficients with respect to the moments[10]. The derivatives of the moments with respect to the atomic positions, however, are more problematic. Formally we can write down the derivative of the p^th moment:

 

$${{\partial\mu^{(p)}_{i\alpha}}\over{\partial\vec r_k}} {{\partial}\over{\partial\vec r_k}}\sum_{j_1\beta_1...j_{p-1}\bet... ..._{i\alpha,j_1\beta_1}H_{j_1\beta_1,j_2\beta_2}...H_{j_{p-1}\beta_{p-1},i\alpha}$$ =

$$\sum_{j_1\beta_1...j_{p-1}\beta_{p-1}}\left\{ {{\partial H_{i\alp... ...\vec r_k}} H_{j_1\beta_1,j_2\beta_2}\dots H_{j_{p-1}\beta_{p-1},i\alpha}\right.$$ =

$$\left. H_{i\alpha,j_1\beta_1}{{\partial H_{j_1\beta_1,j_2\beta_2}} \over{\partial\vec r_k}} \dots H_{j_{p-1}\beta_{p-1},i\alpha}+\dots\right\}.$$ + (34)

Unfortunately this expression is, in general, very slow to evaluate on a computer (though it can be used for very low order moment expansions[32]). The reason for this is that the n^th moment needs to be evaluated n times for each component of force, and there will be 3N components to be considered, where N is the number of atoms in the cluster from which the moment is evaluated.

However, there is a way to greatly accelerate the evaluation of the derivatives, and that is to work with the global moments, rather than the local moments. The global moments ( $\bar\mu^{(n)}$) are defined by:

 

$$\bar\mu^{(p)} \sum_{i\alpha,j_1\beta_1...j_{p-1}\beta_{p-1}} H_{i\alpha,j_1\beta_1}H_{j_1\beta_1,j_2\beta_2}...H_{j_{p-1}\beta_{p-1},i\alpha}$$ =

$$\sum_{i\alpha}\langle i\alpha\mid\hat H^p\mid i\alpha\rangle$$ =

$${\rm Tr}\left\{\hat H^p\right\}.$$ = (35)

Because we can permute matrices inside a trace, the derivative of the global moment is given by:

$${{\partial\bar\mu^{(p)}}\over{\partial\vec r_k}}=p\sum_{i\alp... ...}} H_{j_1\beta_1,j_2\beta_2}...H_{j_{p-1}\beta_{p-1},i\alpha}.$$ (36)

This is very easy to calculate efficiently on a computer[10], though presumably it can only be evaluated in a stable manner for about the first 20 moments. In the global density of states method, the global moments are used to construct recursion coefficients from which the density of states, band energy and atomic forces are evaluated.

In summary, the use of the global density of states leads to a reduced rate of convergence of the energy with number of moments as compared with the local densities of states for inhomogeneous systems. Further, only about 21 moments can be used before the conversion of moments into recursion coefficients becomes unstable, and since all the moments are stored, this method can be quite memory intensive. This has to be weighed against the chief benefit of using the global density of states, namely that the analytical forces are exact derivatives of the energy. This method can be implemented on a parallel machine with almost perfect scaling, as little information needs to be passed between processors.