💥quantumGrid💥

Brought to you by the AMO theory group at Berkeley National Lab.

Note:

The image above is for a scattered wave on a 2D FEM-DVR grid which was created with POV-Ray. The next two figures were created with the example scripts time_dep_two_potential_excitation and two_electron, respectively.

About

quantumGrid is a package for solving a 1-D Schrödinger equation for an arbitrary potential on any interval. The heart of this package is using a Finite Element Method with a Discrete Variable Representation (FEM-DVR) grid to solve the time-dependent or time-independent Schrödinger equation. This grid provides a compact supported foundation for numerically accurate integration and also allows for a natural application of outgoing scattering boundary conditions by adding a complex tail as the last finite element of the FEM-DVR grid, called exterior complex scaling (ECS). Therefore, this grid can be applied to scattering problems where the resonances become square integrable under this complex rotation of the Hamiltonian.

Motivation

This python package was created for a graduate course in time-dependent quantum mechanics at UC Davis. Given the generality and usefulness of a Finite Element Method - Discrete Variable Representation (FEM-DVR) grid for solving the Schrödinger equation and simple scattering problems, we wanted to go open source and provide this numerical tool for others in the name of science!

Contributors ✨

Thanks goes to these wonderful people (emoji key):

Willaim (Bill) McCurdy 💻 🚧 📖	Zachary Streeter💻 🚧 📖	Giuseppe Barbalinardo💻

This project follows the all-contributors specification. Contributions of any kind welcome!

Installation

Recommended

Before going further using conda, you should always update your conda package system:

$ conda update -y conda

First, we always recommend to create a conda environment for using our package. This is the general recommended procedure so there are no dependency issues and systemic issues that could break other parts of your computer. To create a conda environment named “DVRenv” run this command in your terminal:

$ conda create -n DVRenv

Now we want all the rest of the packages to only be installed into this environment so activate it before moving forward:

$ conda activate DVRenv

If you have Anaconda intergrated with your shell, you should see (DVRenv) in the front of your prompt, indicating you are now in the DVRenv environment. If you do not have Anaconda integrated with your shell, then run the following command and confirm you see DVRenv on the next line in your terminal:

$ echo $CONDA_DEFAULT_ENV
$ DVRenv

Now the quantumgrid package is on the PyPI index so we need to install pip to access that index.

(DVRenv) $ conda install pip

(Note that this pip will only be installed in our DVRenv environment!)

Now we can install quantumgrid!

(DVRenv) $ pip install quantumgrid

Now you should be able to use the quantumgrid package in your DVRenv conda environment! You can also run the example scripts simply by executing (see example directory for more details):

(DVRenv) $ ecs_femdvr_time_indep_h2
(DVRenv) $ ecs_femdvr_time_indep_h2 --want_to_plot=true

From sources

The sources for quantumGrid can be downloaded from the Github repo.

You can either clone the public repository:

$ git clone git://github.com/zstreeter/quantumGrid

Or download the tarball:

$ curl -OJL https://github.com/zstreeter/quantumGrid/tarball/master

Once you have a copy of the source, you can install it with:

$ python setup.py install

Background

These notes are an introduction to Discrete Variable Representations (DVRs) using an example that has particularly general applicability. The Finite Element Method with a Discrete Variable Representation (FEM-DVR) provides a way to solve the time-independent or time-dependent Schrödinger equation that, like all DVR methods, is more accurate and faster than finite difference. This method is one of a family of Discrete Variable Representations that are in common use today in chemistry and physics. It can be applied to problems with any potential on any interval. These notes explain the FEM-DVR method using Gauss-Lobatto quadrature, and also outline the Crank-Nicolson propagator for solving the time-dependent Schrödinger equation.

Introduction

Here we describe a method for solving the Schrödinger equation for a particle moving in one dimension with coordinate \(x\) for any potential \(V (x)\) on any interval of \(x\). The variational method of course provides a way to do so, but its application generally poses a practical problem we would like to overcome: If we expand the unknown wave function in \(H |\Psi\rangle = E |\Psi\rangle\) in a finite set of basis functions

\[|\Psi\rangle \approx \sum_{n=1}^N c_n |\varphi_n\rangle\]

substitute it into the Schrödinger equation, and project from the left with \(\langle \varphi_m |\), we come quickly to the familiar matrix representation

\[\begin{split}\mathbf{H}\vec{c} &= E \vec{c} \\ \textrm{with} \quad H_{mn} &= \langle \varphi_m|\hat{T}|\varphi_n \rangle + \langle \varphi_m|\hat{V}|\varphi_n \rangle\end{split}\]

that we can also get from the variational theorem. This is a variational basis representation of the Schrödinger equation.

To construct this matrix eigenvalue problem we need the matrix elements of both the kinetic energy and potential energy operators \(\hat{T}\) and \(\hat{V}\). If we choose a basis for which the kinetic energy matrix elements are easy to evaluate, energy operators T and then try to apply it to solving this problem for various potentials, we generally find that the matrix elements are difficult to perform for many of those potential functions. The DVR is a way of getting around this problem. It is described in an extensive literature on many kinds of DVR that began in the 1980s with seminal work by John Light in the Chemistry Department at the University of Chicago and his coworkers [1]. The central idea is this: We choose a particular form of basis functions, no matter what the potential, that are constructed based on a Gaussian quadrature. Then we use that Gaussian quadrature to do every integral in the problem ­ of both the kinetic and potential energy operators. The result is that the potential energy matrix is diagonal

\[\langle \varphi_m|\hat{V}|\varphi_n\rangle = \delta_{nm}V(x_n)\]

where \(x_n\) is a Gauss quadrature point. In other words, in this basis set method the potential energy matrix is always diagonal, and there are no potential energy integrals to be evaluated. We only have to evaluate the potential on a grid of quadrature points in \(x\). So if we can evaluate the potential energy function as a function of the coordinates, we can reduce the Schrödinger equation directly to a matrix eigenvalue equation. To see how this works we have first to familiarize ourselves with the basics of Gaussian quadrature for performing integrals in general.

Gassian Quadrature

A Gaussian quadrature is a set of points and weights for approximating an integral as a sum. The most basic example is the Gauss-Legendre quadrature summarized together with many other Gauss quadratures in Abramowitz and Stegun [2] (the NBS Handbook of Mathematical Functions) starting on page 887, or in the successor to reference [2], in NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/ chapter 3, section 3.5 on Gauss quadrature (see Eqs. 3.5.15 and 3.5.19)

\[\begin{split}\begin{eqnarray} \int_{-1}^{+1} f(x) dx = \sum_{i=1}^n w_i f(x_i) + R_n \\ R_n \propto f^{(2n)}(\xi) \qquad \textrm{for some} -1 < \xi < 1 \end{eqnarray}\end{split}\]

the points, \(x_i\) , are the zeros of the nth order Legendre polynomial \(P_n(x)\), and the weights, \(w_i\) are given by formulas involving the derivative of that polynomial, in this case, \(w_i = 2/ ((1 - x_i )^2 [P_n (x_i)]^2)\). The error is proportional to the \(2n^{th}\) derivative of the integrand. That means that this quadrature is exact if \(f(x)\) is any polynomial of degree \(2n - 1\) or lower. This property is at the heart of Gauss quadrature and the DVR method. The quadrature can be scaled onto any finite interval

(1)\[\begin{split}\begin{eqnarray} \int_a^b f(y) dy &= \frac{b-a}{2}\sum_{i=1}^n w_i f(y_i) + R_n \\ y_i &= \Big(\frac{b-a}{2}\Big)x_i + \Big(\frac{b-a}{2}\Big) \end{eqnarray}\end{split}\]

by scaling the points and weights. The Legendre polynomials do not have zeros at the endpoints of the interval [-1, 1], so all the quadrature points at any order are interior to that interval. For our DVR we need to have points in the quadrature at the endpoints of the interval, because that is where we have to apply the boundary conditions of the Schrödinger equation. That quadrature is called the Gauss-Lobatto quadrature [2],

\[\begin{split}\begin{eqnarray} \int_{-1}^{+1} f(x) dx = \frac{2}{n(n-1)}[f(-1) + f(1)] + \sum_{i=2}^{n-1} w_i f(x_i) + R_n \\ R_n \propto f^{(2n)}(\xi) \qquad \textrm{for some} -1 < \xi < 1 \end{eqnarray}\end{split}\]

which can of course be written in the form of Eq. (1) with the same scaling. The only difference is that before the scaling the first and last points are explicitly \(-1\) and \(+1\), with the weights \(\frac{2}{n(n -1)}\). There are algorithms for computing the remaining Gauss- Lobatto points and weights at any order, and those are available in C, C++, Fortran and Python libraries. Fixing the two endpoints to be quadrature points lowers the accuracy of the quadrature. It now quadratures polynomials of degree up to \(2n-3\) exactly. We will use the Gauss-Lobatto quadrature here, but there are many other Gauss quadratures. In general they quadrature integrals of the form

\[\int_a^b f(x) W(x) dx = \sum_{i=1}^{n} w_i f(x_i) + R_n\]

with some positive definite weight function \(W(x)\). For example if \(a = -\infty\) and \(b = \infty\), with the weight function \(W = e^{-x^2}\), the quadrature is based on Hermite polynomials and is called Gauss-Hermite quadrature. See Abramowitz and Stegun [2] for a summary of eight common quadratures. Of course, there is a Wikipedia page too.

DISCRETE VARIABLE REPRESENTATION OF THE WAVE FUNCTION AND HAMILTONIAN

The Gauss-Lobatto quadrature points provide a way to construct \(n\) polynomials of degree \(n - 1\) called “interpolating polynomials” that are each zero at all but one of the quadrature points,

(2)\[L_j(x) = \prod_{i \ne j}^n \frac{x-x_i}{x_j-x_i}\]

and are equal to 1 at \(x = x_j\) . So if we have the values of a function at the quadrature points we can interpolate it by effectively fitting it with Nth order polynomials as

\[f(x) \approx \sum_{j=1}^n f(x_j) L_j(x)\]

this fit reproduces the function exactly at the quadrature points, and provides a polynomial interpolation between them. This kind of interpolation, using interpolating poynomials of the form of Eq. (2), the starting idea of the DVR and also the start of the derivation of quadrature rules like Simpson’s rule that involve evenly (or unevenly) spaced quadrature points. To form our DVR basis functions, \(\phi_j(x)\), we know normalize these interpolating polynomials by multiplying by \(1/\sqrt{w_j}\),

(3)\[\phi_j(x) = \frac{1}{w_j^{1/2}}\prod_{i \ne j}^n \frac{x-x_i}{x_j-x_i}\]

These functions are now normalized in the very specific sense that they are orthonormal within the quadrature approximation

\[\int_a^b \phi_i(x) \phi_j(x) dx \approx \sum_{k=1}^n w_k \phi_i(x_k) \phi_j(x_k) = \delta_{ij}\]

The central idea of a DVR using the Gauss-Lobatto quadrature is to

1. Expand our unknown wave function in these Gauss-Lobatto basis functions and then

2. Define every matrix element as its approximation by the underlying Gauss-Lobatto quadrature

The idea of using the Gauss-Lobatto quadrature to define a DVR was introduced by Manolopoulos and Wyatt [3, 4] in 1989. It allows us to apply boundary conditions at the two endponts, a and b of the interval on which the quadrature is defined. We will apply the boundary conditions that the wave function is zero at the endpoints, which we can accomplish by simply excluding from the basis the basis functions that are nonzero at first and last of the quadrature points. We expand the wave function in the DVR basis

\[\Psi(x) = \sum_{m=2}^{n-1} \psi_m\phi_m(x)\]

where the unknown coefficients are related to the values of the wave function \(\Psi(x)\) at the grid points by \(\psi_m = \Psi(x_m)w_m^{1/2}\). Then the matrix of the kinetic energy is

(4)\[\begin{split}\begin{align} T_{i,j} &= -\frac{\hbar^2}{2m}\int_a^b\phi_i(x)\frac{d^2}{dx^2}\phi_j(x)dx \\ &= \frac{\hbar^2}{2m} \int_a^b \frac{d}{dx}\phi_i(x) \frac{d}{dx}\phi_j(x)dx \\ &= \sum_{k=2}^{n-1}w_k\phi_i'(x_k)\phi_j'(x_k) \end{align}\end{split}\]

where we integrated by parts and use the boundary condition we are applying to get an obviously symmetric matrix result. We can see that the integrand in Eq. (4) is of order \(2(n - 2) = 2n - 4\). That’s because the order of the interpolating polynomials in Eq. (2) is \(n-1\), and so their derivatives are of order \(n-2\). Thus the Gauss-Lobatto quadrature, which integrates polynomials of order \(2n-3\), is exact for the kinetic energy.

For the potential energy we have

\[\begin{split}V_{i,j} &= -\frac{\hbar^2}{2m}\int_a^b\phi_i(x)V(x)\phi_j(x)dx \\ &\approx \sum_{k=2}^{n-1} w_k \phi_i(x_k)V(x_k)\phi_j(x_k) = \delta_{i,j}V(x_i)\end{split}\]

The potential matrix is diagonal because the DVR basis functions are at all the quadrature points but one. That has to be the same point for the two basis functions, or else the integral is zero. With this DVR, the time-independent Schrödinger equation has been reduced to a matrix equation

\[\begin{split}\mathbf{H}\vec{\psi} &= E\vec{\psi} \\ \textrm{with} \quad H_{ij} &= T_{ij} + V_{ii} \delta_{ij}\end{split}\]

The kinetic energy matrix is full and is treated exactly. The potential energy matrix elements are approximated by the quadrature. To improve the approximation we increase the order of the quadrature. We only require the values of the potential at the quadrature points, and once we have the coefficients we can scale them so that the function to be normalized.

DVRs are a popular numerical approach in one or more dimensions. The sort of DVR described here can be used for the angular degrees of freedom in spherical polar coordinates, for example, with the change of variable x = cos . To increase the accuracy of the DVR basis set representation of the wave function one only needs to increase the order of the underlying quadrature.

For large intervals in x and wave functions with short wavelengths, we might need many hundreds or thousands of grid points. In this approach, that would mean calculating the points and weights to high precision for a quadrature that might of the order of thousands. To escape that problem, there is a popular variant of this method called the FEM-DVR.

FINITE-ELEMENT METHOD WITH DISCRETE VARIABLE REPRESENTATION

In this method, first proposed by Rescigno and McCurdy [5] in 2000, we divide up the interval in x into sub-intervals call finite elements. With each of those elements we associate a separate Gauss-Lobatto quadrature and a separate DVR basis constructed according to Eq. (3). However, we must pay special attention to the boundaries between the elements. To allow the wave function to be continuous across the finite element boundaries we must include the Lobatto shape functions that we eliminated above, and that are nonzero at the ends of the finite element interval. To do so, we start by labeling the Lobatto-shape functions according to which interval, \(\alpha\), they are associated with

(5)\[\begin{split}L_j^\alpha(x) = \begin{cases} \prod_{i \ne j}^n \frac{x-x_i^\alpha}{x_j^\alpha-x_i^\alpha} & x_1^\alpha \le x \le x_n^\alpha \\ 0 & x < x_1^\alpha, \quad x_n^\alpha < x \end{cases}\end{split}\]

and similarly labeling the DVR basis functions

(6)\[\phi_j^\alpha(x) = \frac{1}{(w_j^\alpha)^{1/2}} L_j^\alpha(x) \qquad 2 \le j \le n-1\]

We can then create bridging functions that connect the elements \(\alpha\) and \(\alpha+1\) according to

(7)\[\phi_n^\alpha(x) = \frac{1}{\sqrt{w_n^\alpha + w_1^{\alpha+1}}}\left(L_n^\alpha(x) + L_1^{\alpha+1}(x) \right)\]

The bridging functions are normalized, because they are integrated by the combination of the quadratures in the two finite elements.

This idea is illustrated in the figures directly above. The bridging functions extend across two elements, while the FEM-DVR basis within each element is nonzero only on that interval.

We can extend this process to create as many elements as we want. So for example, we might use 20th order quadrature in each element and have 100 elements of various sizes cover the entirety of the interval \(a \le x \le b\).

Rescigno and McCurdy~cite{Rescigno_McCurdy2000} give the formulas for the kinetic energy matrix elements for this method, which are still exactly evaluated by the Gauss-Lobatto quadrature(s) over the entire domain of \(x\). The potential matrix elements remain diagonal. The kinetic energy is no longer a full matrix, but has the form of overlapping blocks originating from their finite elements. Schematically we can display the form of the Hamiltonian for four elements as

\[\begin{split}\begin{eqnarray} \\ \mathbf{H} = &\mathbf{T} &+\qquad \qquad \qquad \mathbf{V} \\ \\ \mathbf{H} =& \begin{pmatrix} X&X&0&0&0&0&0&0&0 \\ X&X&X&X&0&0&0&0&0 \\ 0& X&X&X&0&0&0&0&0\\ 0& X&X&X&X&X&0&0&0\\ 0& 0&0&X&X&X&0&0&0\\ 0& 0&0&X&X&X&X&X&0\\ 0& 0&0&0&0&X&X&X&0\\ 0& 0&0&0&0&X&X&X&X\\ 0& 0&0&0&0&0&0&X&X\\ \end{pmatrix} &+ \begin{pmatrix} X&0&0&0&0&0&0&0&0 \\ 0&X&0&0&0&0&0&0&0 \\ 0&0&X&0&0&0&0&0&0 \\ 0&0&0&X&0&0&0&0&0 \\ 0&0&0&0&X&0&0&0&0 \\ 0&0&0&0&0&X&0&0&0 \\ 0&0&0&0&0&0&X&0&0 \\ 0&0&0&0&0&0&0&X&0 \\ 0&0&0&0&0&0&0&0&X \\ \end{pmatrix} \end{eqnarray}\end{split}\]

where there are four blocks in the kinetic energy matrix, the first and last of which are smaller than all the intermediate blocks, because the lack the first and last basis functions that would have been nonzero at the boundaries of the entire interval in \(x\).

The FEM-DVR basis is orthogonal within the overall Gauss-Lobatto quadrature. The kinetic energy matrix is exactly evaluated with the quadrature. The potential energy is only required at the quadrature points,

\[\begin{equation} \langle \phi_i |V|\phi_j \rangle = \delta_{ij}V(x_i) \end{equation}\]

The DVR representation not variational because both the potential and overlap matrix elements are approximated by the quadrature and are not exact, although with moderately dense grids that property barely noticeable. Also, the FEM-DVR enforces only continuity at the finite element boundaries and not continuity of the derivatives of the wave function. Therefore there is first order error at each finite element boundary. Again, with moderately dense grids the discontinuity in the derivative at the boundaries is very slight and causes no numerical pathologies.

Higher dimensions: Extension of both the DVR and FEM-DVR methods to higher dimensions is straightforward. Suppose we had a problem with \(x\) and \(y\) degrees of freedom. In a product basis of separate DVR basis functions \(\phi_i(x)\) and \(\chi_m(y)\) the kinetic and potential energy matrices are

\[\begin{split}\begin{equation} \begin{split} T_{im,jn} &=\langle \phi_i(x) \chi_m(y) |T_x + T_y| \phi_j(x) \chi_n(y) \rangle = \delta_{mn}T^x_{ij} +\delta_{ij}T^y_{mn} \\ V_{im,jn} &= \langle \phi_i(x) \chi_m(y) |V(x,y) | \phi_j(x) \chi_n(y) \rangle = \delta_{ij} \delta_{mn} V(x_i,y_m) \end{split} \end{equation}\end{split}\]

and extentions to higher dimensions are similar. One only needs the potential evaluated at the quadrature points, and the kinetic energy matrix becomes increasingly sparse as the number of dimensions increases. Scalable algorithms for implementing the FEM-DVR in higher dimensions take advantage of that sparsity, both in storage and numerical operations.

DVRs constructed from other orthogonal polynomials and the notion of a “proper DVR”

Orthogonal polynomials tridiagonalize the matrix of the position with the weight with which the orthogonal polynomials are orthogonal

\[\begin{equation} Q_{i,j} = \int_a^b \phi_i(x) \, x \, \phi_j(x) \, W(x) dx = \textrm{tridiagonal matrix} \end{equation}\]

That is an essential property that is derivable from the three-term recursion relations satisfied by all orthogonal polynomials. Every set of orthogonal polynomials defines a quadrature for their values of \(a\), \(b\) and \(W(x)\). For example we can construct a DVR for the interval \(-\infty < x < \infty\) using harmonic oscillator functions that apply the proper bound-state boundary conditions at infinity, and the underlying quadrature will be the Gauss-Hermite quadrature as described in the appendices of reference [6].

If we begin with a real orthogonal basis \(\{\phi_i(x)\}_{i=1}^{N}\) that tridiagonalizes the position operator, \(x\), we can generate a “proper DVR”, in the terminology of the appendices of the review article [6] by Hans Dieter-Meyer and coworkers, as follows.

1. Diagonalize the matrix of the position operator \(\mathbf{Q}\)

\[\begin{split}\begin{equation} Q_{i,j} = \langle \phi_i|x|\phi_j\rangle \\ \mathbf{Q}=\mathbf{U}\mathbf{X}\mathbf{U}^T \end{equation}\end{split}\]

where \(\mathbf{X}\) is the matrix of position eigenvalues \(X_{\alpha,\beta} = x_\alpha \delta_{\alpha,\beta}\)

2. Transform the exactly evaluated matrix of the kinetic energy, \(\mathbf{T}\) to the basis of position eigenfunctions

\[\begin{equation} T^{DVR}=\mathbf{U}^T \mathbf{T} \mathbf{U} \end{equation}\]

3. Construct the DVR representation of the Hamiltonian as

\[\begin{equation} H^{DVR}_{\alpha,\beta} = T^{DVR}_{\alpha,\beta}+V(x_{\alpha})\delta_{\alpha,\beta} \end{equation} \label{eq:DVRham}\]

in which the potential is diagonal, so no matrix elements of it are required.

The DVR basis in which the Hamiltonian is represented by \(H^{DVR}\) is

\[\begin{equation} \chi_\alpha(x) = \sum_{j=1}^N \phi_j(x)U_{j,\alpha} \end{equation}\]

These functions also provide the “interpolating basis” mentioned below.

Because the original basis tridiagonalized the position operator, we can construct the underlying Gauss quadrature, whose abscissas are \(\{x_\alpha\}_{\alpha=1}^N\), and whose weights are given by

\[\begin{equation} w_\alpha^{1/2}=U_{k,\alpha} / \phi_k(x_\alpha) \end{equation}\]

which in this case does not depend on \(k\). If the original basis had not tridiagonalized \(x\) this formula for the weights would have depended on \(k\).

In a proper DVR we can use the properties of the underlying Gauss quadrature to show that the original basis satisfied discrete orthonormality

\[\begin{equation} \sum_{\alpha = 1}^N w_\alpha \phi_j (x_\alpha) \phi_k(x_\alpha) = \delta_{j,k} \end{equation}\]

and discrete completeness

\[\begin{equation} \sum_{j = 1}^N (w_\alpha w_\beta)^{1/2} \phi_j (x_\alpha) \phi_j(x_\beta) = \delta_{\alpha,\beta} \end{equation}\]

The combination of these two properties means that the the DVR and an expansion of the wave function in the original orthogonal polynomials are exactly equivalent.

Moreover the DVR basis in a proper DVR has the important discrete \(\delta\) -property

\[\begin{equation} \chi_\alpha(x_\beta) = w_\alpha^{-1/2}\delta_{\alpha,\beta} \end{equation}\]

This relation tells us that for any function spanned by the original basis we have, using the quadrature,

\[\begin{equation} \langle \chi_\alpha|\psi \rangle = w_\alpha^{1/2} \psi(x_\alpha) \end{equation}\]

So that if we find the eigenvectors of \(\mathbf{H}^{DVR}\) or use it in the DVR representation of linear equations like \((E-\mathbf{H}^{DVR})\boldsymbol{\psi} = \mathbf{d}\), the resulting vectors represent the values of \(\psi\) on the gridpoints

\[\begin{equation} \psi(x) \rightarrow \boldsymbol{\psi}=(w_1^{1/2}\psi(x_1),w_2^{1/2}\psi(x_2),\cdots,w_N^{1/2}\psi(x_N)) \end{equation}\]

Equally important is the fact that the DVR basis provides the interpolating basis whereby we can get \(\psi\) at any value of $x$ via its expansion

\[\begin{equation} \psi(x) = \sum_{\alpha=1}^N \psi_\alpha \chi_\alpha(x) \end{equation}\]

Note: If the original basis does not tridiagonalize the position operator, but instead tridiagonalizes an invertable function of \(x\), like \(1/x\) or \((x_0-x)^2\), a proper DVR can be generated by diagonalizing that operator, and writing the abscissas \(x_\alpha\) in terms of the eigenvalues of that operator.

Time propagation using the FEM-DVR and Crank-Nicolson propagator

Finite difference in time as a route to approximate solutions of the time-dependent Schrödinger equation

The time-dependent Schrödinger equation has been converted by our DVR to a time-dependent matrix equation (\(\hbar = 1\)),

\[\begin{equation} i\frac{\partial}{\partial t} \vec{\psi}(t) = \mathbf{H}(t) \, \vec{\psi}(t) \end{equation}\]

which is of course an initial value problem. The simplest way to treat it is by approximating the derivative by finite difference,

\[\begin{split}\begin{equation} \begin{split} i \frac{\vec{\psi}(t+\Delta) - \vec{\psi}(t)}{\Delta} \approx \mathbf{H}(t)\, \vec{\psi}(t) \\ \textrm{or}\quad \vec{\psi}(t+\Delta) = \left(1-i\Delta \mathbf{H}(t)\right) \vec{\psi}(t) \end{split} \end{equation}\end{split}\]

which would allow us to step forward in time with a single matrix multiplication. Unfortunately this idea, known as “forward Euler”, produces an unstable algorithm, whose error, depending on the spectrum of \(\mathbf{H}\) (its eigenvalues) grows exponentially. There are many ways to remedy this problem at the cost of more matrix multiplications per step, and several of them can be derived by using higher order finite difference approximations to first derivative involving more points in time. Standard methods, like the variable time-step Runge-Kutta method, and various “predictor-corrector” methods can be found in the literature and software libraries.

A propagator that has the particularly important advantage for quantum mechanics of being unitary is due to Crank and Nicolson~ [7] (whose original paper was about heat-conduction equations). We can ” derive” it easily by replacing the right hand side by an average of it evaluated at \(t\) and \(t+\Delta\)

\[\begin{split}\begin{equation} \begin{split} i \frac{\vec{\psi}(t+\Delta) - \vec{\psi}(t)}{\Delta} & \approx \frac{1}{2} \left( \mathbf{H}(t+\Delta)\, \vec{\psi}(t+\Delta) +\mathbf{H}(t)\, \vec{\psi}(t))\right) \\ &\approx \frac{1}{2} \left( \mathbf{H}(t+\Delta/2)\, \vec{\psi}(t+\Delta) +\mathbf{H}(t+\Delta/2)\, \vec{\psi}(t))\right) \\ \textrm{or}\quad \left(1 + i\frac{\Delta}{2} \mathbf{H}(t+\Delta/2)\right) & \vec{\psi}(t+\Delta) = \left(1 - i \frac{\Delta}{2} \mathbf{H}(t+\Delta/2)\right) \vec{\psi}(t) \end{split} \end{equation}\end{split}\]

We can write the Crank-Nicolson propagator result in two equivalent forms,

(8)\[\begin{split}\begin{equation} \boxed{ \begin{split} \left(1 + i\frac{\Delta}{2} \mathbf{H}(t+\Delta/2)\right) \vec{\psi}(t+\Delta) &= \left(1 - i \frac{\Delta}{2} \mathbf{H}(t+\Delta/2)\right) \vec{\psi}(t) \\ \vec{\psi}(t+\Delta) = \left(1 + i\frac{\Delta}{2} \mathbf{H}(t+\Delta/2)\right)^{-1} & \left(1 - i \frac{\Delta}{2} \mathbf{H}(t+\Delta/2)\right) \vec{\psi}(t) \\ \end{split} } \end{equation}\end{split}\]

Implementing the first relation in Eq. (8) to take a time step requires one matrix multiplication plus a single solution of linear equations. It is one of a class of “implicit” methods that require a solution of linear equations in contrast to “explicit” methods that require only matrix-vector multiplications. Solving the linear equations is the most efficient implementation of the Crank-Nicolson propagator if the Hamiltonian is time dependent, because solving one set of linear equations requires much less computation than inverting a matrix. If the Hamiltonian is not time dependent however we can construct the inverse once and for all and use it to take any number of time steps, which then only require one matrix-vector multiplication each.

(9)\[\begin{split}\begin{equation} \begin{split} \mathbf{U} = \left(1 + i\frac{\Delta}{2} \mathbf{H}\right)^{-1} & \left(1 - i \frac{\Delta}{2} \mathbf{H}\right) \\ \vec{\psi}(t+\Delta) &= \mathbf{U}\vec{\psi}(t) \\ \vec{\psi}(t+n \Delta) &= \mathbf{U}^n\vec{\psi}(t) \quad \textrm{for }\mathbf{H} \, \textrm{ time-independent} \end{split} \end{equation}\end{split}\]

A single matrix inversion and matrix-matrix multiplication then allows us to take $n$ time steps with only n matrix-vector multiplications.

Properties of the Crank-Nicolson propagator

1. Unitarity: Crank-Nicolson propagation is unitary, meaning that for each time step the norm of the wave function is conserved. The reason is that the matrix \(\mathbf{U}\) is unitary, including the case that \(H\) is time-dependent. Proof uses the fact that \(H\) is hermitian, \(H^\dagger = H\):

\[\begin{split}\begin{equation} \begin{split} \mathbf{U}&=\left(1+i\frac{\Delta}{2}\mathbf{H}(t+\Delta/2)\right)^{-1} \left(1- i\frac{\Delta}{2}\mathbf{H}(t+\Delta/2)\right) \\ \mathbf{U}^\dagger& = \left(1+ i\frac{\Delta}{2}\mathbf{H}^\dagger \right) \left(1-i\frac{\Delta}{2}\mathbf{H}^\dagger\right)^{-1} = \left(1+ i\frac{\Delta}{2}\mathbf{H} \right) \left(1-i\frac{\Delta}{2}\mathbf{H} \right)^{-1} \\ \mathbf{U} \mathbf{U}^\dagger =& \left(1+i\frac{\Delta}{2}\mathbf{H}\right)^{-1} \left(1- i\frac{\Delta}{2}\mathbf{H}\right) \left(1+ i\frac{\Delta}{2}\mathbf{H}\right) \left(1-i\frac{\Delta}{2}\mathbf{H}\right)^{-1} \\ =& \left(1+i\frac{\Delta}{2}\mathbf{H}\right)^{-1} \left(1+ i\frac{\Delta}{2}\mathbf{H}\right) \left(1- i\frac{\Delta}{2}\mathbf{H}\right) \left(1-i\frac{\Delta}{2}\mathbf{H}\right)^{-1} \textrm{H commutes with itself and 1} \\ & = 1 \end{split} \end{equation}\end{split}\]

So \(\vec{\psi} (t+\Delta)^\dagger \cdot \vec{\psi}(t+\Delta)= \vec{\psi}(t)^\dagger \mathbf{U}^\dagger \cdot \mathbf{U} \vec{\psi}(t) = \vec{\psi}(t)^\dagger \cdot \vec{\psi}(t)\) and the norm of the wave function is preserved at every step.

2. Error is third order: The order of the error in the Crank-Nicolson propagator is \(\Delta^3\), as we see by expanding in powers of \(\Delta\) and comparing with the expansion of \(\exp(-i\mathbf{H}\Delta)\),

\[\begin{split}\begin{equation} \begin{split} \mathbf{U}&=\left(1+i\frac{\Delta}{2}\mathbf{H}\right)^{-1} \left(1- i\frac{\Delta}{2}\mathbf{H}\right) \\ &= \left(1 - i\frac{\Delta}{2}\mathbf{H} +(i\frac{\Delta}{2}\mathbf{H})^2 -(i\frac{\Delta}{2}\mathbf{H})^3\cdots\right) \left(1- i\frac{\Delta}{2}\mathbf{H}\right) \\ &= 1 - i\mathbf{H}\Delta +\frac{1}{2} \left(- i\mathbf{H}\Delta \right)^2 + O(\Delta^3) \\ &= e^{-i\mathbf{H} \Delta} + O(\Delta^3) \end{split} \end{equation}\end{split}\]

3. Stability: The Crank-Nicolson propagator is unconditionally stable, meaning that no matter what the step size and spectrum of \(\mathbf{H}(t)\) the error of the propagation does not increase exponentially as for forward Euler.

References

1

J V. Lill, Gregory Parker, and John Light. The discrete variable–finite basis approach to quantum scattering. The Journal of Chemical Physics, 85:900–910, 07 1986. doi:10.1063/1.451245.

2(1,2,3,4)

Milton Abramowitz and Irene A. Stegun. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York, ninth dover printing, tenth gpo printing edition, 1964.

3

D.E. Manolopoulos and R.E. Wyatt. Quantum scattering via the log derivative version of the kohn variational principle. Chemical Physics Letters, 152(1):23 – 32, 1988. doi:https://doi.org/10.1016/0009-2614(88)87322-6.

4

David E. Manolopoulos, Michael D'Mello, and Robert E. Wyatt. Quantum reactive scattering via the log derivative version of the kohn variational principle: general theory for bimolecular chemical reactions. The Journal of Chemical Physics, 91(10):6096–6102, 1989. doi:10.1063/1.457428.

5

T. N. Rescigno and C. W. McCurdy. Numerical grid methods for quantum-mechanical scattering problems. Phys. Rev. A, 62:032706, Aug 2000. doi:10.1103/PhysRevA.62.032706.

6(1,2)

M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer. The multiconfiguration time-dependent Hartree method: A highly efficient algorithm for propagating wavepackets. prep, 324:1–105, 2000.

7

J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Advances in Computational Mathematics, 6(1):207–226, Dec 1996. doi:10.1007/BF02127704.

Usage

For examples on how to use quantumGrid, please refer to the examples directory.

For time propagation, we recommend using Intel optimized libraries. Using the Crank-Nicolson propagator does enforce unitarity but is relatively slow.

quantumgrid

quantumgrid package

Submodules

quantumgrid.cli module

Console cript for quantumgrid.

quantumgrid.femdvr module

Finite Eement Method – Discrete Variable Representation (DVR) for 1D Schroedinger equation using Gauss-Lobatto quadrature C. William McCurdy, Zachary Streeter, and Giuseppe Barbalinardo – UC Davis

Includes time propagation with Crank-Nicolson propagator and routines for dynamics on two coupled potential curves.

November 2019

time-dependent potentials allowed for two-state propagation vectorized logic in Hamiltonian and Potential builds

March 2020

Exterior Complex Scaling implemented for single state portion both time-independent and time-dependent calculations Vectorized logic for diagonals of Hamiltonian can be commented and alternate logic uncommented for potential functions that don’t vectorize correctly. 2-state routines are not yet implemented for ECS

class quantumgrid.femdvr.FEM_DVR(n_order, FEM_boundaries, Mass=1, Complex_scale=1, R0_scale=0.0)

Bases: object

Constructor method Builds 1D FEM DVR grid and kinetic energy matrix representation in normalized FEM DVR basis of shape functions and bridging functions Finite elements of any sizes determined by FEM_boundaries array, but all with same order DVR

Parameters

• n_order (int) – The DVR order. More dense takes longer but improves acuracy

• FEM_boundaries (ndarray) – An array of ‘double’ or ‘complex’ boundaries for the Finite-Elements.

• Complex_scale (int,optional) – Flag that starts complex scaling at grid boundary closest to R0 if .ne. 1. Defaults to 1

• R0_scale (complex,optional) – Grid point where the complex tail starts. Defaults to 0.0

n_order

Can access the DVR order later

Type

int

FEM_boundaries

Can access the ‘double’ or ‘complex’ Finite-Element boundaries later

Type

ndarray

nbas

Total number of basis functions in the FEM-DVR grid

Type

int

x_pts

Array of ‘double’ or ‘complex’ DVR points

Type

ndarray

w_pts

Array of ‘double’ or ‘complex’ DVR weight

Type

ndarray

KE_mat

Kintetic Energy matrix, dimension is nbas x nbas of ‘double’ or ‘complex’ types

Type

ndarray

i_elem_scale

index of last real DVR point

Type

int

Crank_Nicolson(t_initial, t_final, N_times, Coefs_at_t_initial, potential, Is_H_time_dependent=True)

Crank Nicolson Propagator for time-dependent or time-independent Hamiltonian Both implementations are unitary. A time step is

\[C_t = \Big(1 + i H\big(t-\frac{\Delta t}{2}\big) *\frac{\Delta t}{2}\Big)^{-1} * \Big(1 - i H\big(t-\frac{\Delta t}{2}\big) *\frac{\Delta t}{2}\Big)*C_t-\Delta t\]

The ECS FEM-DVR Hamiltonian is complex symmetric, and is not conjugated in this expression

Parameters

• t_initial (int) – Initial time for Crank_Nicolson propagation

• t_final (int) – Final time for Crank_Nicolson propagation

• N_times (int) – Number of time steps

• Coefs_at_t_initial (ndarray) – Array of coefficients representing the wavefunction in the FEM-DVR grid at initial time

• potential (function) – A caller provided potential that provides the pertubation to the system

• Is_H_time_dependent (bool) – Boolean for specifying if Hamiltonian is

• time-dependent –

• to true (defaults) –

Returns

new array of coefficients after one Crank Nicolson time-step

Return type

Ct (ndarray)

Crank_Nicolson_Two_States(t_initial, t_final, N_times, Coefs_at_t_initial, V_potential_1, V_potential_2, V_coupling, Is_H_time_dependent=True)

Crank Nicolson Propagator for time-dependent or time-independent Hamiltonian Both implementations are unitary. A time step is

\[C_t = \Big(1 + i H\big(t-\frac{\Delta t}{2}\big) *\frac{\Delta t}{2}\Big)^{-1} * \Big(1 - i H\big(t-\frac{\Delta t}{2}\big) *\frac{\Delta t}{2}\Big)*C_t-\Delta t\]

In this routine the wave function has two components propagating on coupled potential curves specified by V_potential_1, V_potential_2 and V_coupling

Time-dependent case can slow. Building the Hamiltonian and M matrices dominates for small grids, but linear equation solve dominates time for large grids. This implementation allows both the diagonal and coupling potentials to be time-dependent and allows the coupling to be complex hermitian.

Parameters

• t_initial (int) – Initial time for Crank_Nicolson propagation

• t_final (int) – Final time for Crank_Nicolson propagation

• N_times (int) – Number of time steps

• Coefs_at_t_initial (ndarray) – Array of coefficients representing the wavefunction in the FEM-DVR grid at initial time

• V_potential_1 (function) – A caller provided potential for the propagation of the first component of the wavefunction

• V_potential_2 (function) – A caller provided potential for the propagation of the second component of the wavefunction

• V_coupling (function) – A caller provided potential that provides coupling between the two potential curves

• Is_H_time_dependent (bool) – Boolean for specifying if Hamiltonian is time-dependent, defaults to true

Returns

new array of coefficients after one Crank Nicolson time-step

Return type

Ct (ndarray)

Hamiltonian(V_potential, time)

Build Hamiltonian given Kinetic energy matrix and potential function Add potential to diagonal – specified by V_potential(x,t)

Parameters

• V_potential (function) – A potential function the caller must provide

• time (int) – Time dependence of the Hamiltonian. This argument is passed to the potential to simulate turning on a field pertubation, for example

Returns

Hamiltonian for the system defined by the caller provided V_potential, size nbas x nbas

Return type

H_mat (ndarray)

Hamiltonian_2D(H_mat_1D, V_coupling_mat, time)

1. Build 2D Hamiltonian from H_mat_1D matrix (nbas by nbas) and the coupling, V_coupling The two dimensions are equivalent with the same grid and same 1D (one-body) Hamiltonian

2. Add coupling potential (two-body potential) to diagonal – specified by V_coupling(r1,r1,t), function passed in as argument.

Parameters

• H_mat_1D (ndarray) – One-electron Hamiltonian (nbas by nbas)

• V_coupling_mat (function) – Coupling potential

• time (int) – Time dependence of the Hamiltonian. This argument is passed to the potential to simulate turning on a field pertubation, for example

Returns

Two-electron Hamiltonian

Return type

H_mat_2D (ndarray)

Hamiltonian_Two_States(V_potential_1, V_potential_2, V_coupling, time)

Build Hamiltonian given Kinetic energy matrix and potential function

Diagonal blocks are FEM-DVR Hamiltonians for V_potential_1 and Vpotential_2 Off diagonal blocks are diagonal representation ov V_coupling(x,t) Potential functions are passed in

Parameters

• V_potential_1 (function) – A potential function for the first state, the caller must provide

• V_potential_2 (function) – A potential function for the second state, the caller must provide

• V_coupling (function) – A potential coupling function providing the off-diagnol blocks of the returned potential, the caller must provide

• time (int) – Time dependence of the Hamiltonian. This argument is passed to the potential to simulate turning on a field pertubation, for example

Returns

Hamiltonian for the two state system defined by the caller provided V_potential_1, V_potential_2, and V_coupling, size nbas x nbas

Return type

H_mat (ndarray)

Kinetic_Energy_FEM_block(n, x, w)

Calculates block of Kinetic energy from one Finite-Element. If FEM-DVR grid has complex tail, then Kmat is complex

Parameters

• n (int) – the dimension of the derivative matrix

• x (ndarray) – array of ‘double’ or ‘complex’ DVR points

• w (ndarray) – array of ‘double’ or ‘complex’ DVR weights

Returns

square matrix of size n x n

Return type

Kmat (ndarray)

Plot_Psi(Psi_coefficients, plot_title_string='Plot of FEM-DVR representation', N_plot_points=500, make_plot=True)

Quick generic plot of function represented by FEM DVR coefficients

Parameters

• Psi_coefficients (ndarray) – Array of type ‘double’ or ‘complex’ coefficients for the representation of \(\Psi\) on the FEM-DVR grid

• plot_title_string (string) – Title of plot, defaults to “Plot of FEM-DVR representation”

• N_plot_points (int) – Number of points to plot, default 500

• make_plot (bool) – Boolean that turns off/on this plotting feature, default true

Returns

x and y coordinates of the graph

Return type

x_Plot, Psi_plot (ndarray, ndarrya)

Plot_Psi_2D(Psi_coefficients, plot_title_string='Plot of FEM-DVR representation', N_plot_points=50, make_plot=True)

Quick generic plot of function in represented by FEM DVR coefficients on 2D grid

Parameters

• Psi_coefficients (ndarray) – Array of type ‘double’ or ‘complex’ coefficients for the representation of \(\Psi\) on the FEM-DVR grid

• plot_title_string (string) – Title of plot, defaults to “Plot of FEM-DVR representation”

• N_plot_points (int) – Number of points to plot, default 500

• make_plot (bool) – Boolean that turns off/on this plotting feature, default true

Returns

x, y, and Psi values (Real parts)

Return type

x_Plot, y_Plot, Psi_plot (ndarray, ndarrya)

Plot_Psi_Two_States(Psi_coefficients1, Psi_coefficients2, plot_title_string='Plot of FEM-DVR representation', N_plot_points=500, make_plot=True)

Quick generic plot of wave function represented by FEM DVR coefficients which plots two states on the same grid that may be components of a coupled channels wave function.

Parameters

• Psi_coefficients1 (ndarray) – Array of type ‘double’ or ‘complex’ coefficients for the representation of first state \(\Psi_1\) on the FEM-DVR grid

• Psi_coefficients2 (ndarray) – Array of type ‘double’ or ‘complex’ coefficients for the representation of first state \(\Psi_2\) on the FEM-DVR grid

• plot_title_string (string) – Title of plot, defaults to “Plot of FEM-DVR representation”

• N_plot_points (int) – Number of points to plot, default 500

• make_plot (bool) – Boolean that turns off/on this plotting feature, default true

Returns

x and y coordinates of the graph

Return type

x_Plot, Psi_plot (ndarray, ndarrya)

Potential_Two_States(V_potential_1, V_potential_2, V_coupling, time)

Build potential function

Diagonal blocks are FEM-DVR Hamiltonians for V_potential_1 and Vpotential_2 Off diagonal blocks are diagonal representation of V_coupling(x,t) Potential functions are passed in

Parameters

• V_potential_1 (function) – A potential function for the first state, the caller must provide

• V_potential_2 (function) – A potential function for the second state, the caller must provide

• V_coupling (function) – A potential coupling function providing the off-diagnol blocks of the returned potential, the caller must provide

• time (int) – Time dependence of the Hamiltonian. This argument is passed to the potential to simulate turning on a field pertubation, for example

Returns

potential for the two state system defined by the caller provided V_potential_1, V_potential_2, and V_coupling, size nbas x nbas

Return type

potential (ndarray)

deriv(n, x, w)

Derivative of Unormalized lobatto shape functions. If FEM-DVR grid has complex tail, then derivative array is complex

Parameters

• n (int) – the dimension of the derivative matrix

• x (ndarray) – array of ‘double’ or ‘complex’ DVR points

• w (ndarray) – array of ‘double’ or ‘complex’ DVR weights

Returns

square matrix of size n x n

Return type

deriv_matrix (ndarray)

hello_world()

psi_evaluate(x, coef_vector)

Evaluate a function represented by a vector of coefficients of the FEM-DVR basis functions at the point x. Array containing vector of coefficients does not contain coefficients for the beginning and end of the FEM-DVR grid where boundary condx. enforce wavefunction = zero.

Find which finite element x is in (or on the boundary of)

Parameters

• x (complex) – FEM-DVR grid point of evalutation of \(\Psi\)

• coef_vector (ndarray) – Function representation in the FEM-DVR basis

Returns

The value of \(\Psi(x)\)

Return type

psi_value (complex)

psi_evaluate_2D(x, y, coef_vector)

def psi_evaluate_2D(self, x, y, coef_vector):

Evaluate a 2D function represented by a vector of coefficients of the FEM-DVR basis functions at the point x, y. Array containing vector of coefficients does not contain coefficients for the beginning and end of the FEM-DVR grid in either dimension where boundary conditions enforce wavefunction = zero.

Order of coefficients in coef_vector is i2d = i + (j-1)*nbas labels the (i,j) grid point x_i, y_j where nbas = number of DVR functions in one dimension. Grids in the x and y dimensions are the same.

Parameters

• x (complex) – FEM-DVR grid point of evalutation of \(\Psi\)

• y (complex) – FEM-DVR grid point of evalutation of \(\Psi\)

• coef_vector (ndarray) – Function representation in the FEM-DVR basis

Returns

The value of \(\Psi(x)\) at the point x, y

Return type

psi_value (complex)

rescale_kinetic(reduced_mass)

initialize KE_mat with the Kinetic energy matrix multiplied by \(\frac{1}{\mu}\), where \(\mu\) is the reduced mass

Parameters

reduced_mass (double) – The reduced mass used to rescale the KE matrix

Returns

rescaled KE matrix

Return type

KE_mat (ndarray)

quantumgrid.potential module

Class for a more OOP interface to several potentials.

class quantumgrid.potential.Potential(file=None)

Bases: object

Constructor method, added vectorized version of all the methods.

Parameters

file (string, optional) – The caller may provide a file to interpolate a potential onto the dvr grid. The file must be in two tab separated columns. Default to None

vectorized_V_morse

Vectorized version of the morse function

Type

ndarray

vectorized_V_Bernstein

Vectorized version of the Bernstein function

Type

ndarray

vectorized_V_c_state

Vectorized version of the c-state interpolated function of the cStateDCalc.csv file, which is NOT provided in the released package

Type

ndarray

vectorized_V_Interpolated

Vectorized version of the

Type

ndarray

interpolated function.

vectorized_V_Coulomb

Vectorized version of the Coulomb

Type

ndarray

function

Todo

Add a general interpolation scheme so any file passed into this class’s constructor will work. Right now, only works for a particular file that was tested in house and not released with package.

V_Bernstein(r: complex, time: int = 0.0) → complex

\(H_2\) potential from T-G. Wiechand R.B. Bernstein, J. Chem. Phys. 46 (1967) 4905. This is an accurate fit to the accurate Kolos and Wolneiwicz potential curve representation is valid \(0.4 <= R\) to infinity used in old ECS calculation in Julia Turner and C. William McCurdy, Chemical Physics 71(1982) 127-133 for resonances in dissociation for \(j .ne. 0\)

Note

ECS contour must begin beyond \(r = 9.5 a_0\) for safe analytic continuation

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• time (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value at the point r at the time t

Return type

pot (complex)

V_Coulomb(r: complex, time: int = 0.0) → complex

Coulomb potential for He or H- one-electron Hamiltonian

Note

Nuclear charge is set to 2.0

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• Znuc (double) – Charge on the residual ion

• t (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value on the Coulomb tail

Return type

pot (complex)

V_Interpolated(r: complex, time: int) → complex

Interpolated values using scipy CubicSpline

Note

Requires a file of potential values to interpolate in this class’s constructor!

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• t (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value at the point r at the time t

Return type

pot (complex)

V_c_state(r: complex, time: int = 0.0) → complex

Interpolate computed values using scipy CubicSpline \(\frac{1}{R^4}\) tail added matching value and finite diff derivative at \(R=5\)

Note

At this point constants are for Lucchese 4/3/2020 calculation: \(c ^4\Sigma_u^-\) state of \(O_2^+\) where the orbitals come from a SA-MCSCF on the ion using an aug-cc-vTZP basis set. This was for a cStateDCalc.csv file, which is NOT provided in the released package. Therefore, this potential is just a scaffold to implement a general interpolation scheme and shouldn’t be experimented with until then.

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• time (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value at the point r at the time t

Return type

pot (complex)

V_colinear_model(r1: complex, r2: complex) → complex

Colinear model for two-electron atom

Parameters

• r1 (complex) – FEM-DVR first point where this potential is evaluated at

• r2 (complex) – FEM-DVR second point where this potential is evaluated at

Returns

potential value at the point r1 and r2

Return type

pot (complex)

V_morse_1(r: complex, time: int = 0.0) → complex

Morse Potential defined by

\[V = d*(y^2 - 2*y)\]

with \(y\) defined by

\[y = e^{(-a*(r-re))}\]

This potential also defines parameters specifically for \(H_2\), De = 4.75 eV

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• t (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value at the point r at the time t

Return type

pot (complex)

V_morse_2(r: complex, time: int = 0.0) → complex

Morse Potential defined by

\[V = d*(y^2 - 2*y)\]

with \(y\) defined by

\[y = e^{(-a*(r-re))}\]

This potential also defines parameters specifically for \(H_2\)

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• t (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value at the point r at the time t

Return type

pot (complex)

V_morse_centrifugal(r: complex, time: int = 0.0) → complex

Morse Potential defined by

\[V = d*(y^2 - 2*y) + \mathrm{Centrifugal potential}\]

with \(y\) defined by

\[y = e^{(-a*(r-re))}\]

This potential also defines parameters specifically for \(H_2\)

Parameters

• r (complex) – FEM-DVR point where this potential is evaluated at

• t (int) – Time dependence of this potential to simulate turning on a field pertubation, for example. Defaults to t=0.0

Returns

potential value at the point r at the time t

Return type

pot (complex)

Module contents

Examples

There are four example scripts that come with quantumGrid: ecs_femdvr_time_indep_h2 for a time independent calculation and ecs_femdvr_time_dep_h2 for a time dependent calculation and two examples that calculate vibrational states of \(H_2\) and CO called femdvr_vib_states_h2 and femdvr_vib_states_co, respectively. Once quantumGrid is installed in your local environment, these scripts can be called without typing “python”. To see command line options for both scripts, just use the help command:

$ ecs_femdvr_time_dep_h2 --help

For example, to turn on plotting for the time independent example run the script with the plotting option:

$ ecs_femdvr_time_indep_h2 --want_to_plot=True

ecs_femdvr_time_indep_h2

For time-independent potential, this example implements Exterior Complex Scaling on the FEM-DVR contour. The value of R0 and the complex scale factor \(e^{I*theta}\) are specified. The representation of the potential must be able to be evaluated on the complex part of the contour.

Example:

Finds all eigenvalues of complex scaled Hamiltonian and plots any one of them, specified by n_Plot

Args:

1. want_to_plot (bool): Optional command that turns on plotting; default is false.

Potentials defined here:

1. Morse potential for \(H_2\)

2. Bernstein fit of Kolos and Wolneiwicz potential with \(\frac{1}{R^6}\), \(\frac{1}{R^8}\), \(\frac{1}{R^{10}}\) asymptotic behavior – Gives near spectroscopic accuracy used in [1], results there are reproduced by this code.

ecs_femdvr_time_dep_h2

Finds all eigenvalues of complex scaled \(H_2\) Hamiltonian for nuclear motion plots any one of them, specified by n_Plot Then starts a Gaussian packet in the well (e.g. with \(j=17\)) and follows it as it separates into a part that is bound in the well and a part that dissociates and vanishes on the ECS contour.

Args:

1. number_of_time_intervals (int): First command line argument. Number of time intervals to perform the Crank-Nicolson propagation; defaults to 300.

2. time_step (int): Second command line argument. Time step in the propagator to relax or restrict the calculation as needed; defaults to 0.1 atu.

3. want_to_plot_animate (bool): Optional command that turns on plotting; default is false.

Potentials defined here:

1. Morse potential for \(H_2\)

2. Bernstein fit of Kolos and Wolneiwicz potential with \(\frac{1}{R^6}\), \(\frac{1}{R^8}\), \(\frac{1}{R^{10}}\) asymptotic behavior – Gives near spectroscopic accuracy used in [1], results there are reproduced by this code.

femdvr_vib_states_h2

\(H_2\) vibrational states using CI singles and doubles potential curve from Psi4. This potential yields a \(n = 0 \rightarrow 1\) excitation energy within a few wavenumbers of the value using the NIST values for constants of diatomic molecules for \(H_2\) in the formula \(E_n = (n+\frac{1}{2})w_e - (n+\frac{1}{2})^2 w_ex_e\), which is \(4158 cm^{-1}\).

Shows how to

1. Read in and interpolate a potential function known at discrete points

2. Use FEMDVR class to build FEM-DVR grid

3. Use FEMDVR class to build Hamiltonian in DVR basis

4. Find eigenvalues and eigenvectors of Hamiltonian

5. Plot eigenfunctions of the Hamiltonian

Args:

1. want_to_plot (bool): Optional command that turns on plotting; default is false.

File read for this example:

1. potcurve_CISD_H2_ccpvTZ.dat

femdvr_vib_states_co

CO vibrational states using CI singles, doubles and triples potential curve from Psi4. This potential gives a dissociation energy of \(~12.2 eV\), not very good by comparison to the \(~11.1 eV\) experimental value. It yields a \(n = 0 \rightarrow 1\) excitation energy of \(2207 cm^{-1}\) compared with the value using the NIST values for constants of diatomic molecules for \(H_2\) in the formula \(E_n = (n+\frac{1}{2})w_e - (n+\frac{1}{2})^2 w_ex_e\), which is \(2143 cm^{-1}\) So not quite spectroscopic accuracy.

Shows how to

1. Read in and interpolate a potential function known at discrete points

2. Use FEMDVR class to build FEM-DVR grid

3. Use FEMDVR class to build Hamiltonian in DVR basis

4. Find eigenvalues and eigenvectors of Hamiltonian

5. Plot eigenfunctions of the Hamiltonian

Args:

1. want_to_plot (bool): Optional command that turns on plotting; default is false.

File read for this example:

1. potcurve_CISDT_CO_ccpvDZ.dat

time_dep_two_potential_excitation

In this example we demonstrate excitation from one potential curve (electronic state) to another in a diatomic molecule by a finite pulse using the FEM-DVR grid. The dipole matrix element between the two states is assumed to be constant as a function of internuclear distance. Potential curves are Morse oscillator functions with different well depths and shapes shifted by 0.15 hartrees. Reduced mass is reduced mass of H2. Note that a denser DVR grid may be necessary for heavier masses. Pulse has Sin^2 envelope with 3 femtosecond duration and is centered at 0.2 hartrees

Output includes an animation of the wave packet, plots of the potentials and initial wave packet, and text output of grid parameters, Hamiltonian eigenvalues, and properties of wave packets during the propagation at the plotting intervals in time.

Plot output, including an mp4 file of the animation is placed in the directory Plot_Output/

Uses the basic FEM-DVR functions in the FEM_DVR() class and also the functions specific to the case of nuclear motion on two potential surfaces:

1. Hamiltonian_Two_States: constructs the two Hamiltonian matrices and coupling

2. Crank_Nicolson_Two_States: propagates the two component wave function with coupling between them

3. Plot_Psi_Two_States: plots both components of the wave function (on the two potentials)

Args:

1. number_of_time_intervals (int): First command line argument. Number of time intervals to perform the Crank-Nicolson propagation; defaults to 300.

2. time_step (int): Second command line argument. Time step in the propagator to relax or restrict the calculation as needed; defaults to 0.1 atu.

3. want_to_plot_animate (bool): Optional command that turns on plotting; default is false.

two_electron

Temkin-Poet (s-wave limit) or colinear model of a two-electron atom: H- anion or He bound and autoionizing states

Finds all eigenvalues of complex scaled 2D Hamiltonian and plots any one of them, specified by n_state_plot. NB Diagonalization does not take advantage of sparse banded nature of the 2D Hamiltonian matrix in the FEM-DVR. Larger scale practical calculations must do so.

As an example the \(2s^2\) autoionizing state of He is chosen and the plot of the wave function shows the localized resonant state in the middle and autoionization decay down the sides parallel to the \(r_1\) and \(r_2\) axes. Numerical check: \(E_{gnd} = -2.8790288\) for He in s-wave limit singlet S autoionizing state \(E_{res} = -0.7228 - 0.001199 i\)

Note that the basis in this example is a simple product basis of DVR functions in \(r_1\) and \(r_2\), so both singlet (symmetric) and triplet (antisymmetric) spatial wave functions appear as eigenfunctions of the 2D Hamiltonian.

Modifying Scripts

The actual names of these four example scripts are ECS_FEMDVR_diatomic_time_indep_vibration_H2.py, ECS_FEMDVR_diatomic_time_dep_vibration_H2.py, H2_vib_states_FEM_DVR.py, CO_vib_states_FEM_DVR.py, Time-dep_excitation_2_potential_curves.py, and Two-electron_ECS.py. If you downloaded the source package from github, then these examples are in the examples directory. If quantumgrid was installed using the conda instruction then the scripts should be in /Path/to/Anaconda/envs/YOUR_ENVIRONMENT_NAME/lib/python3.7/site-packages/quantumgrid_examples. If you are in a Unix environment then you can simply find them with the following command:

$ locate ECS_FEMDVR_diatomic_time_dep_vibration_H2.py

At any rate, once found you can modify your script however you like!

References

1(1,2)

Julia Turner and C.William McCurdy. The application of exterior complex scaling in calculations on resonances in nuclear motion in molecular systems. Chemical Physics, 71(1):127 – 133, 1982. URL: http://www.sciencedirect.com/science/article/pii/0301010482870122, doi:https://doi.org/10.1016/0301-0104(82)87012-2.

Contributing

Contributions are welcome, and they are greatly appreciated! Every little bit helps, and credit will always be given.

You can contribute in many ways:

Types of Contributions

Report Bugs

Report bugs at https://github.com/zstreeter/quantumGrid/issues.

If you are reporting a bug, please include:

• Your operating system name and version.

• Any details about your local setup that might be helpful in troubleshooting.

• Detailed steps to reproduce the bug.

Fix Bugs

Look through the GitHub issues for bugs. Anything tagged with “bug” and “help wanted” is open to whoever wants to implement it.

Implement Features

Look through the GitHub issues for features. Anything tagged with “enhancement” and “help wanted” is open to whoever wants to implement it.

Write Documentation

quantumGrid could always use more documentation, whether as part of the official quantumGrid docs, in docstrings, or even on the web in blog posts, articles, and such.

Submit Feedback

The best way to send feedback is to file an issue at https://github.com/zstreeter/quantumGrid/issues.

If you are proposing a feature:

• Explain in detail how it would work.

• Keep the scope as narrow as possible, to make it easier to implement.

• Remember that this is a volunteer-driven project, and that contributions are welcome :)

Get Started!

Ready to contribute? Here’s how to set up quantumGrid for local development.

1. Fork the quantumGrid repo on GitHub.

2. Clone your fork locally:

   $ git clone git@github.com:your_name_here/quantumGrid.git
   

3. Install your local copy into a virtualenv. Assuming you have virtualenvwrapper installed, this is how you set up your fork for local development:

   $ mkvirtualenv quantumGrid
   $ cd quantumGrid/
   $ python setup.py develop
   

4. Create a branch for local development:

   $ git checkout -b name-of-your-bugfix-or-feature
   

   Now you can make your changes locally.

5. When you’re done making changes, check that your changes pass flake8 and the tests, including testing other Python versions with tox:

   $ flake8 quantumGrid tests
   $ python setup.py test or pytest
   $ tox
   

   To get flake8 and tox, just pip install them into your virtualenv.

6. Commit your changes and push your branch to GitHub:

   $ git add .
   $ git commit -m "Your detailed description of your changes."
   $ git push origin name-of-your-bugfix-or-feature
   

7. Submit a pull request through the GitHub website.

Pull Request Guidelines

Before you submit a pull request, check that it meets these guidelines:

1. The pull request should include tests.

2. If the pull request adds functionality, the docs should be updated. Put your new functionality into a function with a docstring, and add the feature to the list in README.rst.

3. The pull request should work for Python 3.7 and 3.8, and for PyPy. Check https://travis-ci.com/zstreeter/quantumGrid/pull_requests and make sure that the tests pass for all supported Python versions.

Tips

To run a subset of tests:

$ pytest tests.quantumGrid

Deploying

A reminder for the maintainers on how to deploy. Make sure all your changes are committed (including an entry in HISTORY.rst). Then run:

$ bump2version patch # possible: major / minor / patch
$ git push
$ git push --tags

Travis will then deploy to PyPI if tests pass.

Credits

Development Leads

• William (Bill) McCurdy https://chemistry.ucdavis.edu/people/william-mccurdy

• Giuseppe Barbalinardo http://giuseppe.barbalinardo.com/

• Zachary Streeter https://www.linkedin.com/in/zachary-streeter-44a323102/

Want To Help Develop?

Please contact Bill or Zachary if you would like to contribute!

History

0.0.1 (2020-05-31)

• First release on PyPI.

Indices and tables

• Index

• Module Index

• Search Page