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