To verify your implementation of the Kohn-Sham equations, in the next section you will check with finite differences that your code for the gradient of the energy function is correct.
Not infrequently, when implementing a complicated test like this for debugging purposes, one spends much time looking for a bug in perfectly working code because of a bug in the test itself! Thus, it is good programming practice to fist debug the test by running it with code which you know to be working.
To do this for our finite difference test, implement the following code fragment. It is important that the test be put in after u is initialized to random values but before it it normalized so that the gradient formula is tested for non-normalized vectors (the usual case).
u.eval(rnd); // Initialize u to random elements
.
.
// Finite diff test on gradient:
dir.eval(rnd); // Random elements
{
double oE,E;
ograd=get_grad_and_E(oE,u,Vtilde);
for (double alpha=1; alpha>=1e-8; alpha/=10) {
grad=get_grad_and_E(E,u+alpha*dir,Vtilde);
printf("Gradient test: %20.16f\n",(E-oE)/real(dir^ograd)/alpha);
printf(" (precision %20.16f)\n",fabs((1e-15*E)/(E-oE)));
}
}
.
.
u=u*(1./sqrt(u^O(u))); // Normalize u for starting
.
.
// Proceed to rest of preconditioned conjugate-gradient minimization code...
When you run this, you should find (just as you did with the test of
excp earlier in the term) that the ``Gradient test'' approaches a
value of 1, adding one digit for each new iteration. This pattern
should continue for a number of iterations until rounding errors
become significant. Note that, in terms of precision, you can only
trust those digits in the gradient test which are zero
on the ``precision'' line.
Hint: You will probably find that the formula for the gradient
given on HW#10 was missing a factor of 2. This factor canceled out
of all of our calculations on that problem set, but you should now put
it in the gradient formula so that the gradient test actually does approach 1.