The density employed rho (x,y,z) ,is that of the hydrogen electron in the 1s state.The corresponding potential is known phi 1s(x,y,z) and compared with.
c POISSON IN THREE DIM (X,Y,Z) ,calculation of potential for one
c electron density distribution rho=psi1s(x,y,z)**2 , analytical
c solution is phi1s(r)-see below
c phi1s(r)=(4.*a**3/r)*(.25/a**3-exp(-2.*a*r)*(.25*r/a**2+.25/a**3))
c nucleus is located at the origin
equivalence (limx1,x1),(x1,y1,z1),(dx,dy,dz),(limy1,limx1,limz1),
$ (nx,ny,nz)
real limx1,limx2,limy1,limy2,limz1,limz2,lx
dimension phi(0:40,0:40,0:40),phi1(0:40,0:40,0:40)
r(x,y,z)=sqrt(x**2+y**2+z**2)+1.e-4
psi1s(x,y,z)=sqrt(alfa**3/pi)*exp(-alfa*r(x,y,z))
rho(x,y,z)=psi1s(x,y,z)**2
phi1s(x,y,z)=(4.*alfa**3/r(x,y,z))*(.25/alfa**3-
$ exp(-2.*alfa*r(x,y,z))*(.25*r(x,y,z)/alfa**2+.25/alfa**3))
pi=2.*asin(1.)
alfa=1.
c limits of integration x.le.+Lx x.ge.-Lx , same for y and z
lx=6.
nx=20
dx=2.*lx/float(nx)
print*,'nx=',nx
iter=20
x1=-lx
limx1=-lx
c 3 dim. integrations
c Phi at the boundaries = 1./r
c boundary at plane x=-lx ,x=+lx
do 205 j=0,ny
do 205 k=0,nz
x=limx1
y=y1+dy*float(j)
z=z1+dz*float(k)
r1=sqrt(limx1**2+y**2+z**2)
phi(0,j,k)=1./r1
phi(nx,j,k)= 1./r1
phi1(0,j,k)=1./r1
phi1(nx,j,k)= 1./r1
205 continue
c boundary at plane y=limy1
do 206 i=0,nx
do 206 k=0,nz
y=limy1
x=x1+dx*float(i)
z=z1+dz*float(k)
r1=sqrt(x**2+limy1**2+z**2)
phi(i,0,k)=1./r1
phi(i,ny,k)=1./r1
phi1(i,0,k)=1./r1
phi1(i,ny,k)=1./r1
206 continue
c boundary at plane z=limz1
do 207 i=0,nx
do 207 j=0,ny
z=limz1
x=x1+dx*float(i)
y=y1+dy*float(j)
r1=sqrt(x**2+y**2+limz1**2)
phi(i,j,0)=1./r1
phi(i,j,nz)=1./r1
phi1(i,j,0)=1./r1
phi1(i,j,nz)=1./r1
207 continue
c initial value of the potential phi(i,j,k) (at interior points)
do 30 iz=1,nz-1
do 30 iy=1,ny-1
do 30 ix=1,nx-1
x=x1+dx*float(ix)
y=y1+dy*float(iy)
z=z1+dz*float(iz)
c if(r(x,y,z).gt.0.0.and.r(x,y,z).le.1.)then
if(r(x,y,z).eq.0.0)then
phi(ix,iy,iz)=1./sqrt(3.*dx**2)
c phi(ix,iy,iz)= phi1s(x,y,z)
endif
c if(r(x,y,z).gt.1.)then
c phi(ix,iy,iz)=1.
phi(ix,iy,iz)=1./r(x,y,z)
c phi(ix,iy,iz)=phi1s(x,y,z)
c endif
c print*,'x,y,z,phi=',x,y,z,phi(ix,iy,iz)
30 continue
c**************
c integration using finite difference
do 55 it=1,iter
z=z1
do 35 iz=1,nz-1
y=y1
do 25 iy= 1,ny-1
x=x1
do 15 ix=1,nx-1
trho=(1./6.)*4.*pi*dx**2*rho(x,y,z)
sumphi=(1./6.)*( phi(ix+1,iy,iz)+phi(ix-1,iy,iz)+
$ phi(ix,iy+1,iz)+phi(ix,iy-1,iz)+phi(ix,iy,iz+1)+phi(ix,iy,iz-1))
phi1(ix,iy,iz)= trho+sumphi
x=x1+dx*float(ix)
15 continue
y=y1+dy*float(iy)
25 continue
z=z1+dz*float(iz)
35 continue
c relabel phi and phi1
do 45 i=0,nx
do 45 j=0,ny
do 45 k=0,nz
phi(i,j,k)=phi1(i,j,k)
45 continue
55 continue
c check phi vs phi1s
do 65 i=1,nx-1,4
do 65 j=1,ny-1,4
do 65 k=1,nz-1,4
x=x1+dx*float(i)
y=y1+dy*float(j)
z=z1+dz*float(k)
print 130,x,y,z,r(x,y,z)
print 100,phi(i,j,k),phi1s(x,y,z)
print*,' '
65 continue
100 format(2x,'phinum, phi1s=',3x,e10.3,3x,e10.3)
print*,' '
130 format(2x,'x,y,z,r=',3(1x,e11.4),1x,e10.3)
stop
end
RUN
x,y,z,r= -0.3000E+01 -0.5400E+01 -0.3000E+01 0.687E+01
phinum, phi1s= 0.146E+00 0.146E+00
x,y,z,r= -0.3000E+01 -0.5400E+01 -0.6000E+00 0.621E+01
phinum, phi1s= 0.161E+00 0.161E+00
x,y,z,r= -0.3000E+01 -0.5400E+01 0.1800E+01 0.643E+01
phinum, phi1s= 0.155E+00 0.155E+00
x,y,z,r= -0.3000E+01 -0.5400E+01 0.4200E+01 0.747E+01
phinum, phi1s= 0.134E+00 0.134E+00
x,y,z,r= -0.3000E+01 -0.3000E+01 -0.5400E+01 0.687E+01
phinum, phi1s= 0.146E+00 0.146E+00
x,y,z,r= -0.3000E+01 -0.3000E+01 -0.3000E+01 0.520E+01
phinum, phi1s= 0.192E+00 0.192E+00
x,y,z,r= -0.3000E+01 -0.3000E+01 -0.6000E+00 0.428E+01
phinum, phi1s= 0.233E+00 0.233E+00
x,y,z,r= -0.3000E+01 -0.3000E+01 0.1800E+01 0.461E+01
phinum, phi1s= 0.217E+00 0.217E+00
x,y,z,r= -0.3000E+01 -0.3000E+01 0.4200E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= -0.3000E+01 -0.6000E+00 -0.5400E+01 0.621E+01
phinum, phi1s= 0.161E+00 0.161E+00
x,y,z,r= -0.3000E+01 -0.6000E+00 -0.3000E+01 0.428E+01
phinum, phi1s= 0.233E+00 0.233E+00
x,y,z,r= -0.3000E+01 -0.6000E+00 -0.6000E+00 0.312E+01
phinum, phi1s= 0.314E+00 0.318E+00
x,y,z,r= -0.3000E+01 -0.6000E+00 0.1800E+01 0.355E+01
phinum, phi1s= 0.283E+00 0.281E+00
x,y,z,r= -0.3000E+01 -0.6000E+00 0.4200E+01 0.520E+01
phinum, phi1s= 0.193E+00 0.192E+00
x,y,z,r= -0.3000E+01 0.1800E+01 -0.5400E+01 0.643E+01
phinum, phi1s= 0.155E+00 0.155E+00
x,y,z,r= -0.3000E+01 0.1800E+01 -0.3000E+01 0.461E+01
phinum, phi1s= 0.217E+00 0.217E+00
x,y,z,r= -0.3000E+01 0.1800E+01 -0.6000E+00 0.355E+01
phinum, phi1s= 0.283E+00 0.281E+00
x,y,z,r= -0.3000E+01 0.1800E+01 0.1800E+01 0.393E+01
phinum, phi1s= 0.259E+00 0.254E+00
x,y,z,r= -0.3000E+01 0.1800E+01 0.4200E+01 0.547E+01
phinum, phi1s= 0.184E+00 0.183E+00
x,y,z,r= -0.3000E+01 0.4200E+01 -0.5400E+01 0.747E+01
phinum, phi1s= 0.134E+00 0.134E+00
x,y,z,r= -0.3000E+01 0.4200E+01 -0.3000E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= -0.3000E+01 0.4200E+01 -0.6000E+00 0.520E+01
phinum, phi1s= 0.193E+00 0.192E+00
x,y,z,r= -0.3000E+01 0.4200E+01 0.1800E+01 0.547E+01
phinum, phi1s= 0.184E+00 0.183E+00
x,y,z,r= -0.3000E+01 0.4200E+01 0.4200E+01 0.665E+01
phinum, phi1s= 0.150E+00 0.150E+00
x,y,z,r= -0.6000E+00 -0.5400E+01 -0.5400E+01 0.766E+01
phinum, phi1s= 0.131E+00 0.131E+00
x,y,z,r= -0.6000E+00 -0.5400E+01 -0.3000E+01 0.621E+01
phinum, phi1s= 0.161E+00 0.161E+00
x,y,z,r= -0.6000E+00 -0.5400E+01 -0.6000E+00 0.547E+01
phinum, phi1s= 0.183E+00 0.183E+00
x,y,z,r= -0.6000E+00 -0.5400E+01 0.1800E+01 0.572E+01
phinum, phi1s= 0.175E+00 0.175E+00
x,y,z,r= -0.6000E+00 -0.5400E+01 0.4200E+01 0.687E+01
phinum, phi1s= 0.146E+00 0.146E+00
x,y,z,r= -0.6000E+00 -0.3000E+01 -0.5400E+01 0.621E+01
phinum, phi1s= 0.161E+00 0.161E+00
x,y,z,r= -0.6000E+00 -0.3000E+01 -0.3000E+01 0.428E+01
phinum, phi1s= 0.233E+00 0.233E+00
x,y,z,r= -0.6000E+00 -0.3000E+01 -0.6000E+00 0.312E+01
phinum, phi1s= 0.314E+00 0.318E+00
x,y,z,r= -0.6000E+00 -0.3000E+01 0.1800E+01 0.355E+01
phinum, phi1s= 0.283E+00 0.281E+00
x,y,z,r= -0.6000E+00 -0.3000E+01 0.4200E+01 0.520E+01
phinum, phi1s= 0.193E+00 0.192E+00
x,y,z,r= -0.6000E+00 -0.6000E+00 -0.5400E+01 0.547E+01
phinum, phi1s= 0.183E+00 0.183E+00
x,y,z,r= -0.6000E+00 -0.6000E+00 -0.3000E+01 0.312E+01
phinum, phi1s= 0.314E+00 0.318E+00
x,y,z,r= -0.6000E+00 -0.6000E+00 -0.6000E+00 0.104E+01
phinum, phi1s= 0.571E+00 0.717E+00
x,y,z,r= -0.6000E+00 -0.6000E+00 0.1800E+01 0.199E+01
phinum, phi1s= 0.499E+00 0.474E+00
x,y,z,r= -0.6000E+00 -0.6000E+00 0.4200E+01 0.428E+01
phinum, phi1s= 0.240E+00 0.233E+00
x,y,z,r= -0.6000E+00 0.1800E+01 -0.5400E+01 0.572E+01
phinum, phi1s= 0.175E+00 0.175E+00
x,y,z,r= -0.6000E+00 0.1800E+01 -0.3000E+01 0.355E+01
phinum, phi1s= 0.283E+00 0.281E+00
x,y,z,r= -0.6000E+00 0.1800E+01 -0.6000E+00 0.199E+01
phinum, phi1s= 0.499E+00 0.474E+00
x,y,z,r= -0.6000E+00 0.1800E+01 0.1800E+01 0.262E+01
phinum, phi1s= 0.446E+00 0.375E+00
x,y,z,r= -0.6000E+00 0.1800E+01 0.4200E+01 0.461E+01
phinum, phi1s= 0.224E+00 0.217E+00
x,y,z,r= -0.6000E+00 0.4200E+01 -0.5400E+01 0.687E+01
phinum, phi1s= 0.146E+00 0.146E+00
x,y,z,r= -0.6000E+00 0.4200E+01 -0.3000E+01 0.520E+01
phinum, phi1s= 0.193E+00 0.192E+00
x,y,z,r= -0.6000E+00 0.4200E+01 -0.6000E+00 0.428E+01
phinum, phi1s= 0.240E+00 0.233E+00
x,y,z,r= -0.6000E+00 0.4200E+01 0.1800E+01 0.461E+01
phinum, phi1s= 0.224E+00 0.217E+00
x,y,z,r= -0.6000E+00 0.4200E+01 0.4200E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= 0.1800E+01 -0.5400E+01 -0.5400E+01 0.785E+01
phinum, phi1s= 0.127E+00 0.127E+00
x,y,z,r= 0.1800E+01 -0.5400E+01 -0.3000E+01 0.643E+01
phinum, phi1s= 0.155E+00 0.155E+00
x,y,z,r= 0.1800E+01 -0.5400E+01 -0.6000E+00 0.572E+01
phinum, phi1s= 0.175E+00 0.175E+00
x,y,z,r= 0.1800E+01 -0.5400E+01 0.1800E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= 0.1800E+01 -0.5400E+01 0.4200E+01 0.707E+01
phinum, phi1s= 0.141E+00 0.141E+00
x,y,z,r= 0.1800E+01 -0.3000E+01 -0.5400E+01 0.643E+01
phinum, phi1s= 0.155E+00 0.155E+00
x,y,z,r= 0.1800E+01 -0.3000E+01 -0.3000E+01 0.461E+01
phinum, phi1s= 0.217E+00 0.217E+00
x,y,z,r= 0.1800E+01 -0.3000E+01 -0.6000E+00 0.355E+01
phinum, phi1s= 0.283E+00 0.281E+00
x,y,z,r= 0.1800E+01 -0.3000E+01 0.1800E+01 0.393E+01
phinum, phi1s= 0.259E+00 0.254E+00
x,y,z,r= 0.1800E+01 -0.3000E+01 0.4200E+01 0.547E+01
phinum, phi1s= 0.184E+00 0.183E+00
x,y,z,r= 0.1800E+01 -0.6000E+00 -0.5400E+01 0.572E+01
phinum, phi1s= 0.175E+00 0.175E+00
x,y,z,r= 0.1800E+01 -0.6000E+00 -0.3000E+01 0.355E+01
phinum, phi1s= 0.283E+00 0.281E+00
x,y,z,r= 0.1800E+01 -0.6000E+00 -0.6000E+00 0.199E+01
phinum, phi1s= 0.499E+00 0.474E+00
x,y,z,r= 0.1800E+01 -0.6000E+00 0.1800E+01 0.262E+01
phinum, phi1s= 0.446E+00 0.375E+00
x,y,z,r= 0.1800E+01 -0.6000E+00 0.4200E+01 0.461E+01
phinum, phi1s= 0.224E+00 0.217E+00
x,y,z,r= 0.1800E+01 0.1800E+01 -0.5400E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= 0.1800E+01 0.1800E+01 -0.3000E+01 0.393E+01
phinum, phi1s= 0.259E+00 0.254E+00
x,y,z,r= 0.1800E+01 0.1800E+01 -0.6000E+00 0.262E+01
phinum, phi1s= 0.446E+00 0.375E+00
x,y,z,r= 0.1800E+01 0.1800E+01 0.1800E+01 0.312E+01
phinum, phi1s= 0.406E+00 0.318E+00
x,y,z,r= 0.1800E+01 0.1800E+01 0.4200E+01 0.491E+01
phinum, phi1s= 0.211E+00 0.204E+00
x,y,z,r= 0.1800E+01 0.4200E+01 -0.5400E+01 0.707E+01
phinum, phi1s= 0.141E+00 0.141E+00
x,y,z,r= 0.1800E+01 0.4200E+01 -0.3000E+01 0.547E+01
phinum, phi1s= 0.184E+00 0.183E+00
x,y,z,r= 0.1800E+01 0.4200E+01 -0.6000E+00 0.461E+01
phinum, phi1s= 0.224E+00 0.217E+00
x,y,z,r= 0.1800E+01 0.4200E+01 0.1800E+01 0.491E+01
phinum, phi1s= 0.211E+00 0.204E+00
x,y,z,r= 0.1800E+01 0.4200E+01 0.4200E+01 0.621E+01
phinum, phi1s= 0.162E+00 0.161E+00
x,y,z,r= 0.4200E+01 -0.5400E+01 -0.5400E+01 0.872E+01
phinum, phi1s= 0.115E+00 0.115E+00
x,y,z,r= 0.4200E+01 -0.5400E+01 -0.3000E+01 0.747E+01
phinum, phi1s= 0.134E+00 0.134E+00
x,y,z,r= 0.4200E+01 -0.5400E+01 -0.6000E+00 0.687E+01
phinum, phi1s= 0.146E+00 0.146E+00
x,y,z,r= 0.4200E+01 -0.5400E+01 0.1800E+01 0.707E+01
phinum, phi1s= 0.141E+00 0.141E+00
x,y,z,r= 0.4200E+01 -0.5400E+01 0.4200E+01 0.803E+01
phinum, phi1s= 0.125E+00 0.125E+00
x,y,z,r= 0.4200E+01 -0.3000E+01 -0.5400E+01 0.747E+01
phinum, phi1s= 0.134E+00 0.134E+00
x,y,z,r= 0.4200E+01 -0.3000E+01 -0.3000E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= 0.4200E+01 -0.3000E+01 -0.6000E+00 0.520E+01
phinum, phi1s= 0.193E+00 0.192E+00
x,y,z,r= 0.4200E+01 -0.3000E+01 0.1800E+01 0.547E+01
phinum, phi1s= 0.184E+00 0.183E+00
x,y,z,r= 0.4200E+01 -0.3000E+01 0.4200E+01 0.665E+01
phinum, phi1s= 0.150E+00 0.150E+00
x,y,z,r= 0.4200E+01 -0.6000E+00 -0.5400E+01 0.687E+01
phinum, phi1s= 0.146E+00 0.146E+00
x,y,z,r= 0.4200E+01 -0.6000E+00 -0.3000E+01 0.520E+01
phinum, phi1s= 0.193E+00 0.192E+00
x,y,z,r= 0.4200E+01 -0.6000E+00 -0.6000E+00 0.428E+01
phinum, phi1s= 0.240E+00 0.233E+00
x,y,z,r= 0.4200E+01 -0.6000E+00 0.1800E+01 0.461E+01
phinum, phi1s= 0.224E+00 0.217E+00
x,y,z,r= 0.4200E+01 -0.6000E+00 0.4200E+01 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= 0.4200E+01 0.1800E+01 -0.5400E+01 0.707E+01
phinum, phi1s= 0.141E+00 0.141E+00
x,y,z,r= 0.4200E+01 0.1800E+01 -0.3000E+01 0.547E+01
phinum, phi1s= 0.184E+00 0.183E+00
x,y,z,r= 0.4200E+01 0.1800E+01 -0.6000E+00 0.461E+01
phinum, phi1s= 0.224E+00 0.217E+00
x,y,z,r= 0.4200E+01 0.1800E+01 0.1800E+01 0.491E+01
phinum, phi1s= 0.211E+00 0.204E+00
x,y,z,r= 0.4200E+01 0.1800E+01 0.4200E+01 0.621E+01
phinum, phi1s= 0.162E+00 0.161E+00
x,y,z,r= 0.4200E+01 0.4200E+01 -0.5400E+01 0.803E+01
phinum, phi1s= 0.125E+00 0.125E+00
x,y,z,r= 0.4200E+01 0.4200E+01 -0.3000E+01 0.665E+01
phinum, phi1s= 0.150E+00 0.150E+00
x,y,z,r= 0.4200E+01 0.4200E+01 -0.6000E+00 0.597E+01
phinum, phi1s= 0.168E+00 0.167E+00
x,y,z,r= 0.4200E+01 0.4200E+01 0.1800E+01 0.621E+01
phinum, phi1s= 0.162E+00 0.161E+00
x,y,z,r= 0.4200E+01 0.4200E+01 0.4200E+01 0.727E+01
phinum, phi1s= 0.138E+00 0.137E+00