Problems in solving hydrogen atom energy levels

General Discussion
1

Hi, it’s me again.

Problem

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

image
image1448×232 15.9 KB

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

image
image1226×1510 57.2 KB

Then defined the problem via

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, but it is too slow. This is why I turn to FEM.

Code

Here is what I wrote:

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 (253.3 KB)

2

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

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

image
image1574×734 85.4 KB

That means there is no minus sign on the Laplace operator:
image
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 ==

image
image1144×1254 92.5 KB

image
image1144×1254 95.1 KB

image
image1144×1254 111 KB

image
image1144×1254 107 KB

image
image1144×1254 96.2 KB

Code

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]);
}
3

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

image
image2164×2066 399 KB

Eigen_1.zip (2.1 MB)

4

Okay…

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