 # Problems in solving hydrogen atom energy levels

Hi, it’s me again.

# Problem

I want to solve the hydrogen atom using FF. The eigenvalue problem is like

I built a cylinder mesh (r = 15A, h = 30A) as the FE space:

Then defined the problem via

``````varf a(u1, u2) =
// ka and kc are coefficients before nabla^2 and 1/r
// + 0.0001 to avoid the singularity of 1/r at r=0
- (kc/sqrt(x^2 + y^2 + z^2 + 0.0001)) * u1 * u2);
varf b(u1, u2) =
int3d(Th) (u1 * u2);
``````

Since the ground state locates at -13.6eV, I set `sigma = -10.0`, but the code ended up with no solutions.

What should I do the make it work as expected?

P.S. When I ran it via `ff-mpirun -np 5`, 5 processes were running, but each process consumed only 50% CPU occupation.

In fact, I’ve implemented a finite difference version here, but it is too slow. This is why I turn to FEM.

# Code

Here is what I wrote:

``````load "PETSc"
include "macro_ddm.idp"

int nEV = getARGV("-nev", 10);
real sigma = -10;

real ka = 3.8099821161548593; // hbar^2 / 2m_e
real kc = 14.39964547842567; // e^2 / (4pi epsilon_0)

real r = 15.0;
border C(t=0, 2*pi) { x = r*cos(t); y = r*sin(t); label=1; }
mesh Baseh = buildmesh(C(60));

int[int] rup=[0, 1], rdown=[0, 2], rmid=[1,3];

func zmin=-15;
func zmax= 15;

int nlayer = 30;
mesh3 Th = buildlayers(Baseh, nlayer,
coef=1.0,
zbound = [zmin, zmax],
labelmid = rmid,
reffaceup = rup,
reffacelow = rdown
);

plot(Th, wait=1);

fespace Vh(Th, P1);
Vh u1, u2;

macro Grad(u) [dx(u), dy(u), dz(u)] // EOM

varf a(u1, u2) =
- (kc/sqrt(x^2 + y^2 + z^2 + 0.0001)) * u1 * u2);
varf b(u1, u2) =
int3d(Th) (u1 * u2);

matrix<real> A = a(Vh, Vh);
matrix<real> B = b(Vh, Vh);

Mat DistA(A);
Mat DistB(B);

real[int] EigenVAL(0);
Vh<real>[int] EigenVEC(1);

string ssparams =
" -eps_nev " + nEV +
" -eps_type krylovschur" +
" -eps_target " + sigma +
" -eps_view " +
" -eps_gen_hermitian "
;

int k = EPSSolve(
DistA,
DistB,
vectors = EigenVEC,
values  = EigenVAL,
sparams = ssparams
);

k = min(k, nEV);

if(mpirank) {
cout << "  Found " << k << " solutions" << endl;
}

Vh<real> Temp;

for (int i=0; i<k; ++i) {
if(!mpirank) {
cout << "  Eigenvalue \$" + i + " = " + EigenVAL[i] << endl;
}
Temp = EigenVEC[i];
if (!mpirank) {
plot(Temp, cmm = "Psi("+i+")  EV = " + EigenVAL[i],
wait=true, value=1, ps="eigen"+i+".eps");
}
}
``````

# Log

dump.log (253.3 KB)

OK, I’ve partly solved it by reading User:Tclamb/FEniCS - Wikiversity .

I need to transform the original PDE into a weak variational form.

That means there is no minus sign on the Laplace operator: ` int3d(Th) ( (hbar^2/2m) * Grad(u1)' * Grad(u2) )` is OK for it.

Then I ran it and got the results below. The eigenvalues seem to be consistent with the FDM I used. If I use a denser mesh on the center, the ground state energy would be much closer to the analytical result.

# Result

`````` ---- 0 -9.60472 ==
Plot::  Sorry no ps version for this type of plot 6
---- 1 -3.03749 ==
Plot::  Sorry no ps version for this type of plot 6
---- 2 -3.03749 ==
Plot::  Sorry no ps version for this type of plot 6
---- 3 -2.8008 ==
Plot::  Sorry no ps version for this type of plot 6
---- 4 -2.60847 ==
Plot::  Sorry no ps version for this type of plot 6
---- 5 -1.4506 ==
Plot::  Sorry no ps version for this type of plot 6
---- 6 -1.4506 ==
Plot::  Sorry no ps version for this type of plot 6
---- 7 -1.40484 ==
Plot::  Sorry no ps version for this type of plot 6
---- 8 -1.40484 ==
Plot::  Sorry no ps version for this type of plot 6
---- 9 -1.38226 ==
---- 10 -1.36496 ==
---- 11 -1.36496 ==
---- 12 -1.26952 ==
---- 13 -1.24226 ==
---- 14 -0.869513 ==
---- 15 -0.859284 ==
---- 16 -0.857945 ==
---- 17 -0.855488 ==
---- 18 -0.848229 ==
---- 19 -0.84818 ==
``````

# Code

``````load "PETSc"

int nn = 30;

real ka = 3.8099821161548593; // hbar^2 / 2m_e
real kc = 14.39964547842567; // e^2 / (4pi epsilon_0)

mesh3 Th = cube(nn, nn, nn, [30*x-15, 30*y-15, 30*z-15]);
plot(Th, wait=1);

fespace Vh(Th, P1);

cout << "Th :  nv = " << Th.nv << " nt = " << Th.nt << endl;

real sigma = -13;

macro Grad(u) [dx(u), dy(u), dz(u)] // EOM
varf a(u1, u2) = int3d(Th) (
- 1.0 / sqrt(x^2 + y^2 + z^2)  * u1 * u2 * kc
- sigma * u1 * u2
)
//+ on(1, 2, 3, 4, 5, 6, u1=0.0)
;

varf b([u1], [u2]) = int3d(Th) ( u1 * u2 ) ;

matrix A = a(Vh, Vh);
matrix B = b(Vh, Vh);

int nev = 20;
real[int] ev(nev);
Vh[int] eV(nev);

int k = EigenValue(A, B, sym=true, sigma=sigma, value=ev, vector=eV, tol=1e-10);

//string ssparams = "-eps_nev " + nev
//+ " -eps_type krylovschur"
//+ " -eps_target " + sigma
//+ " -eps_view "
//+ " -eps_gen_hermitian";

//Mat DistA(A);
//Mat DistB(B);

//int k = EPSSolve(
//DistA, DistB,
//vectors = eV,
//values = ev,
//sparams = ssparams
//);

k = min(k, nev);

for (int i=0; i<k; ++i) {
cout << " ---- " << i << " " << ev[i] << " == " << endl;
plot(eV[i], cmm="#" + i + " EigenValue=" + ev[i], wait=true);
savevtk("Eigen_" + i + ".vtk", Th, eV[i], dataname="EigenValue=" + ev[i]);
}
``````

Well, are there any tutorials about how to view the 3D eigenvectors using ParaView? I used `savevtk` but could view nothing inside it …

Eigen_1.zip (2.1 MB)

Okay…

I’ve solved it by adding `order=` in `savevtk`