getn(), getE(), H()

function n=getn(psi,f)

Input:

• psi: expansion coefficients of $N_s$ orthonormal wave functions
• f: fillings of the orbitals (typically f=2).

Output:

• n: sample values of electron density, $n_i=f \sum_k \psi_{ik}^* \psi_{ik}$

Both getE() and H() compute the electron density with a call to I() on the wave functions. Rather than updating two nearly identical blocks of code, produce a function of the above prototype which calculates the density with an explicit loop over the columns of psi rather than using calls to diagouter(). Note that you should not need to use any global variables for this function.

Finally, replace the current expressions which calculate the density n in getE() and H() with calls to getn(). Note that, depending on your implementation, you may have to modify getE() and H() somewhat to first obtain the orthonormal wave functions from the input W. The following statement will do this for you,

Y=W*inv(sqrtm(W'*O(W))); %# Orthonormal wave functions

Debugging: Rerun ``Hatoms.m'' to convince yourself that the output is identical. (You need not rerun the entire calculation; because we always seed the random number generator with the same number, identical results through the end of the finite-difference test should be sufficient.)