# Problems in solving hydrogen atom energy levels

**URL:** <https://community.freefem.org/t/problems-in-solving-hydrogen-atom-energy-levels/1110>\
**Category:** General Discussion\
**Created:** [July 20, 2021, 3:59pm UTC](https://community.freefem.org/t/problems-in-solving-hydrogen-atom-energy-levels/1110 "2021-07-20T15:59:02Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ionizing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ionizing/32/648_2.png) [@Ionizing](https://community.freefem.org/u/Ionizing)\
**Post date:** [July 20, 2021, 3:59pm UTC](https://community.freefem.org/t/problems-in-solving-hydrogen-atom-energy-levels/1110/1 "2021-07-20T15:59:02Z")

</div>

Hi, it’s me again.

# Problem

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

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/68c45964b4e6588e5bec61dc92974094698e1bce.png)

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

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/3af4f103343af3dd6a3fa25adc1ac546966cb7ca.png)

Then defined the problem via

```freefem
varf a(u1, u2) = 
               // ka and kc are coefficients before nabla^2 and 1/r
    int3d(Th) ( -ka * Grad(u1)' * Grad(u2)
               // + 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](https://gist.github.com/Ionizing/4a58e5d94b4a57b01fca9e13479f3ae9), but it is too slow. This is why I turn to FEM.

# Code

Here is what I wrote:

```freefem
load "PETSc"
load "msh3"
load "medit"
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) = 
    int3d(Th) ( -ka * Grad(u1)' * Grad(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](https://community.freefem.org/uploads/short-url/biE95oQneQt7nSwxmSsM7tL2TYa.log) (253.3 KB)

---

<div class="post-metadata">

**Author:** ![Ionizing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ionizing/32/648_2.png) [@Ionizing](https://community.freefem.org/u/Ionizing)\
**Post date:** [July 22, 2021, 7:34am UTC](https://community.freefem.org/t/problems-in-solving-hydrogen-atom-energy-levels/1110/2 "2021-07-22T07:34:40Z")

</div>

OK, I’ve partly solved it by reading [User:Tclamb/FEniCS - Wikiversity](https://en.wikiversity.org/wiki/User:Tclamb/FEniCS#Quantum_Harmonic_Oscillator) .

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

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/d42e08085602a997d1836a3bb4e079c451d271bb.png)

That means there is no minus sign on the Laplace operator:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/b32d3be6c4991c0014ff4fa254c03ba7f57d67b9.png)  
` 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

```auto
 ---- 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 ==

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/49bcfe01015375f17b34e0658e7fb2246c0f393d.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/e0a404666a1025de073342578343c0a2eb5f27d9.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/2020a048f11ca9db32332df5cc6efe5e2353d0ef.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/9a8a41c94059c3c2047937fc42037c9dffab7327.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/b86cfe58b1c133efdc37f1b810f0b86fecc2a6c3.png)

# Code

```auto
load "PETSc"
load "msh3"
load "iovtk"

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) (
    Grad(u1)' * Grad(u2) * ka
    - 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]);
}

```

---

<div class="post-metadata">

**Author:** ![Ionizing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ionizing/32/648_2.png) [@Ionizing](https://community.freefem.org/u/Ionizing)\
**Post date:** [July 22, 2021, 7:41am UTC](https://community.freefem.org/t/problems-in-solving-hydrogen-atom-energy-levels/1110/3 "2021-07-22T07:41:17Z")

</div>

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

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/c0108725940961d5666b69867c2fb6deaa5af6cf.jpeg)

[Eigen\_1.zip](https://community.freefem.org/uploads/short-url/zs9Phqg8OocnhD3NCRHPtu0kCQi.zip) (2.1 MB)

---

<div class="post-metadata">

**Author:** ![Ionizing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ionizing/32/648_2.png) [@Ionizing](https://community.freefem.org/u/Ionizing)\
**Post date:** [July 22, 2021, 10:49am UTC](https://community.freefem.org/t/problems-in-solving-hydrogen-atom-energy-levels/1110/4 "2021-07-22T10:49:43Z")

</div>

Okay…

I’ve solved it by adding `order=[1]` in `savevtk` …