Chapter 3. Numerical Methods
3.4 Discretization of the Three-dimensional Time independent Schrödinger
for an Yukawa like soft-coulomb potential.
( ) ( ) ( )
It will be a general form for TDSE with single-active-electron assumption.
3.4 Discretization of the Three-dimensional Time independent Schrödinger Equation
For a generalized eigenvalue problem, the QZ [57] algorithm is the most popular
one to solve it. The QZ algorithm code or program can be easily found in many
numerical mathematical libraries, such as LAPACK and IMSL. Unfortunately, the QZ
algorithm is not efficient for large matrix problem. The computational cost is
proportional to cube of matrix size. For a large matrix (>3000), another efficient
algorithm is implemented in this problem. The Jocobi-Davidson algorithm [58] is an
iterative algorithm, which solves the generalized eigenvalue problem iteratively. It was
found that it is very efficient for large matrix problem for the first few selected
eigenstates. The Jocobi-Davidson algorithm is only capable of dealing with the
symmetric matrix problem. If we use the FVM to discretize equation (3-1), equation
(3-1) in matrix form is not a symmetric matrix because of the potential term. In this
thesis, we have applied finite-element-method to discretize equation (3-1) to form a
symmetric matrix, which the Jocobi-Davidson algorithm can handle efficiently.
We use the finite-element method (FEM) to discretize equation (3-1) on structured
non-uniform grid and the details are described as follows. A typical two-dimensional
projected grid is shown in Figure 1. First, the wave function is approximated as
( )
i i( )
( i
( )
1There are eight grid points within a hexahedron. Numbering and coordinating are
shown schematically in Figure 2. We employ the Lagrange polynomials which can
satisfy the general properties of shape function as described earlier. The shape functions
of a hexahedron can be written as follows:
( ) [ ][ ][ ]
a shape functions denotes the local node number in a typical element as shown in Figure2.
By substituting equation (3-5) into equation (3-1), we can obtain:
( )
12 2( ) ( )
i i yukawa i i
i i
E c N r V r c N r
⋅∑ = − ∇ + ∑ (3-7)
Define the Galerkin weighted residual function as
( )
1 2( ) ( )
equation (3-8) for each hexahedron element in the computational domain, we can obtain the following: .Further, by applying the divergence theorem, equation (3-9) is reduced to
( ) ( )
By applying equation (3-11) to each hexahedron in the computational domain, we
can construct a generalized eigenvalue matrix equation into the form as Ax=λB x, which can then be solved numerically. In the present study, initial spatial distribution of
wave function for the time-dependent Schrödinger equation is obtained by numerically
solving the generalized eigenvalue matrix equation as mentioned in the above using the
Jocobi-Davidson algorithm [58].
3.5 Discretization of the Three-dimensional Time Dependent Schrödinger Equation
As introduced in chapter 1 that very few studies have focused on direct real-space
discretization of the TDSE, except those using FDM [33] and FEM by Collin’s group
[33-37]. It is generally agreed that for solving PDEs the memory storage and
computational time by using FEM would be higher as compared to the cell-centered
FVM for achieving the same solution accuracy, as mentioned in chapter 1. In this thesis,
we solve the 3D TDSE directly on real-space coordinates using cell-centered FVM,
which is much simpler in practical implementation and faster in simulation speed as
compared to those using FEM.
By first dividing the volume of interest into several discrete cells and applying the
standard finite-volume method [59] by taking volume integration to the TDSE, equation
(3-2), in each discrete cell, we can obtain
(
,)
3 1 2(
,) (
,)
3theorem, equation (3-11) is reduced to
(
,)
3 1(
,) (
,) (
,)
3(3-12) using Cartesian-grid based non-uniform hexahedral cells (or termed
“non-conformal” mesh), which a typical sketch of mesh projected in two-dimensional
plane is shown in Figure 3. Then, in each cell, equation (3-12) can be simply
respectively. In addition, Ns and ∆VΩ represents number of surfaces and volume of the
cell under consideration. Note the gradient terms at the cell interface can be further approximated by central difference scheme using values of ψ at centroids of
neighboring cells.
The time propagation term on the left-hand side of equation (3-13) is approximated
using an explicit stagger-time algorithm following the idea presented by Visscher [60]
for 1D time-dependent Schrödinger equation, in which the algorithm was shown to be
2nd-order accuracy in time for an uniform grid. The ideas are redescribed here for
( )
the wave function ψ, respectively. By separating these two terms, equation (3-14) can
be further reduced to
Then, we can propagate the equation (3-15a) and equation (3-15b) alternatively in time
using a leap-frog like explicit scheme, termed as “explicit stagger-time scheme”, as the
following:
interpolation of either R or I. In addition, absorbing type boundary conditions [61] are
employed at the outer boundaries of computational domain, which is similar to previous
studies in this regard. Note larger domain size is often necessary to delay the wave
reflection from the numerical boundaries, which indeed is worthwhile of studying in the
future. The present TDSE solver is designed to easily set up the arbitrary number of
regions having different cell sizes centered around the molecule, where most refined
cells are clustered in this region (as shown in Figure 3). In generating the
non-conformal grid, we always set the ratio of cell length between two neigboring cells
at the interface refined regions as two, which further ensures numerical accuracy using
this kind of mesh.