Debugging:

Your Gaussian functions should contain one unit of charge. To debug GaussR, add the following code fragment (with appropriate initializations) to main() in physics.C

  rho1.evalr(GaussR,1.);
  rho2.evalr(GaussR,2.);
  rho=rho1-rho2;
  printf("Integral of rho1: (%f,%f)\n",sum(rho1)*V/N);
  printf("Integral of rho2: (%f,%f)\n",sum(rho2)*V/N);
  printf("Integral of rho: (%f,%f)\n",sum(rho)*V/N);

With an FFT box of size n=32 and grid spacing of d=1 (bohrs), you should get the correct total charges to better than 0.000001 (electron charges).