2
$\begingroup$

I'm computing the Hartree potentials of atoms by solving the Poisson equation and I use hydrogen atom as a test case. The Poisson equation for hydrogen atom in atomic units is given by $$\nabla^2 V_H = -4 \exp(-2 r)$$ where $r = \sqrt{x^2+y^2+z^2}$. The numerical solution of $V_H(x,y,z)$ with $z=0$ is illustrated in the following figure: Graph of the Hartree potential

The numerical solution is computed with the Conjugate Gradient method so that the Laplacian is computed with stencil $$ \nabla^2 f(x,y,z) \approx \frac{a+b+c-6 f(x,y,z)}{h^2} $$ where $$a = f(x+h,y,z)+f(x-h,y,z),$$ $$b = f(x,y+h,z)+f(x,y-h,z),$$ and $$c = f(x,y,z+h)+f(x,y,z-h).$$

As the right-hand side of the Poisson equation is spherically symmetric the Poisson equation takes the form $$\frac{d^2 v_H}{dr^2} + \frac{2}{r} \frac{d v_H}{dr} = -4 \exp(-2 r)$$ where $v_H$ is the Hartree potential as a function of $r$. The solution of this equation is $$v_H(r) = - \frac{r + 1}{r} \exp(-2r),$$ which is illustrated in the following figure: Graph of the solution of the differential equation (here $r=\sqrt{x^2+y^2+z^2}$ and $z=0$).

The Octave code computing the Laplacian is here:

function aLap = StencilLap3d( a, rStep )
  rH2 = rStep * rStep;
  vSize = size( a );
  nXSize = vSize( 1 );
  nYSize = vSize( 2 );
  nZSize = vSize( 3 );
  nXDim = ( nXSize - 1 ) / 2;
  nYDim = ( nYSize - 1 ) / 2;
  nZDim = ( nZSize - 1 ) / 2;
  assert( nXDim == round( nXDim ) );
  assert( nYDim == round( nYDim ) );
  assert( nZDim == round( nZDim ) );
  aLap = zeros( nXSize, nYSize, nZSize );
  rXP = 0;
  rXM = 0;
  rYP = 0;
  rYM = 0;
  rZP = 0;
  rZM = 0;
  for nX = (-nXDim):nXDim
    for nY = (-nYDim):nYDim
      for nZ = (-nZDim):nZDim
        if ( nX < nXDim )
          rXP = a( nXDim + 1 + nX + 1, nYDim + 1 + nY, nZDim + 1 + nZ );
        else
          rXP = 0.0;
        endif
        if ( nX > -nXDim )
          rXM = a( nXDim + 1 + nX - 1, nYDim + 1 + nY, nZDim + 1 + nZ );
        else
          rXM = 0.0;
        endif
        if ( nY < nYDim )
          rYP = a( nXDim + 1 + nX, nYDim + 1 + nY + 1, nZDim + 1 + nZ );
        else
          rYP = 0.0;
        endif
        if ( nY > -nYDim )
          rYM = a( nXDim + 1 + nX, nYDim + 1 + nY - 1, nZDim + 1 + nZ );
        else
          rYM = 0.0;
        endif
        if ( nZ < nZDim )
          rZP = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ + 1 );
        else
          rZP = 0.0;
        endif
        if ( nZ > -nZDim )
          rZM = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ - 1 );
        else
          rZM = 0.0;
        endif
        n0 = a( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ );
        aLap( nXDim + 1 + nX, nYDim + 1 + nY, nZDim + 1 + nZ ) = ...
          ( rXP + rXM + rYP + rYM + rZP + rZM - 6 * n0 ) / rH2;
      endfor
    endfor
  endfor
endfunction

Obviously the solution of the differential equation should be the same as the numerical solution, so there is something wrong here. Based on my previous computations I suspect that the error is in the differential equation. Can somebody tell what is wrong?

$\endgroup$
  • $\begingroup$ May I ask which tool are you using for the plots? $\endgroup$ – Algo Aug 13 at 14:52
  • $\begingroup$ I use Octave. IMHO it is very good software. It contains a symbolic math package, too. $\endgroup$ – Tommi Höynälänmaa Aug 13 at 15:21
  • $\begingroup$ The sign of the numerical and exact solutions are opposite. I would suspect a bug in your code. $\endgroup$ – cfdlab Aug 14 at 5:53
  • 1
    $\begingroup$ You are solving the Poisson eqn in $R^3$ and the exact solution is in $L^2(R^3)$ but not in $H^1(R^3)$. $\endgroup$ – cfdlab Aug 14 at 6:19
  • $\begingroup$ I have checked the numerical result by substituting it in the equation $\nabla^2 V_H = -4 \exp(-2r)$. What is the significance of the exact solution not being in $H^1(R^3)$? $\endgroup$ – Tommi Höynälänmaa Aug 14 at 8:56
1
$\begingroup$

Problem solved: I just have to include the solution of the homogeneous equation into the solution of the differential equation.

The solution to the differential equation is $$V_H(r) = -\frac{K}{r} - \frac{r+1}{r} \exp(-2r).$$ By requiring that $\lim_{r \to 0} V_H(r)$ is real we get $K=-1$ and $$V_H(r) = \frac{1}{r} - \frac{r+1}{r} \exp(-2r).$$

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Could you expand this answer by giving the complete exact solution. $\endgroup$ – cfdlab Aug 14 at 13:15
  • $\begingroup$ Thanks, it is clear now. You can accept your own answer. $\endgroup$ – cfdlab Aug 15 at 5:03
  • $\begingroup$ What do you mean by accepting the answer? $\endgroup$ – Tommi Höynälänmaa Aug 15 at 10:13
  • $\begingroup$ There is a tick mark beside your answer, hovering on it says "Accept this answer if it solved ...". You click that to accept your own answer and close this post. $\endgroup$ – cfdlab Aug 15 at 13:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.