How to get more exact solution of 3D Schrodinger equation

I use following code to solve 3D Schrodinger equation of Hydrogen atom.

real me=9.10938356e-31;
real planckh=6.62607004e-34;
real hbar=1.0545718e-34;
real kgm2Ds2ToeV=6.24150636309e18;
real a0=5.2917721039e-11;//bohr radius
real qe=1.6021766208e-19;//electron mass
real varepsilon0=8.854187817e-12;

load "msh3"

//网格
int k=5;
int n=40;
//[-k*a0, k*a0]^3
mesh3 Th = cube(n, n, n,[-k*a0+2*k*a0*x, -k*a0+2*k*a0*y, -k*a0+2*k*a0*z]);
fespace Vh(Th,P1);
Vh Psi, v;

int nev=3;
real tilE=-2e-18;
Vh[int] EigenPsi(nev);
real[int] EigenVal(nev);

varf Schrodinger(Psi,v) =
	int3d(Th)( hbar^2/(2*me)*(dx(Psi)*dx(v) + dy(Psi)*dy(v) + dz(Psi)*dz(v)))
	- int3d(Th)( qe^2/(4*pi*varepsilon0*sqrt(x^2+y^2+z^2))*Psi*v)
	- int3d(Th)(tilE*Psi*v) + on(1,2,3,4,5,6,Psi=0);
varf RHS(Psi,v) = int3d(Th)(Psi*v);
matrix A = Schrodinger(Vh,Vh,solver="SPARSESOLVER");
matrix B = RHS(Vh,Vh,solver=CG,eps=1e-20);

// 求解接近于sigma的特征值Ax=lambda*B*x,左边为A-sigma*B,右边为B。
int num = EigenValue(A,B,sym=true,sigma=tilE,value=EigenVal,
vector=EigenPsi,tol=1e-20,maxit=2000,ncv=0);

for(int i = 0; i < nev; i++){
	cout << "Eigenvalue["+i+"] = "<<EigenVal[i]<<" , Energy level: " << EigenVal[i]*kgm2Ds2ToeV << " eV" << endl;
}

The result is -13.20eV, -1.1697eV, -1.1697eV which are different from exact solution -13.6eV, -3.4eV, -1.51eV.

I thought the problem may due to the magnititude of variable. But I don’t know how to correctly scale the equation.

My code is based on QuantumHarmonicOscillator

Use SLEPc, much easier to finely tune than ARPACK. Additionally, the shift (tilE) should not be in the Schrodinger variational formulation.

What do you mean by “should not be in”? May I just remove that int3d term?

Yes, at least that is what you should do if using SLEPc.

I am not sure about if your question is about the quantitative values, or is it qualitatively…

The fact that you have the same value twice ( -1.1697eV, -1.1697eV ) makes sense. Actually you should find four of them with approximately the same energy level.

Remember that the energy levels are degenerate. so you should find:

  • 1 state with lowest eigenvalue (the ground state 1s)
  • 4 states with a higher EV [ 1 for the 2s singlet + 3 for the 2p triplet]
    -…
    Of course, the numerics will slightly lift the degeneracy…

I hope this help

I try to set nev=10. Now there is one -12 eV and four approximate -3.2 eV and five -1.4 eV which coincides with what you said. Thank you.

I want to get result that matching analytical solution from FreeFem++. But increasing mesh range(k) and use more segments(n) doesn’t help much. There are always some gaps.

I have check examples/hpddm which contains some SLEPc examples. But I find it’s too complex and I can’t understand yet. I would like to try it after I’m more familar with FreeFEM++.

You can look at Pierre Jolivet introduction to parallel FreeFEM part 1 - YouTube to figure out what is going on.

Thank you, I will learn that.

I managed to write a SLEPc version that using EPSSolve. But the scale of spatial region is very tiny (about 1e-10), if I directly use EPSSolve, there will be error zero pivot. It looks like that SLEPc just ignore the problem.I have to rewrite th equation to scale.

How ever the original problem can be solved by ARPACK directly. Is there any way to scale the matrix in SLEPc automatically without rewrite the equation? I tried to search in SLEPc documentation but didn’t find any clue.

The code to solve original problem using SLEPc is( there will be error):

//Run: freefem++-mpi Schdroger3D.edp 

real me=9.10938356e-31;  //electron mass
real planckh=6.62607004e-34;
real hbar=1.0545718e-34;
real kgm2Ds2ToeV=6.24150636309e18;
real a0=5.2917721039e-11;//bohr radius
real qe=1.6021766208e-19;//electron charge
real varepsilon0=8.854187817e-12;

load "UMFPACK64"
load "msh3"
load "TetGen"
load "medit"
include "MeshSurface.idp"
load "PETSc"
macro dimension()3 // EOM
include "macro_ddm.idp"
//load "MUMPS_mpi"

int k=17;
int n=20;
real hs1 = a0; //mesh size on sphere
real hs2 = a0; //mesh size on sphere
func Pk=P1;

//正方形

int[int] nn = [n, n, n];
//Box: [-k*a0, k*a0]^3
real [int,int] Box = [[-k*a0, k*a0], [-k*a0, k*a0], [-k*a0, k*a0]]; //Bounding box: [xmin, xmax]
int [int,int] L = [[1, 2], [3, 4], [5, 6]]; //label: left,right, front, back, down, right
meshS ThH = SurfaceHex(nn, Box, L, 1);

meshS ThS1 = Sphere(1.5*a0, hs1, 7, 1);
meshS ThS2 = Sphere(5*a0, hs2, 8, 1);

meshS ThHS = ThH + ThS1 +ThS2;
//medit("Hex-Sphere", ThHS);
 
//Mesh
real[int] domain = [0, 0, 0,     1, hs1^3/6, 
					2.5*a0, 0,0, 2, hs2^3/6, 
					10*a0,0,0,   3, (2*k*a0/nn[1])^3/6];
mesh3 Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=3, regionlist=domain);
//Th=change(Th,fregion=nuTriangle%mpisize); 
//medit("Cube with ball", Th);

fespace Vh(Th,Pk);
Vh Psi, v;
//real tilE=-2e-18;

varf Schrodinger(Psi,v) =
	int3d(Th)( hbar^2/(2*me)*(dx(Psi)*dx(v) + dy(Psi)*dy(v) + dz(Psi)*dz(v)))
	- int3d(Th)( qe^2/(4*pi*varepsilon0*sqrt(x^2+y^2+z^2))*Psi*v);
	//+ on(1,2,3,4,5,6,Psi=0);
	//- int3d(Th, mpirank)(tilE*Psi*v) 
	
	//+ on(1,2,3,4,5,6,Psi=0);
varf RHS(Psi,v) = int3d(Th)(Psi*v);

matrix locA = Schrodinger(Vh,Vh,solver="SPARSESOLVER");
matrix locB = RHS(Vh,Vh,solver="SPARSESOLVER",eps=1e-20);

Mat A;
createMat(Th, A, Pk)
A=locA;
Mat B(A, locB);

// 求解接近于sigma的特征值Ax=lambda*B*x,左边为A-sigma*B,右边为B。
int nev=10;
Vh[int] EigenPsi(nev);
real[int] EigenVal(nev);
int num = EPSSolve(A,B, sparams = " -st_type sinvert -eps_nev "+nev
			+ "  -eps_gen_hermitian"
			+ " -eps_target -16.0",
		 values = EigenVal, vectors = EigenPsi);

//打印能级
for(int i = 0; i < nev; i++){
	cout << "Eigenvalue["+i+"] = "<<EigenVal[i]
	<<" , Energy level: " << EigenVal[i]*kgm2Ds2ToeV << " eV" << endl;
}

Scaled version is:

//verbosity=0;

include "constants.edp";

load "UMFPACK64"
load "msh3"
load "TetGen"
load "medit"
include "MeshSurface.idp";
load "PETSc"
macro dimension()3 // EOM
include "macro_ddm.idp";

//网格
real R=20;   //Original range: [-20a0, 20a0]^3
real a=8.0; //2/(a0*k)=a
	real k=2.0/(a0*a);  //Scaled range : [-k*20*a0, k*20*a0]^3
	real l=k*R*a0;
int n=25;
real R1=0.1*l;
real R2=0.4*l;
real hs1= 0.04*R1;
real hs2=0.1*R2;
func Pk=P1;


int[int] nn = [n, n, n];
//Box: [-k*a0, k*a0]^3
real [int,int] Box = [[-l, l], [-l, l], [-l, l]]; //Bounding box: [xmin, xmax]
int [int,int] L = [[1, 2], [3, 4], [5, 6]]; //label: left,right, front, back, down, right
meshS ThH = SurfaceHex(nn, Box, L, 1);

meshS ThS1 = Sphere(R1, hs1, 7, 1);
meshS ThS2 = Sphere(R2, hs2, 8, 1);

meshS ThHS = ThH + ThS1 +ThS2;
//medit("Hex-Sphere", ThHS);
 
//Mesh
real[int] domain = [0, 0, 0,     1, hs1^3/6, 
					(R1+R2)/2, 0,0, 2, hs2^3/6, 
					(R2+l)/2,0,0,   3, (2*l/nn[1])^3/6];
mesh3 Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=3, regionlist=domain);
medit("All", Th);

fespace Vh(Th,Pk);
Vh Psi, v;

varf Schrodinger(Psi,v) =
	int3d(Th)(dx(Psi)*dx(v) + dy(Psi)*dy(v) + dz(Psi)*dz(v) - a/sqrt(x^2+y^2+z^2)*Psi*v  ) 
	 //- int3d(Th)(tilE*Psi*v) 
	 //+ on(1,2,3,4,5,6,Psi=0)
	 ;
varf RHS(Psi,v) = int3d(Th)(Psi*v);

matrix locA = Schrodinger(Vh,Vh,solver="SPARSESOLVER");
matrix locB = RHS(Vh,Vh,solver="SPARSESOLVER",eps=1e-20);

Mat A;
createMat(Th, A, Pk)
A=locA;
Mat B(A, locB);

// 求解接近于sigma的特征值Ax=lambda*B*x,左边为A-sigma*B,右边为B。
int nev=10;
Vh[int] EigenPsi(nev);
real[int] EigenVal(nev);
int num = EPSSolve(A,B, sparams = " -st_type sinvert -eps_nev "+nev
			+ "  -eps_gen_hermitian"
			+ " -eps_target -16.0",
		 values = EigenVal, vectors = EigenPsi);

for(int i = 0; i < nev; i++){
	cout << "Eigenvalue["+i+"] = "<<EigenVal[i] 
	<<" , Energy level: " << EigenVal[i]*hbar^2*k^2/(2*me)*kgm2Ds2ToeV << " eV" << endl;
}

Finally I come up of a parallel SLEPc version. My idea it to create a very fine mesh in center part and rescale the original Schrodinger equation.
constants.edp:

real me=9.10938356e-31;
real pluckh=6.62607004e-34;
real hbar=1.0545718e-34;
real kgm2Ds2ToeV=6.24150636309e18;
real a0=5.2917721039e-11;//bohr radius
real qe=1.6021766208e-19;//electron mass
real varepsilon0=8.854187817e-12;

Schrodinger3Ds.edp:

// mpiexec -np 4 freefem++-mpi Schrodinger3Ds.edp
// 需要安装MSMPI和Intel MPI

include "constants.edp";

load "UMFPACK64"
load "msh3"
load "TetGen"
load "medit"
include "MeshSurface.idp";
load "PETSc"
macro dimension()3 // EOM
include "macro_ddm.idp";

//设置
real R=17;    //Original range: [-20a0, 20a0]^3
real a=8.0;   //2/(a0*k)=a
	real k=2.0/(a0*a);  //Scaled range : [-k*R*a0, k*R*a0]^3
	real l=k*R*a0;
int n=25;
real R1=0.15*l;
real R2=0.4*l;
real hs1= 0.03*R1;
real hs2=0.07*R2;
func Pk=P1;

real sigma=-16.0;


int[int] nn = [n, n, n];
//Box: [-k*a0, k*a0]^3
real [int,int] Box = [[-l, l], [-l, l], [-l, l]]; //Bounding box: [xmin, xmax]
int [int,int] L = [[1, 2], [3, 4], [5, 6]]; //label: left,right, front, back, down, right
meshS ThH = SurfaceHex(nn, Box, L, 1);

meshS ThS1 = Sphere(R1, hs1, 7, 1);
meshS ThS2 = Sphere(R2, hs2, 8, 1);

meshS ThHS = ThH + ThS1 +ThS2;
 
//Mesh
real[int] domain = [0, 0, 0,     1, hs1^3/6, 
					(R1+R2)/2, 0,0, 2, hs2^3/6, 
					(R2+l)/2,0,0,   3, (2*l/nn[1])^3/6];
mesh3 Th=cube(3, 4, 5);

int[int][int] intersection;   // Local-to-neighbors renumbering
real[int] D;                  // Partition of unity
{
    Th = tetg(ThHS, switch="pqaAAYYQ", nbofregions=3, regionlist=domain);
    build(Th,       // The local mesh
            1,        // Refinement factor
            intersection, // local-to-neighbors renumbering
            D,        // partition of unity
            Pk,           // FE-space
            mpiCommWorld // Communicator
         );
}
 
//medit("All", Th);

fespace Vh(Th,Pk);
Vh Psi, v;

varf Schrodinger(Psi,v) =
	int3d(Th)(dx(Psi)*dx(v) + dy(Psi)*dy(v) + dz(Psi)*dz(v) 
	- (a/sqrt(x^2+y^2+z^2) )*Psi*v  ) 
	;
varf RHS(Psi,v) = int3d(Th)(Psi*v);

matrix locA = Schrodinger(Vh,Vh,solver="SPARSESOLVER");
matrix locB = RHS(Vh,Vh,solver="SPARSESOLVER",eps=1e-20);

Mat A(locA, intersection, D);
Mat B(A, locB);

// 求解接近-16的特征值,因为是变换后的,-16的特征值对应-13.6eV的能级。
int nev=10;
Vh[int] EigenPsi(nev);
real[int] EigenVal(nev);
int num = EPSSolve(A,B, sparams = " -eps_type krylovschur" 
			+ " -eps_target "+sigma
			+ " -st_type sinvert"
			+ " -eps_nev "+nev
			+ " -eps_gen_hermitian",
		 values = EigenVal, vectors = EigenPsi);

if(!mpirank){
	for(int i = 0; i < nev; i++){
		cout << "Eigenvalue["+i+"] = "<<EigenVal[i] 
		<<" , Energy level: " << EigenVal[i]*hbar^2*k^2/(2*me)*kgm2Ds2ToeV << " eV" << endl;
	}
}

The result is [-13.56 eV, -3.38 eV, -3.37 eV ] which is very close to exact enery level of hydrogen atom.

After some experiments, I find that using P2 with a coarser mesh can achieve more accurate result than P1 with a fine mesh:

real hs1= 0.08*R1; 
real hs2=0.1*R2; 
func Pk=P2;

Now the result will be [-13.62, 3.41, 3.40] which is more close to exact result.

I tried P3, but it requires unacceptable memory.

You should not use solver="SPARSESOLVER", that is why you script is currently requiring a lot of memory.