NIC-CAGE: An Open-Source Software Package for Predicting Optimal Control Fields in Photo-Excited Chemical Systems
We present an open-source software package, NIC-CAGE (Novel Implementation of Constrained Calculations for Automated Generation of Excitations), for predicting quantum optimal control fields in photo-excited chemical systems. Our approach utilizes newly derived analytic gradients for maximizing the transition probability (based on a norm-conserving Crank-Nicolson propagation scheme) for driving a system from a known initial quantum state to another desired state. The NIC-CAGE code is written in the MATLAB and Python programming environments to aid in its readability and general accessibility to both users and practitioners. Throughout this work, we provide several examples and outputs on a variety of different potentials, propagation times, and user-defined parameters to demonstrate the robustness of the NIC-CAGE software package. As such, the use of this predictive tool by both experimentalists and theorists could lead to further advances in both understanding and controlling the dynamics of photo-excited systems.
Introduction<!>Theory<!>Numerical Methodology<!>Grid-Based Discretization of the Time-Independent Hamiltonian<!>Time Evolution<!>Numerical Optimization<!>Algorithm 1: NIC-CAGE<!>Numerical Examples and Results<!>Morse Potential<!>Asymmetric Double-Well Potential<!>Conclusions
<p>Nicolas Cage once said that "Every great story seems to begin with a snake." This work, unfortunately, is not one of those stories. Rather, the interaction between light and chemical/material systems is central to a multitude of technological applications, including energy-efficient photovoltaic materials, [1,2,3] plasmon-induced energy transfer, [4,5,6,7] and sustainable photocatalysis. [8,9,10] Our capability to fully harness these light-induced systems has tremendous potential to grow further as we build our understanding on how to control the excited-state quantum dynamical processes in these systems. While the majority of theoretical and experimental studies have focused on how these systems interact with a given electromagnetic radiation source (i.e, a laser pulse or an external light source), fewer studies have investigated the inverse problem to ask "If we desire a specific behavior in a chemical/material system, can we construct a light source to achieve this, and what does its functional form look like?". [11,12,13] To help address this "inverse" question, this work provides an open-source software package, NIC-CAGE (Novel Implementation of Constrained Calculations for Automated Generation of Excitations), for predicting optimized electric fields to control photo-excited chemical systems. With a pre-computed electronic potential energy, V(x) (which can be routinely obtained with most electronic structure software packages), the NIC-CAGE code calculates optimal electric fields, (t), that can drive a system from a known initial vibrational eigenstate to another desired state, as depicted in Fig. 1. In particular, we provide analytic gradients for maximizing the transition probability based on a norm-conserving Crank-Nicolson propagation scheme, which have not been previously reported in the literature. These analytic gradients allow us to efficiently construct a control field, (t), that gives a large (≈ 1) transition probability for (constrained) excitations from an initial to a specified final quantum state. The NIC-CAGE code is provided in the MATLAB and Python programming environments as two separate but fully self-contained implementations. We have chosen to implement this code in these high-level programming environments to allow general researchers to easily understand how these approaches and tools are used in practice to control photo-excited mechanisms in a variety of systems.</p><p>It is worth noting that there are other software implementations of optimal control for quantum systems, such as the QuTiP package. [14] However, our formulation differs from the former in two important ways: (1) the QuTiP package is designed specifically for the optimal control of spinor systems and does not currently treat real-space potentials (i.e., potentials as a function of a real-space coordinate x, which are the focus of the NIC-CAGE code), and (2) the NIC-CAGE code is formulated with a Crank-Nicolson time-propagation scheme and its analytic gradients, whereas the GRAPE algorithm [15] in the QuTiP package uses a full matrix exponential to approximate the full time evolution operator. In a very recent study, Gharibnejad and co-workers carried out extensive timings of various numerical solutions to the one-dimensional time-dependent Schrödinger equation and found that higher-order time-propagation methods (such as split-operator or matrix exponential methods) allow larger time-steps to be taken; however, the number of floating-point operations in these methods increase significantly such that this trade-off does not lead to a performance advantage. In particular, the study by Gharibnejad et al. concluded that the CrankNicolson approach was more efficient than several other higher-order time-propagation methods when small time steps were used. Since the numerical propagation scheme is still the most time consuming part of any quantum control algorithm, the Crank-Nicolson-based approaches and analytic gradients presented in this work could be useful for the control of larger quantum systems in which the split-operator or matrix exponential propagation schemes become too costly. [16] Intermolecular Distance (x)</p><p>Figure 1: Schematic of an arbitrary potential energy well, V(x), with various allowed energy levels (E 0 , E 1 , and E 2 ) and their associated probability densities, |ψ(x)| 2 (drawn in blue). The optimized control fields, (t), enable transitions between the desired vibrational eigenstates.</p><p>The present paper is organized as follows: Section 2 briefly outlines the basic theoretical concepts needed to understand the general idea for controlling the time evolution of a photo-excited system. Section 3 gives a more detailed derivation of our Crank-Nicolson analytical gradients and describes the numerical implementation for carrying out this time-evolution operation. With the basic concepts and mathematical formulation in place, Section 4 provides several examples and outputs on a variety of different potentials, propagation times, and user-defined parameters to demonstrate the robustness of the NIC-CAGE software package. Finally, Section 5 concludes with a brief discussion and perspective look at potential future applications of our approach.</p><!><p>Within this section, we outline the basic theoretical concepts for optimizing/controlling the time evolution of a molecular system under the influence of electromagnetic radiation (i.e, a laser pulse or an external light source). A more detailed description of our numerical implementation and approach for carrying out this time-evolution optimization is given in Section 3. Although there have been previous studies on optimizing quantum control fields (each with their own purposes and advantages) [11,12,13,17], we focus on the grid-based Crank-Nicolson discretization utilized in this work since (1) the analytic gradients for this norm-conserving approach have not been previously reported in the scientific literature, and (2) the resulting numerical scheme is more straightforward to implement and understand compared to other approaches.</p><p>The temporal dynamics for describing the motion of nuclei in a molecular system are governed by the timedependent Schrödinger equation, which, in atomic units is given by:</p><p>where the time-dependent Hamiltonian operator, H(x, t), is</p><p>In the expression above, x denotes the reduced coordinate along a reaction path (cf. Fig. 1), m is the effective mass associated with the molecular motion along the reaction path, and V(x, t) is a time-dependent potential energy function (additional details on both x and V(x, t) are given later). The effective mass, m, is closely related to Wilson's G-matrix formalism [18], and previous work by us [19] gives details and an open-source code for calculating this molecular parameter along the reaction path. The complex-valued molecular wavefunction, ψ(x, t), in Eq. ( 1) represents the probability amplitude for the motion of the nuclei along the reduced coordinate path, and |ψ(x, t)| 2 can be interpreted as a probability density with the normalization condition</p><p>which holds for all t. In this work, we utilize the dipole approximation [5,4,20] in which the time-dependent potential can be written as</p><p>where V(x) is the time-independent portion of the potential that contains the Born-Oppenheimer electronic energy of the molecule along the reaction path (which can be obtained from a standard quantum chemistry calculation). The dipole moment function, µ(x), can also be readily computed from an electronic structure calculation at every point, x, along the reaction path as described in Ref. [21]. Finally, (t) is the time-dependent external electric field whose functional form is iteratively optimized using the NIC-CAGE software package in this work. At this point, it is worth mentioning a few choices for the reduced coordinate, x, given in Eqs. ( 2) and (4), respectively. One of the most common choices for x is the set of mass-weighted coordinates along the minimum energy path (also known as the intrinsic reaction coordinate), which is defined as the steepest descent path from a transition state toward reactants or products [22,23]. However, obtaining this quantity in even one dimension can be computationally expensive since it typically requires a full computation of the molecular Hessian at each point along the reaction path. It is also worth mentioning that while algorithms for computing the intrinsic reaction coordinate in one dimension are available in many quantum chemistry packages [22,23], no such algorithm is available for calculating minimum energy paths in two or higher dimensions. This practical difficulty arises since the couplings between the molecular motions cannot be removed to satisfy the Eckart conditions in two or higher dimensions [24] (the Eckart condition is automatically satisfied along the one-dimensional intrinsic reaction coordinate [25,26]). Second, there have been several contemporary experimental studies [27,28] and a recent review [29] showing that quantum control calculations on one-dimensional slices of the multi-dimensional potential energy surface still capture most of the dynamical effects of these systems. For these reasons, we only focus on quantum control calculations along one-dimensional potentials since (1) two-or higher-dimensional intrinsic reaction paths (or their associated dipole moment surfaces, which are required for quantum control) are not readily computable and (2) much of the relevant dynamics can still be gleaned from the one-dimensional slices of the potential energy surface. To this end, prior work by us [21,30,31,32,33,34] and other researchers [26,35] have shown that an approximate path (which differs from the intrinsic/minimum path) can be parameterized with a single internal coordinate such as a bond length [21], valence bend angle [32,33,34], or dihedral angle [30,31] that can accurately describe reactions involving bond dissociation, isomerization, or internal rotation, respectively. Once a suitable reduced coordinate, x, is chosen, both V(x) and µ(x) can be readily computed in most quantum chemistry packages such as Gaussian [36], Q-Chem [37], GAMESS [38], or NWChem [39] by carrying out a relaxed potential energy scan.</p><p>With x and V(x) properly chosen/computed, Eqs. ( 1) and (2) allow us to mathematically answer the question: "Given an electric field (t), how does an initial state, ψ 0 (x, t = 0), evolve after some final time T has elapsed?" However, as mentioned in the Introduction, we instead seek the answer to the following "inverse" question: "If we want to reach a desired final state ψ N−1 (x, t = T ) at time T (after N − 1 propagation steps), what does the functional form of (t) look like?" Providing accurate and efficient answers to this inverse question is the ultimate goal of the NIC-CAGE software package described in this work. To be more mathematically precise, we seek the functional form of an external electric field, (t), that maximizes the transition probability, P ψ N−1 (x) , given by</p><p>where ψ f is a known desired final target wavefunction (given by the user and discussed further below), and ψ N−1 is obtained after applying N − 1 successive propagation steps of the time-dependent Schrödinger equation (i.e., Eq. 1). It should be noted that the transition probability is essentially a measure of the similarity of the final target and the propagated wavefunction. As written, the maximization of Eq. ( 5) is an unconstrained optimization problem to numerically solve for the time-dependent, external electric field (t). However, to make a direct connection to realistic experiments that may have limited power constraints for generating the optimal (t), we can augment Eq. ( 5) with a fluence penalty term:</p><p>where the fluence, F [ ], acts as a penalty to prevent unphysically large values of the electric field and is given by:</p><p>where α is a positive constant to be chosen by the user. For the remainder of this paper, we focus on the constrained maximization of J ψ N−1 (x), given by Eq. ( 6).</p><p>With the quantum optimal control problem now properly defined, one can next choose the desired initial and final molecular wavefunctions to be used as initial conditions in the numerical quantum control procedure. In principle, both the initial and final wavefunctions can take on any form and are not even required to be stationary states of the Schrödinger equation (i.e., the numerical optimization techniques used in the NIC-CAGE software package can be used equally well for stationary or non-stationary states). In practice, however, one typically desires a control field that drives the system from a known initial vibrational eigenstate, ψ i=ν i (x), to another desired vibrational eigenstate, ψ f =ν f (x), of the system (ν i and ν f denote the vibrational quantum numbers of the initial and final wavefunction, respectively). As such, to obtain these stationary wavefunctions as input, the NIC-CAGE software package contains a numerical routine for solving the time-independent Schrödinger equation; i.e., [− 1 2m dψ(x)</p><p>Put more concretely, we choose our initial state ψ 0 (x, t = 0) to be a user-specified stationary vibrational eigenstate ψ i=ν i (x) of the time-independent Schrödinger equation, and the desired final state ψ N−1 (x, t = T ) is chosen to be a different vibrational eigenstate ψ f =ν f (x). To make our notation more succinct in the remainder of this paper, ψ i=ν i (x) and ψ f =ν f (x) are abbreviated as ψ i (x) and ψ f (x), respectively.</p><!><p>With the basic theoretical concepts discussed in the previous section, we now give a detailed description of the various numerical approaches used in the NIC-CAGE code for calculating optimized control fields that drive the system from ψ i (x) to ψ f (x). Section 3.1 gives the numerical scheme used to discretize the time-independent portion of the Hamiltonian across a spatial grid. With the time-independent Hamiltonian properly discretized, Section 3.2 describes the Crank-Nicolson approach used in the NIC-CAGE code to numerically propagate the full time-dependent Schrödinger equation. Finally, Section 3.3 describes the gradient ascent approach for maximizing an objective functional of the transition probability P. Throughout the following subsections, we use the following notation: scalars are denoted by italic letters (x), vectors by bold lower case letters (x) (unless specified otherwise), matrices by bold upper case letters (X), the element in the k th row and l th column of matrix X by [X] kl , the identity matrix by I, where I (m±) denotes a matrix with entries of 1 on the mth diagonal either above (+) or below (-) the main diagonal, a diagonal matrix having the entries x along the diagonal by diag(x), the Schur (element-wise) product by the symbol, and the conjugate and conjugate transpose of X by X * and X † , respectively.</p><!><p>We commence by discretizing the kinetic energy portion of the Hamiltonian operator across a grid of L equidistant points on the interval [x min , x max ]. Using a fourth-order five-point stencil finite difference approximation gives</p><p>where ∆x = x max −x min L−1 . The notation of Eq. ( 8) can be considerably simplified by defining the ith grid point as x i = x min + i∆x (where i = 0, . . . , L − 1) and ψ(x i ) = ψ i . With this simplified notation, the discretized time-independent Schrödinger equation can be written as</p><p>The notation in Eq. ( 9) can be mathematically interpreted as a matrix-vector equation, and the time-independent Hamiltonian can be written as the following matrix</p><p>where I is an L × L identity matrix with entries of 1 on the main diagonal, I (1±) are L × L matrices with entries of 1 on the 1st diagonal either above (1+) or below (1-) the main diagonal, I (2±) are L × L matrices with entries of 1 on the 2nd diagonal either above (2+) or below (2-) the main diagonal, and V is an L × L diagonal matrix with entries [V] i j = V(x i )δ i j , where δ i j is the Kronecker delta. As mentioned in Section 2, the time-independent Hamiltonian matrix in Eq. ( 10) is diagonalized to obtain the initial and final states,</p><p>respectively, as input to the gradient ascent optimization in the NIC-CAGE software package.</p><!><p>With the time-independent Schrödinger equation properly discretized in Section 3.1, we can now proceed with the general case for a time-dependent Hamiltonian and its discretized form. By combining explicit and implicit Euler schemes as in Ref. [40], one obtains the following Crank-Nicolson scheme for a time-dependent Hamiltonian:</p><p>where τ = T N−1 is the time step across a grid of N equidistant points on the interval [0, T ]. It is important to note that the Crank-Nicolson scheme in Eq. ( 11) is unitary and, therefore, preserves the norm of the wavefunction. We are now in a position to discretize time in a similar manner as we have previously discretized space in Section 3.1. Using a similar notation scheme as before gives ψ(x i , t j ) = ψ(x min + i∆x, jτ) = ψ i j where i = 0, . . . , L − 1 and j = 0, . . . , N − 1. With this simplified notation, in conjunction with the matrix notation of Eq. ( 10), the Crank-Nicolson expression in Eq. ( 11) can be written as the following matrix equation:</p><p>where ψ j is vectorized in space (x) and evaluated at time t j = jτ, where j = 0, . . . , N − 1 (i.e., ψ j is a column vector with entries ψ i j ). In addition, V j+1/2 is a diagonal matrix (cf. Eq. ( 10)) evaluated at time t j+1/2 = ( j + 1/2)τ, where j = 0, . . . , N − 2. It should also be noted that Eq. ( 12), by definition, is only valid for values of j = 0, . . . , N − 2, and ψ 0 , therefore, represents the initial state.</p><p>As will be shown in Section 3.3, it is helpful to write Eq. ( 12) in a more compact notation. To this end, we denote the operator in parentheses on the left-hand side of Eq. ( 12) as the L × L matrix U j+1/2 . Similarly, we represent the operator on the right-hand side of Eq. ( 12) by the L × L matrix W j+1/2 , which results in the following matrix equation:</p><p>where U j+1/2 and W j+1/2 are banded pentadiagonal L × L matrices of the form</p><p>Upon closer examination of Eq. ( 12), one finds that the upper and lower diagonals of Eq. ( 14) are all constants, while the diagonal entries are functions of the time-dependent potential</p><p>), as summarized below:</p><p>Since the matrix-vector product on the right-hand side of Eq. ( 13) is computed several times for each time step, significant computational savings are utilized in the NIC-CAGE software package by replacing the matrix-vector product with a single pre-computed vector with the following elements:</p><p>where</p><!><p>The numerical procedure presented in the previous section gives the time-evolution of the wavefunction for a given j+1/2 where j = 0, . . . , N − 2. In this section, we derive analytic gradients and present a numerical maximization method (using a gradient ascent algorithm) that back-propagates the gradient of the objective functional, J ψ N−1 (x), :</p><p>where we have used the definition of J ψ N−1 (x), from Eq. ( 6). To make our discussion of this numerical optimization more concrete (and to remind the reader of the specific time steps used in this calculation), Fig. 2 shows a schematic of how the electric field optimization process can be conceptualized. As mentioned in Section 3.2, the time interval [0, T ] is discretized across a grid of N equidistant points (endpoints inclusive); however, the electric field amplitudes j+1/2 are only defined at times t j+1/2 = ( j + 1/2)τ, where j = 0, . . . , N − 2. Within the small duration τ = T N−1 , the electric field amplitudes, j+1/2 , are approximated to be constant. The vertical arrows in Fig. 2 represent the gradients dJ d j+1/2 that indicate how each amplitude j+1/2 should be modified in the next iteration to maximize the objective functional J. It is worth mentioning that our approach for representing/optimizing the electric field is similar in spirit to the GRAPE algorithm used to design pulse sequences in NMR spectroscopy [15]. In the original GRAPE algorithm, a full matrix exponential is computed for each τ to approximate the full time evolution operator (which is typically carried out with a fully dense algorithm since the exponential of a matrix is nearly always dense [41]). However, as mentioned in the Introduction, our formulation differs algorithmically from the former since the NIC-CAGE code is formulated with a Crank-Nicolson time-propagation scheme and its analytic gradients, which have been shown to be significantly faster than other split-operator or matrix exponential methods for small time intervals [16]. To enable this numerical optimization across a grid of points, we first approximate the transition probability in Eq. ( 6) as</p><p>where ψ f is a known final target wavefunction (given by the user, as described in Section 2), and ψ N−1 is obtained after applying N − 1 iterations of Eq. ( 13). For clarification, Eq. ( 19) uses the same notation as Section 3.2, where ψ N−1 is a column vector with entries ψ i N−1 and ψ f † is the complex conjugate transpose of ψ f (i.e., ψ f † is a row vector with entries (ψ i f ) * ). With these definitions, it follows that P ≤ 1. Similarly, the fluence in Eq. ( 7) can be approximated as</p><p>It is worth noting that the expression for the fluence in Eq. ( 20) could be generalized to other constrained solutions by introducing a time dependence for the penalty factor, α (i.e., to produce electric field shapes that are experimentally realizable). For example, one could choose α(t) = α 0 s(t) , where α 0 is a constant and s(t) serves as a shape function (such as a Gaussian profile) for the laser pulse. The effect of such a shape function would force the value of α(t) to go to infinity when s(t) approaches zero. That is, the penalty factor in those temporal regions would be very high, resulting in field intensities to be very low, which would allow for a smooth on/off switching behavior of the designed field. Within the NIC-CAGE code, we only focus on the situation where α is a positive constant (which also prevents a novice user from inputting functional forms of s(t) that would have zeroes and cause numerical overflows in the code), although the implementation of this general constraint could serve as interesting modifications to explore in future work.</p><p>Returning to the problem at hand, to evaluate Eq. ( 18), we commence by first obtaining expressions for the derivatives of P with respect to each j+1/2 value using the chain rule of Wirtinger's derivatives [42].</p><p>where</p><p>As such, our gradient ascent algorithm culminates in finding an expression for dψ N−1 d j+1/2 . To obtain a numerically tractable expression for this derivative, we note that ψ j+1 can be propagated from ψ j by inverting Eq. ( 13); i.e.,</p><p>where we recall that the above expression, by definition, is only valid for values of j = 0, . . . , N − 2 (cf. Eq. ( 12)). It follows that repeated iteration of Eq. ( 23) gives</p><p>We have specifically chosen to write Eq. ( 23) in this suggestive form since only the diagonal entries of U j+1/2 and W j+1/2 are a function of j+1/2 (cf., Eqs. ( 13), (15), and ( 16)). Therefore:</p><p>where diag iτ 2 µ denotes a diagonal matrix having the entries iτ 2 µ along the diagonal, and N−2 k= j+1 U −1 k+1/2 W k+1/2 denotes a "time-ordered product," where matrices at earlier times are to the right of matrices at later times. Eq. ( 25) can be simplified even further by using the following identity for the derivative of a matrix inverse:</p><p>where the symbol denotes the Schur (element-wise) product (i.e., µ ψ j+1 + ψ j is a column vector with entries</p><p>. It should also be mentioned that dψ N−1 d j+1/2 , by definition, is only defined for values of j = 0, . . . , N − 2 (see Eqs. ( 12), (23), and recall that the "latest" electric field that ψ N−1 can depend on is the previous N−3/2 , not the future N−1/2 ). However, it should be pointed out that Eq. ( 26) cannot be directly evaluated for j = N − 2 since the lower index of the product symbol then becomes larger than the upper index. To evaluate dψ N−1 d j+1/2 for the special case when j = N − 2, we can use the original definition of ψ N−1 to obtain</p><p>Using the same mathematical techniques and identities as in Eqs. ( 25) -( 26) finally gives</p><p>For convenience, we summarize the closed form expressions for dψ N−1 d j+1/2 for all possible values of j:</p><p>Now that we have a closed form expression for dψ N−1 d j+1/2 , the only remaining quantity needed to evaluate the gradient of the objective functional (i.e., dJ[ψ N−1 (x), ]</p><p>in Eq. ( 18)) is the derivative of the fluence term (Eq. ( 20)), which fortunately has a relatively simple form:</p><p>With these results, we can now numerically solve for the optimized, time-dependent, external electric field at each l th iteration step, (l) j+1/2 , to give (l)</p><p>where we have made use of Eqs. ( 18), ( 21), (22), and (30) in the above expression. Finally, dψ N−1 d j+1/2 can be calculated using Eq. ( 29), and γ in Eq. ( 31) is the learning rate of the gradient ascent algorithm, which we calculate using the bisection line-search approach [43] described in Ref. [44]. In short, the gradient ascent algorithm proceeds until the transition probability defined in Eq. ( 19) reaches a specified tolerance or exceeds a maximum number of iterations (set to 0.99 and 100, respectively, in the NIC-CAGE code). It is also important to note that this algorithm (as with most other gradient ascent approaches) can only converge to a local maximum, which may not necessarily be the global maximum. We conclude this section with the following algorithm flowchart that summarizes the various inputs and algorithmic steps used in the NIC-CAGE software package to calculate the optimized electric field (t):</p><!><p>Input: spatial interval [x min , x max ], grid spacing ∆x, time interval [0, T ], time step τ, mass m, dipole moment function µ(x), potential energy function V(x), initial state number i, desired final state number f Output: initial wavefunction ψ i (x), desired final wavefunction ψ f (x), final propagated wavefunction ψ N−1 (x), optimized electric field (t), power spectrum of optimized electric field 1. Diagonalize time-independent Schrödinger Equation in Eq. ( 10) to obtain ψ i (x) and ψ f (x) 2. Initialize j+1/2 = 0 for j = 0, . . . , N − 2 3. Initialize P = 0 and iter = 0 while P < 0.99 & iter < 100 do for j = 1 to N − 1 do Calculate ψ j from Eq. ( 12) end Update J and P using Eq. ( 6) for j = N − 2 to 0 do Calculate dJ d j+1/2 from Eq. ( 18) end Calculate γ with bisection linesearch method Update j+1/2 for j = 0, . . . , N − 2 using Eq. ( 31) Set iter = iter + 1 end return j+1/2</p><!><p>As mentioned in Section 2, the time-independent portions of Eq. ( 4) [i.e., V(x) and µ(x)] can be obtained from a quantum chemistry calculation or take any arbitrary form. To demonstrate the flexibility and capabilities of the NIC-CAGE software package, we give a few representative results of running the NIC-CAGE code on two example potentials that mimic a bond stretching/dissociation process in a photocatalytic reaction. All of the parameters for both potentials (expressed in atomic units) are included in the easy-to-use demo1.m, demo2.m, and demo3.m scripts.</p><!><p>For our first example, we perform a quantum control optimization using a Morse potential that mimics the photoinduced stretching of an O-H bond [45]. The spatial dependence, V(x), of the Morse potential that we used in Eq. ( 4) has the following form:</p><p>where the well depth, width, and equilibrium bond distance were taken from Ref. [45] and are given by D e = 0.1994, β = 1.189, and x 0 = 1.821, respectively. In addition, the dipole moment function that we used in Eq. ( 4) has the functional form</p><p>where µ 0 = 3.088 and x * = 0.6. The value of m that we used in Eq. ( 2) is the reduced mass of oxygen and hydrogen.</p><p>A snippet from the included demo1.m MATLAB script given below provides values of the parameters used in this section: These values correspond to the variables described previously in this paper as follows:</p><p>• x_min, x_max, dx: spatial range specified by [x min , x max ], grid spacing ∆x • T, tau: propagation time T , time step τ • m: effective mass m associated with the molecular motion along the reaction path • mu: dipole moment function µ(x), cf. Eq. ( 33)</p><p>• alpha: fluence penalty pre-factor α, cf. Eq. ( 7)</p><p>The quantum control optimization commences with an initial guess of (t) = 0 with α = 1. It is worth mentioning that this specific choice of initial conditions works for all the examples presented in this paper; however, larger energy transitions may require a non-zero initial guess. Throughout this section, we use the notation mentioned at the end of Section 2 where transitions between the initial and final vibrational eigenstates, ψ i (x) and ψ f (x), are simply denoted by ν i and ν f , respectively, where ν = 0 represents the zero-point ground state. Fig. 3 depicts the values of both the objective functional, J, and transition probability, P, as a function of the number of iterations (iteration zero corresponds to the initial values of J or P with an initial guess of (t) = 0). For the transition ν i = 0 → ν f = 1, only 3 steps are required to reach a transition probability of 99.2%. Since the vertical axis in Fig. 3 is scaled linearly, it would be incorrect to infer that the improvement at iteration number one is minimal. We can fully appreciate the progress made at the first iterative step by computing the percentage change in J from the previous step. For instance, at the first iteration, the value of J increases from 1.03 × 10 −27 to 6.19 × 10 −9 , which translates to an increase by 6.01 × 10 20 % -the highest among all other iterations. This overall trend is consistent with all the numerical examples discussed in this paper. Fig. 4 summarizes the various outputs and post-processed results of the NIC-CAGE software package for the ν i = 0 → ν f = 1 transition in a Morse potential. The panel in Fig. 4a shows the norm-squared vibrational eigenstate, |ψ i (x)| 2 , with ν i = 0 (dashed red line), superimposed on the Morse potential energy function (solid blue line). The optimized electric field in Fig. 4b is nearly a pure sinusoid, which is also reflected in the frequency domain as a sharp peak at around 112 THz in Fig. 4c. It is worth noting that the peak frequency matches the resonant frequency required for the transition between the ground and first excited state. As a further check on the fidelity of the results, Fig. 4d shows that the norm-square of the propagated wavefunction, |ψ N−1 | 2 , closely matches the desired target |ψ f | 2 .</p><p>Since the optimized electric field in this example is nearly mono-chromatic, it is natural to ask whether a simple π pulse could achieve the same desired target state. To test this possibility, we carried out several tests in the Supplementary Information that compare propagated wavefunctions obtained from analytical π pulses vs wavefunctions resulting from numerically optimized electric fields predicted by the NIC-CAGE program. In particular, the first test in the Supplementary Information shows that the propagated wavefunction obtained with the π pulse attains a transition probability of 0.991, which indicates that the electric field is not unique and can differ within a few tenths of a percent with vastly different mechanisms. Subsequent tests with the π pulse, however, also show this approach is not completely general and can even give sub-optimal results for simple modifications to the potential (cf. Fig. 2a in the Supplementary Information). In contrast, the numerically optimized electric fields obtained from the NIC-CAGE program provide a more general approach for obtaining these desired transitions, particularly in situations where the electric field may be nontrivial (as demonstrated by the other examples in this work). Finally, it should be mentioned that the propagated wavefunction, ψ N−1 , is complex-valued in general and has both real and imaginary parts (not shown in Fig. 4d); nevertheless, the transition probability</p><p>is real-valued and close to unity, emphasizing the ability of the NIC-CAGE software package to directly optimize the control field (t) in this system. In our next example, we attempt a more difficult case with the NIC-CAGE software package that involves a transition between non-adjacent energy levels (included in the demo2.m MATLAB script). Fig. 5 shows the convergence behavior of both the objective functional, J, and transition probability, P, as a function of the number of iterations for the ν i = 3 → ν f = 5 transition. After only two iteration steps, over 88% of the maximum objective functional value is achieved, and the subsequent steps give an incremental improvement before terminating at the 44th step with 99.1% probability. It is interesting to note that terminating the optimization prematurely after only two iterations results in an electric field that is almost a pure sinusoid at approximately 188 THz. As such, the subsequent iterations clearly show that the fluence term (i.e., Eq. ( 20)) significantly impacts the overall objective functional value with the algorithm effectively penalizing the dominant frequency component, allowing for multi-photon transitions to occur (discussed further below). Fig. 6 summarizes the various outputs and post-processed results of the NIC-CAGE software package for the ν i = 3 → ν f = 5 transition in a Morse potential. The panel in Fig. 6a shows the norm-squared initial vibrational eigenstate, |ψ i (x)| 2 , with ν i = 3 (dashed red line), superimposed on the Morse potential energy function (solid blue line). The optimized electric field in Fig. 6b is signficantly more complicated than the ν i = 0 → ν f = 1 case discussed previously. The power spectrum of the optimized electric at the end of 44th iteration step, shown in Fig. 6c, now contains three peaks that correspond to the ν = 3 → 4, ν = 4 → 5, and ν = 3 → 5 transitions. While the single-photon transition has a higher intensity than the peaks corresponding to the ν = 3 → 4 → 5 transition, this example shows that using a single "trivial" frequency corresponding to the energy difference between the ground and excited state may not be sufficient, and a fully numerical optimization as implemented in the NIC-CAGE code is necessary. Finally, as a check on the fidelity of the results, Fig. 6d</p><!><p>To test the robustness of the NIC-CAGE software package, we now perform a quantum control optimization on a different potential energy function, namely an asymmetric double-well potential:</p><p>The functional form of this double-well potential captures all the basic physics of a donor-acceptor reactant/product photochemical system. Specifically, the initial ground state with ν i = 0 is localized in the deeper "electron donor reactant" left well, and the desired vibrationally-excited final state with ν f = 1 is primarily confined to the "electron acceptor product" well on the right (cf. Figs. 8a and 8d). In addition to examining a completely new potential energy function, this particular example investigates the effects of a significantly shorter pulse on the optimization. For this specific example, we chose a unit mass with a dipole moment given by µ(x) = x. The following code snippet from the included demo3.m MATLAB script shows only those initialization variables that differ from the previous section: For the ν i = 0 → ν f = 1 transition in the asymmetric double-well potential, we examined the effect of using two different propagation times: T = 1000 and 100 a.u. As illustrated in Fig. 7a for T = 1000, we observed a rather quick convergence to a transition probability of 1.000 within 5 iteration steps. The corresponding optimized electric field for this case, not shown for brevity, is essentially a pure sinusoid at 1,041 THz, the resonant frequency associated with the desired transition. In contrast, for T = 100, the transition probability reaches 0.975 only after a lengthy 100 iteration steps (cf. Fig. 7b). Since the solutions of the asymmetric double-well potential for T = 1000 are relatively trivial, we only present the various outputs and post-processed results of the NIC-CAGE software package for the T = 100 case in Fig. 8. The panel in Fig. 8a shows the norm-squared initial vibrational eigenstate, |ψ i (x)| 2 , with ν i = 0 (dashed red line), superimposed on the asymmetric double-well potential energy function (solid blue line). The optimized electric field in Fig. 8b is clearly a superposition of numerous frequencies, which is more complicated than either of the cases discussed previously. This complexity is reflected in the power spectrum of the optimized electric field shown in Fig. 8c, which contains a multitude of peaks. Finally, as a check on the fidelity of the results, Fig. 8d shows that the norm-square of the propagated wavefunction, |ψ N−1 | 2 , closely matches the desired target |ψ f | 2 . Taken together, both this example and the Morse potential discussed previously demonstrate the robustness of the NIC-CAGE software package for handling different potentials, propagation times, and various user-defined parameters.</p><!><p>In this work, we have fully described, documented, and provided the NIC-CAGE software package -an opensource code for predicting optimized electric fields, (t), to control photo-excited chemical systems . Specifically, the NIC-CAGE code utilizes analytic gradients with a grid-based (norm-conserving) Crank-Nicolson propagation scheme to optically drive a chemical system from a known initial vibrational eigenstate to another desired state. Our approach is implemented in the high-level MATLAB and Python programming languages to allow researchers to get a detailed "look under the hood" to understand how our approaches are numerically incorporated in practice. Most importantly, we have tested the NIC-CAGE software package on a variety of different potentials, propagation times, and user-defined parameters to demonstrate its robustness in obtaining transition probabilities with extreme (typically over 97%) accuracy.</p><p>Looking forward, we anticipate that the NIC-CAGE software package could be used in a multitude of applications to control photo-excited chemical systems. For instance, as demonstrated in our asymmetric double-well example, the optimal electric field, (t), can be strongly dependent on the total propagation time T . As such, the flexibility in choosing T provides an opportunity to use the NIC-CAGE program to examine how this parameter can control the underlying dynamics. For example, if a low-intensity electric field is desired, then a long propagation time T should be chosen. Alternatively, one may desire a short propagation time to overcome rapid relaxation mechanisms in a system (such as intramolecular energy transfer) [46] or other experimental physical limitations. In the same vein of exploring parameter space, one can also use the NIC-CAGE program to explore the effect of altering the parameter, α, in the fluence penalty term (cf. Eq. ( 7)). For instance, if a low-power electric field is desired, a larger value α should be chosen to penalize high pulse intensities and meet the experimental limitations of the user. Similarly, as discussed in Section 3.3 (cf. Eq. ( 20)), one can also modify α to have the form of a shape function such as a Gaussian pulse shape, which may be more experimentally accessible.</p><p>Finally, the NIC-CAGE software package could also be used as a predictive tool by experimentalists to give rigorous bounds on the wavelengths of light that would lead to the desired products in a photo-induced reaction/experiment. More concretely, an experimentalist could run this easy-to-use code to predict the specific range of optical frequencies that would be needed to populate a desired photocatalytic product. We envision that further interactions of this kind between quantum control theorists and experimentalists could lead to further advances in site-specific excitations and bond-selective control in a multitude of photocatalytic and photo-excited systems.</p>