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$
7
  • $\begingroup$ May I ask which tool are you using for the plots? $\endgroup$
    – Algo
    Aug 13 '20 at 14:52
  • $\begingroup$ I use Octave. IMHO it is very good software. It contains a symbolic math package, too. $\endgroup$ Aug 13 '20 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 '20 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 '20 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$ Aug 14 '20 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).$$

$\endgroup$
4
  • $\begingroup$ Could you expand this answer by giving the complete exact solution. $\endgroup$
    – cfdlab
    Aug 14 '20 at 13:15
  • $\begingroup$ Thanks, it is clear now. You can accept your own answer. $\endgroup$
    – cfdlab
    Aug 15 '20 at 5:03
  • $\begingroup$ What do you mean by accepting the answer? $\endgroup$ Aug 15 '20 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 '20 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.