v0.16.3
Loading...
Searching...
No Matches
SCL-13: Schrödinger equation

Table of Contents

Note
Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh


Note
Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • solving eigen problem
  • Use of form integrators
  • How to push the developed UDOs to the Pipeline

Introduction

Schematic of the problem

In this tutorial, we solve a basic quantum-mechanical problem: the Schrödinger equation in a 2D infinite potential well. Below is the interactive visualization of the electron in an infinite potential well.

Electron in an Infinite Potential Well

👈 Drag to adjust n 👉
En = 1 E1

Figure 1: Interactive simulation of the wave function and probability distribution inside the potential well.

Problem Statement

The following equation describes the dynamics of the spatial probability distribution of electrons:

\[ \begin{equation} H(\mathbf{r},t)\psi(\mathbf{r},t) = \mathrm{i}\hbar\frac{\partial}{\partial t}\psi(\mathbf{r},t), \label{eq:time_dep_sch} \end{equation} \]

where \( \hbar \) is the reduced Planck constant. The wave function \( \psi \) describes an electron, and its squared absolute value gives the spatial probability density. The Hamiltonian operator \( H \) is defined as

\[ \begin{equation} H(\mathbf{r},t)=-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{r},t), \label{eq:time_H} \end{equation} \]

where \( m \) is the electron mass. The first term on the right-hand side of Eq. \( \eqref{eq:time_H} \) represents the kinetic energy, and the second term \( V(\mathbf{r},t) \) represents the potential energy from the electrical field. Given a steady state when \( V \) does not evolve with time, \( H \) is purely space-dependent. Consequently, the wave function \( \psi \) can be separated into a purely time-dependent part and a space-dependent part. In the time-dependent phase part, a constant \( E \) is used to describe the rate of change of the wave function, which is the system energy:

\[ \begin{equation} \psi(\mathbf{r},t)=\psi(\mathbf{r})e^{-\mathrm{i}Et/\hbar}. \label{eq:wave_func_sep} \end{equation} \]

Here \(E\) takes the same role as angular frequency in classical mechanics, as in quantum mechanics the energy of the (quasi-)particle is related to angular frequency by the constant \(\hbar\). Replacing \(H(\mathbf{r},t)\) with \(H(\mathbf{r})\) and substituting Eq. \eqref{eq:wave_func_sep} into Eq. \eqref{eq:time_dep_sch}, we get the time-independent Schrödinger equation:

\[ \begin{equation} H(\mathbf{r})\psi(\mathbf{r})=E\psi(\mathbf{r}), \label{eq:time_ind_sch} \end{equation} \]

where \( H(\mathbf{r})=-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{r}) \). Therefore, \( E \) is the eigenvalue of \( H \), and \( \psi \) is the corresponding eigenfunction. Thus, the problem of solving for the wave function reduces to solving an eigenvalue problem.

For a 2D square infinite potential well, the potential on the boundary is infinitely large for electrons to exist there (see the 1D case shown in Figure 1 where the potential is infinite at both sides of the well). Therefore, the probability density of electrons must vanish on the boundary. The strong form is

\[ \begin{align} \label{eq:schro_strong_form_domain} \left(-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{r})\right)\psi(\mathbf{r}) &= E\psi(\mathbf{r}), && \text{in } \Omega, \\[2mm] \label{eq:schro_strong_form_boundary} \psi(\mathbf{r}) &= 0, && \text{on } \Gamma. \end{align} \]

With the strong form given in Eq. \( \eqref{eq:schro_strong_form_domain} \) and Eq. \( \eqref{eq:schro_strong_form_boundary} \), the weak form can be derived:

find \( \psi \in \mathrm{H}^1_0(\Omega) \) such that

\[ \begin{equation} -\frac{\hbar^2}{2m} \int_{\Gamma}\delta\psi\nabla\psi\cdot\vec{n}\,\mathrm{d}\Gamma + \frac{\hbar^2}{2m} \int_{\Omega}\nabla(\delta\psi)\nabla\psi\,\mathrm{d}\Omega + \int_{\Omega}\delta\psi V(\mathbf{r})\psi\,\mathrm{d}\Omega = E\int_{\Omega}\delta\psi\psi\,\mathrm{d}\Omega \quad \forall \delta\psi \in \mathrm{H}^1_0(\Omega). \label{eq:schro_weak} \end{equation} \]

Here \( \mathrm{H}^1_0(\Omega) \) denotes a subspace of \( H^1(\Omega) \) consisting of functions that are set to zero on the Dirichlet boundary. Therefore, the boundary integral in Eq. \( \eqref{eq:schro_weak} \) vanishes. To compare with a case where analytical solutions are available, this tutorial assumes \( V(\mathbf{r})=0 \).

Discrete Form

To get the discrete form, we approximate the test and trial functions by the following basis functions:

\[ \begin{equation} \delta\psi=\sum_{i=1}^{n}N_i\delta\bar{\psi}_i, \end{equation} \]

\[ \begin{equation} \psi=\sum_{j=1}^{n}N_j\bar{\psi}_j. \end{equation} \]

Therefore, the discrete form is

\[ \begin{equation} \sum_{j=1}^{n} \frac{\hbar^2}{2m} \int_{\Omega}\nabla N_i\nabla N_j\,\mathrm{d}\Omega\;\bar{\psi}_{j} = \sum_{j=1}^{n} E\int_{\Omega}N_iN_j\,\mathrm{d}\Omega\;\bar{\psi}_{j}, \quad \forall i=1,2,\ldots,n. \end{equation} \]

By substituting

\[ \begin{equation} H_{ij}=\frac{\hbar^2}{2m} \int_{\Omega}\nabla N_i\nabla N_j\,\mathrm{d}\Omega, \label{eq:H_matrix} \end{equation} \]

\[ \begin{equation} M_{ij}=\int_{\Omega}N_iN_j\,\mathrm{d}\Omega, \label{eq:M_matrix} \end{equation} \]

the problem reduces to the generalized eigenvalue problem

\[ \begin{equation} [H]\bar{\psi}=E[M]\bar{\psi}, \label{eq:matrix_total} \end{equation} \]

where the scalar \( E \) is the eigenvalue.

Implementation

Physical Constants and Parameters

The involved physical constants and parameters are clarified below.

double hbar = 1.054571817e-34; // Reduced Planck Constant [J*s]
double m0 = 9.1093837015e-31; // Free electron mass [kg]
double scale = 1e-9; // Length scaling factor: nanometer → meter [m]
double effMass = 1; // Effective mass in units of m0 [a.u.]
double q = 1.602176634e-19; // Elementary charge [C]
double potential = 0 * q; // Potential energy inside the potential well [J]
effMass
double effMass
Definition schrod_eig.cpp:46
q
double q
Definition schrod_eig.cpp:47
m0
double m0
Definition schrod_eig.cpp:44
hbar
double hbar
[Physical constants and parameters]
Definition schrod_eig.cpp:43
potential
double potential
Definition schrod_eig.cpp:48
scale
double scale
Definition plastic.cpp:123

Boundary Condition

The infinite-potential-well boundary condition is applied by removing the DOFs associated with boundary entities using removeBlockDOFsOnEntities, which is equivalent to enforcing a zero value of the field on the boundary.

MoFEMErrorCode Example::boundaryCondition() {
MoFEMFunctionBegin;
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simple->getProblemName(), "BOUNDARY", std::string("PHI"),
true); // Dirichlet BC. (PHI=0 on the boundaries).
MoFEMFunctionReturn(0);
}
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition dynamic_first_order_con_law.cpp:536
Example::simple
Simple * simple
Definition schrod_eig.cpp:63
Example::mField
MoFEM::Interface & mField
Reference to MoFEM interface.
Definition plastic.cpp:226
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:450
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93

Assembling the Matrices

The \([H]\) and \([M]\) matrices are assembled by OpDomainGradGrad and OpDomainMass respectively according to Eq. \eqref{eq:H_matrix} and Eq. \eqref{eq:M_matrix}. Here pipeline_mng->getDomainLhsFE()->B = H and pipeline_mng->getDomainLhsFE()->A = M are used to store the assembled matrices into \([H]\) and \([M]\), respectively. As SLEPC's EPS solver is not deeply integrated with PETSc as KSP solvers, the two matrices need to be explicitly passed to EPS solver later.

MoFEMErrorCode Example::assembleSystem() {
MoFEMFunctionBegin;
auto *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-effMass", &effMass,
PETSC_NULLPTR);
auto get_kinetic_coe = [&](const double, const double, const double) {
return hbar * hbar /
(2.0 * effMass * m0 *
(scale * scale)); // Coefficient for kinetic energy term
};
auto get_potential = [](const double, const double, const double) {
return potential;
};
auto get_mass_coe = [](const double, const double, const double) {
return 1.0;
};
auto dm = simple->getDM();
CHKERR DMCreateMatrix_MoFEM(dm, H);
M = matDuplicate(H, MAT_SHARE_NONZERO_PATTERN);
auto calculate_Hamiltonian = [&]() {
MoFEMFunctionBegin;
pipeline_mng->getDomainLhsFE().reset();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("PHI", "PHI", get_kinetic_coe));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("PHI", "PHI", get_potential));
auto integration_rule = [](int, int, int approx_order) {
return 2 * (approx_order - 1);
};
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
pipeline_mng->getDomainLhsFE()->B = H;
CHKERR MatZeroEntries(H);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
MoFEMFunctionReturn(0);
};
auto calculate_mass_matrix = [&]() {
MoFEMFunctionBegin;
pipeline_mng->getDomainLhsFE().reset();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("PHI", "PHI", get_mass_coe));
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR MatZeroEntries(M);
pipeline_mng->getDomainLhsFE()->B = M;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
MoFEMFunctionReturn(0);
};
CHKERR calculate_Hamiltonian();
CHKERR calculate_mass_matrix();
MoFEMFunctionReturn(0);
}
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
Definition child_and_parent.cpp:54
integration_rule
auto integration_rule
Definition free_surface.cpp:207
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1188
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition helmholtz.cpp:25
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:162
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition PetscSmartObj.hpp:236
approx_order
static constexpr int approx_order
Definition prism_polynomial_approximation.cpp:14
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition dynamic_first_order_con_law.cpp:639
Example::M
SmartPetscObj< Mat > M
Definition schrod_eig.cpp:73
Example::H
SmartPetscObj< Mat > H
Definition schrod_eig.cpp:74
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:799

Solving the Eigenvalue Problem

An EPS solver is created by EPSCreate(comm, &eps), with \(H\) as the left-hand-side stiffness matrix and \([M]\) as the right-hand-side mass matrix to solve the problem. The option EPS_GHEP specifies we are solving a generalized Hermitian eigenvalue problem. By 'generalized' it means an additional mass matrix on the right-hand side, and by 'Hermitian' it means the left-hand-side matrix \([H]\) is equal to its complex conjugate transpose. EPS_SMALLEST_MAGNITUDE instructs the solver to compute the eigenvalues starting from the smallest. Parameter settings for the solver are printed out for reference during execution. From its we can see the number of iterations taken by the solver to converge, and nev indicates the number of eigenvalues that will be computed within the specified tolerance tol. type indicates the method being used in the solver, which in this tutorial is Krylov-Schur method.

MoFEMErrorCode Example::solveSystem() {
MoFEMFunctionBegin;
auto create_eps = [](MPI_Comm comm) {
EPS eps;
CHKERR EPSCreate(comm, &eps);
return SmartPetscObj<EPS>(eps);
};
auto setup_eps = [&]() {
MoFEMFunctionBegin;
CHKERR EPSSetProblemType(eps, EPS_GHEP);
CHKERR EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE);
CHKERR EPSSetFromOptions(eps);
PetscInt nev = 20;
EPSSetDimensions(eps, nev, PETSC_DEFAULT, PETSC_DEFAULT);
MoFEMFunctionReturn(0);
};
auto print_info = [&]() {
MoFEMFunctionBegin;
ST st;
EPSType type;
PetscReal tol;
PetscInt nev, maxit, its;
// Optional: Get some information from the solver and display it
CHKERR EPSGetIterationNumber(eps, &its);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
" Number of iterations of the method: %d", its);
CHKERR EPSGetST(eps, &st);
CHKERR EPSGetType(eps, &type);
MOFEM_LOG_C("EXAMPLE", Sev::inform, " Solution method: %s", type);
CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
MOFEM_LOG_C("EXAMPLE", Sev::inform, " Number of requested eigenvalues: %d",
nev);
CHKERR EPSGetTolerances(eps, &tol, &maxit);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
" Stopping condition: tol=%.4g, maxit=%d", (double)tol, maxit);
MoFEMFunctionReturn(0);
};
// Create eigensolver context
eps = create_eps(mField.get_comm());
CHKERR EPSSetOperators(eps, H, M);
// Setup EPS
CHKERR setup_eps();
// Solve problem
CHKERR EPSSolve(eps);
// Print info
CHKERR print_info();
MoFEMFunctionReturn(0);
}
type
std::string type
Definition JsonConfigManager.cpp:28
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:319
tol
double tol
Definition mesh_smoothing.cpp:23
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition dynamic_first_order_con_law.cpp:876
Example::eps
SmartPetscObj< EPS > eps
Definition schrod_eig.cpp:75
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0

Postprocessing

The eigenpairs are retrieved by EPSGetEigenpair, where nn indicates the index of the eigenpair, eigr/eigi indicates the real/imaginary part of the eigenvalue, and D receives the real part of the corresponding eigenvector. Since in this tutorial the wave functions are composed of real-valued sine functions, only the real parts of the eigenvectors are stored. Besides, as \([H]\) is Hermitian, the eigenvalues are real as well, so only eigr is stored. The wave function PHI is obtained from D, and the square of the wave function is calculated by OpSquare operator.

MoFEMErrorCode Example::outputResults() {
MoFEMFunctionBegin;
auto *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto phi_ptr = boost::make_shared<VectorDouble>();
auto square_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHI", phi_ptr));
post_proc_fe->getOpPtrVector().push_back(new OpSquare(
square_ptr,
phi_ptr)); // Calculate square of wave function for checking probability density.
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"PHI", phi_ptr}, {"SQUARE", square_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{})
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
auto dm = simple->getDM();
auto D = createDMVector(dm);
PetscInt nev, nconv, n_output = 0;
CHKERR EPSGetDimensions(eps, &nev, PETSC_NULLPTR, PETSC_NULLPTR);
CHKERR EPSGetConverged(eps, &nconv);
n_output = std::min(nconv, nev);
if (nconv < nev) {
MOFEM_LOG_C("EXAMPLE", Sev::warning,
" Only %" PetscInt_FMT " of %" PetscInt_FMT
" requested eigenpairs converged",
nconv, nev);
}
PetscScalar eigr, eigi;
for (PetscInt nn = 0; nn < n_output; nn++) {
CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, D, PETSC_NULLPTR);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
" Eigenpair = %" PetscInt_FMT " Eigen Energy = %.8g eV", nn,
eigr / q); // Convert the unit Joule to eV for output.
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
CHKERR pipeline_mng->loopFiniteElements();
post_proc_fe->writeFile("out_schrod_" +
boost::lexical_cast<std::string>(nn) + ".h5m");
}
MoFEMFunctionReturn(0);
}
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
MoFEM::createDMVector
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
D
double D
Definition initial_diffusion.cpp:45
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:40
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition dynamic_first_order_con_law.cpp:1159
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition PostProcBrokenMeshInMoabBase.hpp:717

Definition of Operator

In this operator we calculate the square of the wave function per Gauss point, which is the probability density of electrons.

struct OpSquare : public ForcesAndSourcesCore::UserDataOperator {
OpSquare(boost::shared_ptr<VectorDouble> squareField,
boost::shared_ptr<VectorDouble> field)
: ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
squareField(squareField), field(field) {}
MoFEMErrorCode doWork(int side, EntityType type,
DataForcesAndSourcesCore::EntData &data) {
MoFEMFunctionBegin;
const size_t nb_gauss_pts = getGaussPts().size2();
squareField->resize(nb_gauss_pts, 0);
auto t_squareField = getFTensor0FromVec(*squareField);
auto t_field = getFTensor0FromVec(*field);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
t_squareField = t_field * t_field;
++t_field;
++t_squareField;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<VectorDouble> squareField;
boost::shared_ptr<VectorDouble> field;
};
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition HookeElement.hpp:75
NOSPACE
@ NOSPACE
Definition definitions.h:83

Results

In quantum mechanics, the analytical solutions for electrons in a 2D square infinite potential well of length of \(L\) are: \begin{equation} \label{eq:E_theo} E_{n_x,n_y}=\frac{\hbar^2 \pi^2}{2mL^2}(n_x^2+n_y^2),\quad n_x,n_y=1,2,3,... \end{equation} \begin{equation} \label{eq:phi_theo} \psi_{n_x,n_y}(x,y)=\frac{2}{L}\text{sin}\left(\frac{n_x\pi x}{L}\right)\text{sin}\left(\frac{n_y\pi y}{L}\right),\quad n_x,n_y=1,2,3,... \end{equation} With \(L=2\) nm, \(E\) takes the form \(0.09401*(n_x^2+n_y^2)\) eV. Each \((n_x,n_y)\) quantum number pair corresponds to a unique quantum state. As indicated by Eq. \eqref{eq:E_theo}, degenerate states exist, which means a single energy level may accommodate electrons in multiple quantum states. In Table 1, the numerical eigenenergies computed by MoFEM are compared with the analytical solutions in Eq. \eqref{eq:E_theo}. The degeneracy describes the number of linearly independent eigenfunctions corresponding to the same eigenvalue, namely the number of quantum states with the same energy, and non-degenerate states have a degeneracy of 1. Similarly, Figure 2 shows the space-resolved wave functions from MoFEM and compares them with the theoretical expressions given in Eq. \eqref{eq:phi_theo}. The numerical results for both eigenenergies and eigenfunctions are consistent with what is predicted by quantum mechanical theory. The reason why only the wave functions corresponding to non-degenerate energy levels are shown in Figure 2 is that exact closed-form analytical expressions are available only in these cases. In degenerate cases, the eigenvectors merely span a subspace, and any orthonormal basis of that subspace is mathematically valid.

Table 1: the 20 smallest eigenenergies
Quantum State (n_x,n_y) Degeneracy Analytical results (eV) MoFEM (eV)
(1,1)10.188020.18802
(1,2),(2,1)20.470040.47005, 0.47005
(2,2)10.752060.75211
(1,3),(3,1)20.940080.94019, 0.94020
(2,3),(3,2)21.222101.22236, 1.22237
(1,4),(4,1)21.598131.59871, 1.59872
(3,3)11.692141.69278
(2,4),(4,2)21.880151.88103, 1.88112
(3,4),(4,3)22.350192.35195, 2.35120
(1,5),(5,1)22.444202.44621, 2.44625
(2,5),(5,2)22.726222.72892, 2.72898
(4,4)13.008243.01188


(a) Wave function for the first eigenmode

(b) Error of the first eigenmode

(d) Wave function for the fourth eigenmode

(e) Error of the fourth eigenmode

(d) Wave function for the eleventh eigenmode

(e) Error of the eleventh eigenmode

(d) Wave function for the twentieth eigenmode

(e) Error of the twentieth eigenmode
Figure 2: The wave functions and their corresponding errors.

Source code

Main Program (cpp)

/**
* \file schrod_eig.cpp
* \example schrod _eig.cpp
*
* Calculate wave functions of Schrödinger equation in 2d problems.
*
*/
#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 2
#endif
#include <MoFEM.hpp>
#include <schrod_eig.hpp>
#undef EPS
#include <slepceps.h>
using namespace MoFEM;
using namespace schrod_eig;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
using DomainEle = FaceElementForcesAndSourcesCore;
};
template <> struct ElementsAndOps<3> {
using DomainEle = VolumeElementForcesAndSourcesCore;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; // Space dimension of problem, mesh
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
PETSC>::BiLinearForm<GAUSS>::OpGradGrad<1, 1, SPACE_DIM>;
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
PETSC>::BiLinearForm<GAUSS>::OpMass<1, 1>;
//! [Physical constants and parameters]
double hbar = 1.054571817e-34; // Reduced Planck Constant [J*s]
double m0 = 9.1093837015e-31; // Free electron mass [kg]
double scale = 1e-9; // Length scaling factor: nanometer → meter [m]
double effMass = 1; // Effective mass in units of m0 [a.u.]
double q = 1.602176634e-19; // Elementary charge [C]
double potential = 0 * q; // Potential energy inside the potential well [J]
//! [Physical constants and parameters]
int order = 2;
static char help[] = "...\n\n";
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEM::Interface &mField;
Simple *simple;
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
SmartPetscObj<Mat> M; // Mass matrix
SmartPetscObj<Mat> H; // Hamiltonian matrix
SmartPetscObj<EPS> eps;
};
//! [Run problem]
MoFEMErrorCode Example::runProblem() {
MoFEMFunctionBegin;
CHKERR readMesh();
CHKERR setupProblem();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
MoFEMFunctionReturn(0);
}
//! [Run problem]
//! [Read mesh]
MoFEMErrorCode Example::readMesh() {
MoFEMFunctionBegin;
simple = mField.getInterface<Simple>();
MOFEM_LOG("EXAMPLE", Sev::inform)
<< "Read mesh for problem in " << EXECUTABLE_DIMENSION;
CHKERR simple->getOptions();
CHKERR simple->loadFile();
MoFEMFunctionReturn(0);
}
//! [Read mesh]
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
MoFEMFunctionBegin;
// Add field
CHKERR simple->addDomainField("PHI", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
CHKERR simple->setFieldOrder("PHI", order);
CHKERR simple->setUp();
MoFEMFunctionReturn(0);
}
//! [Set up problem]
//! [Applying essential BC]
MoFEMErrorCode Example::boundaryCondition() {
MoFEMFunctionBegin;
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
simple->getProblemName(), "BOUNDARY", std::string("PHI"),
true); // Dirichlet BC. (PHI=0 on the boundaries).
MoFEMFunctionReturn(0);
}
//! [Applying essential BC]
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
MoFEMFunctionBegin;
auto *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-effMass", &effMass,
PETSC_NULLPTR);
auto get_kinetic_coe = [&](const double, const double, const double) {
return hbar * hbar /
(2.0 * effMass * m0 *
(scale * scale)); // Coefficient for kinetic energy term
};
auto get_potential = [](const double, const double, const double) {
return potential;
};
auto get_mass_coe = [](const double, const double, const double) {
return 1.0;
};
auto dm = simple->getDM();
CHKERR DMCreateMatrix_MoFEM(dm, H);
M = matDuplicate(H, MAT_SHARE_NONZERO_PATTERN);
auto calculate_Hamiltonian = [&]() {
MoFEMFunctionBegin;
pipeline_mng->getDomainLhsFE().reset();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainGradGrad("PHI", "PHI", get_kinetic_coe));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("PHI", "PHI", get_potential));
auto integration_rule = [](int, int, int approx_order) {
return 2 * (approx_order - 1);
};
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
pipeline_mng->getDomainLhsFE()->B = H;
CHKERR MatZeroEntries(H);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
MoFEMFunctionReturn(0);
};
auto calculate_mass_matrix = [&]() {
MoFEMFunctionBegin;
pipeline_mng->getDomainLhsFE().reset();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pipeline_mng->getOpDomainLhsPipeline(), {H1});
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpDomainMass("PHI", "PHI", get_mass_coe));
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR MatZeroEntries(M);
pipeline_mng->getDomainLhsFE()->B = M;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
MoFEMFunctionReturn(0);
};
CHKERR calculate_Hamiltonian();
CHKERR calculate_mass_matrix();
MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
//! [Solve]
MoFEMErrorCode Example::solveSystem() {
MoFEMFunctionBegin;
auto create_eps = [](MPI_Comm comm) {
EPS eps;
CHKERR EPSCreate(comm, &eps);
return SmartPetscObj<EPS>(eps);
};
auto setup_eps = [&]() {
MoFEMFunctionBegin;
CHKERR EPSSetProblemType(eps, EPS_GHEP);
CHKERR EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE);
CHKERR EPSSetFromOptions(eps);
PetscInt nev = 20;
EPSSetDimensions(eps, nev, PETSC_DEFAULT, PETSC_DEFAULT);
MoFEMFunctionReturn(0);
};
auto print_info = [&]() {
MoFEMFunctionBegin;
ST st;
EPSType type;
PetscReal tol;
PetscInt nev, maxit, its;
// Optional: Get some information from the solver and display it
CHKERR EPSGetIterationNumber(eps, &its);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
" Number of iterations of the method: %d", its);
CHKERR EPSGetST(eps, &st);
CHKERR EPSGetType(eps, &type);
MOFEM_LOG_C("EXAMPLE", Sev::inform, " Solution method: %s", type);
CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
MOFEM_LOG_C("EXAMPLE", Sev::inform, " Number of requested eigenvalues: %d",
nev);
CHKERR EPSGetTolerances(eps, &tol, &maxit);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
" Stopping condition: tol=%.4g, maxit=%d", (double)tol, maxit);
MoFEMFunctionReturn(0);
};
// Create eigensolver context
eps = create_eps(mField.get_comm());
CHKERR EPSSetOperators(eps, H, M);
// Setup EPS
CHKERR setup_eps();
// Solve problem
CHKERR EPSSolve(eps);
// Print info
CHKERR print_info();
MoFEMFunctionReturn(0);
}
//! [Solve]
//! [Postprocess results]
MoFEMErrorCode Example::outputResults() {
MoFEMFunctionBegin;
auto *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto phi_ptr = boost::make_shared<VectorDouble>();
auto square_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("PHI", phi_ptr));
post_proc_fe->getOpPtrVector().push_back(new OpSquare(
square_ptr,
phi_ptr)); // Calculate square of wave function for checking probability density.
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"PHI", phi_ptr}, {"SQUARE", square_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{})
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
auto dm = simple->getDM();
auto D = createDMVector(dm);
PetscInt nev, nconv, n_output = 0;
CHKERR EPSGetDimensions(eps, &nev, PETSC_NULLPTR, PETSC_NULLPTR);
CHKERR EPSGetConverged(eps, &nconv);
n_output = std::min(nconv, nev);
if (nconv < nev) {
MOFEM_LOG_C("EXAMPLE", Sev::warning,
" Only %" PetscInt_FMT " of %" PetscInt_FMT
" requested eigenpairs converged",
nconv, nev);
}
PetscScalar eigr, eigi;
for (PetscInt nn = 0; nn < n_output; nn++) {
CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, D, PETSC_NULLPTR);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
" Eigenpair = %" PetscInt_FMT " Eigen Energy = %.8g eV", nn,
eigr / q); // Convert the unit Joule to eV for output.
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
CHKERR pipeline_mng->loopFiniteElements();
post_proc_fe->writeFile("out_schrod_" +
boost::lexical_cast<std::string>(nn) + ".h5m");
}
MoFEMFunctionReturn(0);
}
//! [Postprocess results]
//! [Check]
MoFEMErrorCode Example::checkResults() {
MoFEMFunctionBegin;
PetscBool test_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test", &test_flg,
PETSC_NULLPTR);
if (test_flg) {
PetscInt nconv;
CHKERR EPSGetConverged(eps, &nconv);
if (nconv < 1) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"No eigenpairs converged");
}
PetscScalar eigr, eigi;
CHKERR EPSGetEigenpair(eps, 0, &eigr, &eigi, PETSC_NULLPTR, PETSC_NULLPTR);
constexpr double regression_value =
0.188; // Check the result for ground energy level [eV]
if (fabs(eigr / q - regression_value) > 1e-3) {
PetscPrintf(PETSC_COMM_WORLD,
"Calculated ground energy: %.4g, expected %.4g\n",
(double)eigr / q, regression_value);
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Regression test faileed; wrong eigen value. Try higher order or "
"finer mesh.");
}
}
MoFEMFunctionReturn(0);
}
//! [Check]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
SlepcInitialize(&argc, &argv, param_file, help);
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
CATCH_ERRORS;
SlepcFinalize();
CHKERR MoFEM::Core::Finalize();
}
help
static char help[]
Definition activate_deactivate_dofs.cpp:13
main
int main()
Definition adol-c_atom.cpp:46
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition child_and_parent.cpp:35
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
H1
@ H1
continuous field
Definition definitions.h:85
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
order
constexpr int order
Definition dg_projection.cpp:38
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::PETSC
@ PETSC
Standard PETSc assembly.
Definition FormsIntegrators.hpp:201
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
MoFEM::BcManager::removeBlockDOFsOnEntities
MoFEMErrorCode removeBlockDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true, bool is_distributed_mesh=true)
Remove DOFs from problem based on block entities.
Definition BcManager.cpp:72
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
schrod_eig
Definition schrod_eig.hpp:14
schrod_eig.hpp
Calculate the square of wave function.
ElementsAndOps
[Operators_definition]
Definition child_and_parent.cpp:19
Example
[Example]
Definition plastic.cpp:216
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition dynamic_first_order_con_law.cpp:463
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition dynamic_first_order_con_law.cpp:1181
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition plastic.cpp:254
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition plastic.cpp:273
MoFEM::BcManager
Boundary condition manager for finite element problem setup.
Definition BcManager.hpp:384
MoFEM::BcScalarMeshsetType
Template specialization for scalar field boundary conditions.
Definition BcManager.hpp:180
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:83
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldValues
Specialization for double precision scalar field values calculation.
Definition UserDataOperators.hpp:157
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:711
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, ScalarDataPtr > DataMapVec
Definition PostProcBrokenMeshInMoabBase.hpp:716
MoFEM::PipelineManager
PipelineManager interface.
Definition PipelineManager.hpp:24
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Get domain left-hand side finite element.
Definition PipelineManager.hpp:837
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:735
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
schrod_eig::OpSquare
[Calculate square of wave function]
Definition schrod_eig.hpp:17
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13

Operator (hpp)

/**
* @file schrod_eig.hpp
* @brief Calculate the square of wave function
* @date 2025-11-14
*
* @copyright Copyright (c) 2025
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
namespace schrod_eig {
//! [Calculate square of wave function]
struct OpSquare : public ForcesAndSourcesCore::UserDataOperator {
OpSquare(boost::shared_ptr<VectorDouble> squareField,
boost::shared_ptr<VectorDouble> field)
: ForcesAndSourcesCore::UserDataOperator(NOSPACE, OPLAST),
squareField(squareField), field(field) {}
MoFEMErrorCode doWork(int side, EntityType type,
DataForcesAndSourcesCore::EntData &data) {
MoFEMFunctionBegin;
const size_t nb_gauss_pts = getGaussPts().size2();
squareField->resize(nb_gauss_pts, 0);
auto t_squareField = getFTensor0FromVec(*squareField);
auto t_field = getFTensor0FromVec(*field);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
t_squareField = t_field * t_field;
++t_field;
++t_squareField;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<VectorDouble> squareField;
boost::shared_ptr<VectorDouble> field;
};
//! [Calculate square of wave function]
} // namespace schrod_eig
schrod_eig::OpSquare::field
boost::shared_ptr< VectorDouble > field
Definition schrod_eig.hpp:45
schrod_eig::OpSquare::squareField
boost::shared_ptr< VectorDouble > squareField
Definition schrod_eig.hpp:44
schrod_eig::OpSquare::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition schrod_eig.hpp:23
Generated by Doxygen 1.9.8 and hosted at University of Glasgow