Hi, I’m looking to solve the 3D Schrodinger equation using SLEPc w/ FreeFEM. As part of my code testing I’ve simulated the first 4 eigenvalues/vectors for the “particle in a box” model, with a 10x10x10 box. The eigenvalues I expect are (in units hbar = m =1) 0.049 with 3 fold degeneracy, and 0.197. The states corresponding should be of the form sin(pix/10),sin(piy/10),sin(pi*z/10) and sin(2pix/10). The last state is degenerate with the corresponding y and z states; I’m not sure how FreeFEM would decide which state to pick. However when I attempt to solve for this model I instead find eigenvalues 14.97, 30,1, 30.1, 30.36. Furthermore the eigenstates appear to be incorrect in their shape, and at the very least unnormalised. However, when I remove all boundary conditions I find what appears to be the correct eigenstates for Neumann boundary conditions on the box, correctly normalised, and the correct eigenvalues, with the exception of a spurious mode with a negative eigenvalue ~10^-17, and a corresponding spurious eigenstate. Introducing an isotropic quadratic potential I find all eigenvalues and states to be correct with Neumann BC. My questions are as follows;
- How can I get the correct eigenvalues/eigenvectors for the Dirichlet BC case?
- What is causing the spurious mode in the Neumann case and how can I prevent it?
- In line 198 of my .edp file I attempt to calculate the density from the exact ground state for 3 non-interacting particles in a box w/Neumann conditions. However this integral gives 0, and I’m not sure why?
I’m receiving an error that new users cannot upload files, so please find below a copy of my code;
load “msh3” //Converts mesh files to FreeFem format
load “medit” //.mesh Medit format, use if .msh isnt working
load “mshmet” //Used for mesh truncation/adaptation
load “parmmg” //Used with mshmet
load “gmsh” //Read .msh files for input
load “PETSc” //Library for parallel processing
macro dimension()3// //Tells FreeFem we’re in 3D
include “macro_ddm.idp” //Set of macros for partitioning mesh
macro def(i)i// //Macro for field definition
macro init(i)i//
macro grad(u) [dx(u), dy(u), dz(u)] // //Macro to calculate gradient scalar field
//“C:\Program Files\Microsoft MPI\Bin\mpiexec.exe” -n 16 FreeFem+±mpi.exe AnalyticSchrodinger.edp - v 0 -wg
//string meshFileName = “CubeSmall”; //Testing w/ medit format
//mesh3 Mesh = readmesh3(meshFileName);
//mesh3 Mesh = gmshload3(“CubeSmall.msh”); //Testing w/ gmsh format
//Loads mesh from gmsh .msh file
//DM plex(“CubeSmallv4.msh”, communicator = mpiCommSelf);
//mesh3 Mesh = plex;
mesh3 Mesh = cube(10,10,10);
real vacPerm = 0.008854188; //Vacuum permittivity in farads/nanometre
real elecCharge = 1.6022*(10^(-19)); //Electron charge in Coulombs
real kineticFactor = 3.80998*(10^(-20)); //hbar^2/2m in units eV(nm^2)
real[int] GaAs = [13.1, 0.063, 1.521];
real[int] AlGaAs = [12.11, 0.0904, 1.726];
Mesh = trunc(Mesh, 1, split=2, label = 11);
//Refines mesh by conveting from n nodes to n*split nodes
//Increasing split gives worse eigenvalues
//Increasing cores gives better eigenvalues
if(mpirank == 0)
cout << labels(Mesh) << endl;
mesh3 meshGlobal = Mesh;
//Copy of mesh saved in order to work with functions over the whole finite element space
int[int] n2omesh;
macro MeshN2O()n2omesh//
//Macro keeps global numbering of mesh in array n2omesh
if(mpirank == 0)
cout << labels(Mesh) << endl;
int[int][int] intersection; //Renumbering for local mesh to global
real[int] D; //Partition of unity to convert from local to global
func Pk = P1; //Finite element order for psi
int s = 1; //refinement factor local mesh
build(Mesh,
s,
intersection,
D,
Pk,
mpiCommWorld
) //Builds local distributed mesh
//buildDmesh(Mesh)
//Partitions mesh into -n subdomains with ghost elements for BCs, called local mesh
fespace Vh(Mesh, Pk); //Local finite element space for Psi
fespace Wh(Mesh, P0); //Local finite element space for effective mass
Wh Vtrial = (10^100)*(z >= 5);
//Wh elecMass = (0<=z<=5)(GaAs[1])+(5<z<=70)(AlGaAs[1])+(5<z<=10)(GaAs[1]); //Electron mass as a function of z
//Wh bandgap = (0<=z<=5)(GaAs[2])+(5<z<=7)(AlGaAs[2])+(7<z<=10)(GaAs[2]); //Band gap as a function of z
varf vSchrodinger(psi,v) = int3d(Mesh)(grad(psi)'grad(v)0.5) //variational form Schrodinger eq kineticFactor/elecMass(10^(20))
//+int3d(Mesh)(psiv)//bandgapvpsi(10^(20))
//+int3d(Mesh)(0.5*((x-5)^2+(y-5)^2+(z-5)^2)vpsi); //Potential term
//+on(11, psi=0.0) //Wavefunction vanishes on device surface
//+on(15, psi=0.0) //No boundary conditions gives right eigenvalues
//+on(9, psi=0.0) //Probably 1e-20 giving right BC?
//+on(16, psi=0.0)
//+on(8, psi=0.0)
//+on(13, psi=0.0)
+on(6, psi=0.0)
+on(5, psi=0.0)
+on(4, psi=0.0)
+on(3, psi=0.0)
+on(2, psi=0.0)
+on(1, psi=0.0);
//+on(17, psi=0.0)
//+on(14, psi=0.0);
//+on(7, psi=0.0)
//+on(10, psi=0.0);
varf RHS(psi,v) = int3d(Mesh)(psi*v); //Right hand side variational form
matrix A = vSchrodinger(Vh,Vh, tgv=-2); //FreeFEM matrices for LHS and RHS
matrix B = RHS(Vh,Vh,eps=1e-20);
Mat DistA(A, intersection, D, clean = true); //Ditributed PETSc matrices over processes
Mat DistB(DistA, B, clean = true);
real[int,int] space(DistA.n,1);
space = 1;
real[int] eigenval(0); //Array for eigenvalues
Vh[int] def(eigenvec)(1); //Array for eigenvectors
int nEV = 4;
real relTol = 1e-15;
int maxIt = 20000;
string ssparams = //Parameters for eigenvalue solver
" -eps_nev " + nEV + //Number of eigenvalues nEV
" -eps_type krylovschur" + //Krylov Schur decomposition
" -eps_target "+ 0.0 + //Shift value
" -st_type sinvert " + //Using shift invert method
" -st_pc_type lu " + //Cholesky preconditioner
" -eps_view " + //Outputs eigenvalues to terminal
" -eps_gen_nonhermitian " + //Hermitian matrix
" -eps_tol " + relTol + //Tolerance on egenvalue calculation
" -eps_max_it " + maxIt //Maximum no. iterations for solver
;
int k = EPSSolve(DistA,DistB, //SLEPc Parallel eigenvalue solver
vectors = eigenvec, //Array to store eigenfunctions
values = eigenval, //Array to store eigenvalues
sparams = ssparams
);
k=min(k,nEV); //Ignores converged eigienvalues outside range
Vh Temp;
for(int i=0;i<k;i++){
if(!mpirank) cout << " Eigenvalue #“+i+” = “+eigenval[i]<<endl;
Temp = eigenvec[i];
plotMPI(Mesh, // The local mesh
Temp, // The local solution
Pk, // Local FE-space
def, // Macro for field definition
real, // Type
cmm = “Psi(”+i+”) EV = "+eigenval[i]
)
}
//Eigenvectors are normalised, can test using this loop
Vh densityTotal;
for(int j=0;j<k;j++){
Vh scaledPsi = eigenvec[j];
for[j,v : DistA.D] scaledPsi[j] = v;
densityTotal = densityTotal + scaledPsieigenvec[j];
real localNorm = int3d(Mesh)(scaledPsi*eigenvec[j]); //Calculates square integral
real globalNorm;
mpiAllReduce(localNorm, globalNorm, mpiCommWorld, mpiSUM); //Accounts for ghost element multiplicity in norm
if(mpirank == 0)
cout << globalNorm << endl; //Shows eigenvectors are normalised
}
//for(int j=0;j<k;j++){
//Vh scaledPsi = eigenvec[j];
//for[j,v : DistA.D] scaledPsi[j] = v;
//real localNorm = int3d(Mesh)(scaledPsieigenvec[j]); //Calculates square integral
//real globalNorm;
//mpiAllReduce(localNorm, globalNorm, mpiCommWorld, mpiSUM); //Accounts for ghost element multiplicity in norm
//if(mpirank == 0)
//cout << reduce << endl; //Shows eigenvectors are normalised
//}
for(int j=0;j<k;j++){
Vh scaledPsi = eigenvec[j];
for[j,v : DistA.D] scaledPsi[j] = v;
Vh density = scaledPsieigenvec[j];
plotMPI(Mesh, // The local mesh
density, // The local solution
Pk, // Local FE-space
def, // Macro for field definition
real, // Type
cmm = “|Psi|^2(”+j+") EV = "+eigenval[j]
)
}
plotMPI(Mesh, // The local mesh
densityTotal, // The local solution
Pk, // Local FE-space
def, // Macro for field definition
real, // Type
cmm = “n”
)
Vh exactPlaneWaveDensity3 = ((1/500)^(0.5))((cos(pix/10))^2+(cos(piy/10))^2+(cos(piz/10))^2);
real localElectronNumber = int3d(Mesh)(densityTotal);
real globalElectronNumber;
mpiAllReduce(localElectronNumber, globalElectronNumber, mpiCommWorld, mpiSUM);
if(mpirank == 0)
cout << "globalElectronNumber = " + globalElectronNumber << endl; //Checks density
densityTotal = densityTotal * nEV/globalElectronNumber; //Corrects for compounding floating point error in density
//real correctedLocalElectronNumber = int3d(Mesh)(densityTotal);
//real correctedGlobalElectronNumber;
//mpiAllReduce(correctedLocalElectronNumber, correctedGlobalElectronNumber, mpiCommWorld, mpiSUM);
//if(mpirank == 0)
// cout << "correctedGlobalElectronNumber = " + correctedGlobalElectronNumber << endl; //Checks density has been corrected
real localError = int3d(Mesh)(exactPlaneWaveDensity3 - densityTotal);
real globalError;
mpiAllReduce(localError, globalError, mpiCommWorld, mpiSUM); //Accounts for ghost element multiplicity in norm
if(mpirank == 0)
cout << "globalError = " + globalError << endl; //Shows eigenvectors are normalised
real localExactNumber = int3d(Mesh)(exactPlaneWaveDensity3);
real globalExactNumber;
mpiAllReduce(localExactNumber, globalExactNumber, mpiCommWorld, mpiSUM); //Accounts for ghost element multiplicity in norm
if(mpirank == 0)
cout << "globalExactNumber = " + globalExactNumber << endl; //Shows eigenvectors are normalised
//int[int] fforder(1);
//fforder = 1; //1 for P1 or higher, 0 for P0
//savevtk(“pretty_test_quantum.vtu”, Mesh, densityTotal, order = fforder);
