Quantum effects on the dynamics and properties of soft materials

Abstract

The quantum effects of nuclear and electronic motion play an important role in the structure, dynamics, and function of soft materials, yet they are difficult to capture with conventional classical simulations or static electronic–structure methods. In this work several complementary approaches for treating quantum effects in polymeric and soft–matter systems are demonstrated, with a focus being on the hydrogen-bonded networks, ion and charge transport, and photoactive chromophores. The proton transfer, tunneling, and isotope effects are captured within the reduced-dimensionality models by implementing grid-based nuclear quantum dynamics in terms of the discrete variable and Fourier bases. The nuclear quantum dynamics is extended to larger systems by employing the quantum trajectories and quantum–thermal bath schemes combined with on-the-fly electronic structure, enabling the description of high-dimensional polymeric environments at feasible cost. The dynamics in the electronic degrees of freedom, simulating the optical response in large chromophores such as chlorophylls, is performed using the real-time time-dependent density functional theory implemented in the real-space multigrid (RMG) code. These approaches are demonstrated on case studies of the proton and hydroxide transport in hydrated polymer membranes, charge transfer in conjugated polymers, and the optical spectra of chlorophyll chromophores relevant to polymerized chlorophyll materials and chlorophyll–polymer hybrids. The reviewed methods and applications highlight practical routes of including quantum effects in simulations of soft functional materials.

---

Sophya Garashchuk

Sophya Garashchuk received her PhD in chemical physics from the University of Notre Dame, Indiana, USA, followed by a postdoctoral appointment at the James Frank Institute at the University of Chicago, Illinois, USA. In 2008 she joined the Department of Chemistry and Biochemistry at the University of South Carolina, USA, as a professor of theoretical chemistry. Her research interests include the nuclear quantum effects on reactivity and properties of molecular systems, method development capturing such effect in the trajectory framework, and computational studies of polymers and hybrid materials.

Jacek Jakowski

Jacek Jakowski received his PhD in theoretical chemistry from the University of Warsaw, Poland, in 2002. He carried out postdoctoral research at the University of Utah, Indiana University, and Emory University. He is currently a staff scientist in the Computational Sciences and Engineering Division at Oak Ridge National Laboratory, USA. His research is focused on the development and application of electronic-structure, first-principles molecular dynamics, and quantum dynamics methods, with emphasis on electronic excitation, non-equilibrium electron dynamics, and nuclear quantum effects in molecules, materials, and soft–matter systems.

Vitaly A. Rassolov

Vitaly Rassolov received his PhD in chemistry from the University of Notre Dame, Indiana, USA, which was followed by a postdoctoral appointment at the Northwestern University, Illinois, USA. He is a professor of theoretical chemistry in the Department of Chemistry and Biochemistry at the University of South Carolina, USA, since 2001. His research interests include conceptual problems of electronic and electron-nuclear correlation and the development of the geminal theory of electronic structure.

---

1 Introduction

Quantum effects in both nuclear and electronic motion play an important role in the properties of molecules and molecular assemblies. The nuclear quantum effects (NQEs) can be loosely divided into ‘hard’ effects which include tunneling and interference involving large-amplitude nuclear motion, and ‘soft’ effects, such as the zero-point energy and electric dipole moment, associated with the effects of delocalized nuclear wavefunctions on the energetics of the chemical bonds and transition states, and on the response properties. In addition, transitions between the electronic states and the ensuing nonequilibrium electronic dynamics are also inherently quantum processes even when the single state dynamics is well approximated by the classical dynamics of the nuclei. Recent research demonstrates significance of the NQEs on the processes involving heavier nuclei. Selected examples include heavy atom tunneling, which accounts for both the unusually large H₂/D₂ isotope effect (on the order of 20) and unexpectedly fast (compared to classical models) decay rate of the singlet oxygen in water;¹ the ring expansion of fluorenylazirines in argon matrix affected by the remote substituents;² critical role of the quantum vibronic coupling on the electronic properties and conductivity mechanism in amorphous carbon,³ and on the hot carrier dynamics in metal halide perovskites.⁴ The anharmonic effects on the zero-point energy are seen in the structural and elastic properties of silicon⁵ and in cubic silicon carbide.⁶ More often, of course, the NQEs are associated with the light species such as hydrogen and helium, with numerous experimental and theoretical studies ranging from properties of water and ice,^7,8 to various inorganic, organic and biological systems,^9–13 to superconductivity.¹⁴ In many of these cases, isotope substitution provides a direct experimental handle on the NQEs.

The most common isotope substitution – deuteration – is a well established technique in characterization of materials (e.g. neutron scattering), and in deducing the reaction mechanisms (e.g. enzymatic reactions). Moreover, deuteration has emerged as a way to tune the properties of complex molecular systems and materials, rather than only to probe them. Opportunities for using D₂O in protein research are discussed in ref. 13, and the role of deuteration in polymeric materials is reviewed in ref. 15. Given the current advances in computational capabilities and methods, including machine-learning based potentials, theoretical modeling and simulations play an increasingly important role in interpreting experiments and isotopic effects, yielding fundamental understanding of reaction mechanisms and enabling predictions of the structural, physical and electronic properties of molecular systems.

Incorporation of the NQEs into molecular dynamics simulations is a long-standing challenge and a very active field of research. Thanks to continuous progress in high-performance computing and methodological developments the NQEs ‘enter mainstream’¹⁶ through path integral molecular dynamics (PIMD) accelerated through ring-polymer contraction, generalized Langevine equation, high-order path integrals and numerous other techniques. The PIMD-type approaches combined with machine-learning techniques accelerating the ab initio electronic energy and force calculations, sampling and data analysis, and with enhanced sampling of rare events enabled studies of the NQEs in biological systems, solid state, liquids and at interfaces.^17,18 Thus, accelerated PIMD simulating the NQEs on the static equilibrium properties of distinguishable particles, though not ‘routine’, is deemed affordable for most systems. Simulation of the dynamical properties, such as tunneling and quantum coherences, however, remains an outstanding challenge. The time-evolution of multidimensional systems, governed by a quantum-mechanical Hamiltonian with non-linear coupling, is characterized by the exponential scaling of the wavefunction complexity with the system size. Recasting the Schrödinger equation into alternative forms, such as those operating with the electronic density or with the correlated trajectory ensemble shifts the exponential scaling to the density functional of the electronic structure theory or to the quantum potential in the Bohmian formulation of quantum dynamics, respectively,¹⁹ and approximations are needed for practical applications. Therefore, development of approximate methods based on restricted wavefunction ansatz and/or restricted interactions, combined with different theory levels – quantum, semiclassical/semiempirical, empirical/classical – to describe various modes of motion or degrees of freedom (DOFs), remains an active area of research.

An overview of the vast field of quantum molecular dynamics is beyond the scope of this article, and we direct an interested reader to a recent special issue ‘Algorithms and software for open quantum system dynamics’ of the Journal of Chemical Physics.²⁰ In the remainder of this paper we review some theoretical approaches and case studies from our groups focused on the nuclear and electronic quantum effects in polymeric and other soft–matter systems. Section 2 describes theoretical methods for the nuclear quantum dynamics, approximate multidimensional quantum treatments, and real-time electronic-structure simulations. Section 3 presents their application to specific experimental systems, including the proton and hydroxide transport in hydrated membranes, charge transfer in conjugated polymers, and the optical response of chlorophyll chromophores. Section 4 provides a summary and an outlook.

2 Theoretical and computational approaches

2.1 Methods overview and preliminaries

In simulations of soft materials, a hierarchy of molecular-dynamics (MD) and sampling schemes is available, differing by how electronic and nuclear degrees of freedom are treated and by the approximations used to trade accuracy for accessible time and length scales (Fig. 1). Classical MD and Monte Carlo simulations based on empirical or coarse-grained force fields as well as machine learned force fields remain the workhorses for exploring structure, thermodynamics, and transport over long times and large system sizes.^21–23 Monte Carlo methods stochastically sample configuration space according to a chosen statistical ensemble, without integrating equations of motion, and are therefore particularly useful for equilibrium properties and free energies.

Fig. 1 Schematic overview of common MD approaches used for soft materials. Classical MD with force fields treats nuclei as classical particles moving on parametrized potentials. In ab initio MD (AIMD/BOMD), classical nuclei are propagated on the ground-state electronic PES obtained on-the-fly from electronic-structure calculations. Electronic quantum dynamics (rt-TDDFT) evolves the electronic state in real time, typically with fixed or classically moving nuclei, and is used to describe excited-state and non-equilibrium electronic processes. Nuclear quantum dynamics approaches treat selected nuclear DOFs as quantum using grid-based or trajectory formulations, coupled to an electronic-structure description or force fields for the remaining DOFs. See text for full description.

When an explicit electronic–structure description is required, ab initio molecular dynamics (AIMD), such as Born–Oppenheimer or Car–Parrinello MD, propagates classical nuclei on a ground-state potential-energy surface computed on-the-fly from electronic-structure methods, most often density-functional theory (DFT).^24–26 Enhanced-sampling techniques, including metadynamics MD, umbrella sampling, and related biasing schemes, are routinely combined with both classical MD and AIMD to accelerate rare events and reconstruct free-energy landscapes along selected collective variables.^27–29 In metadynamics MD, for example, a history-dependent bias potential is built along chosen collective variables to discourage revisiting already explored regions of phase space while gradually converging the underlying free-energy surface.

These approaches treat the nuclei as classical point particles and the electrons as remaining in a single adiabatic electronic state. This framework is adequate for many problems, but it cannot capture nuclear quantum effects or explicitly time-dependent electronic processes. To describe tunnelling, isotope effects, zero-point motion, and other quantum nuclear phenomena, the nuclear degrees of freedom must be treated quantum-mechanically at least along selected coordinates, for example using grid-based wavefunction methods or quantum-trajectory schemes.^1,30,31 Similarly, to simulate non-equilibrium electronic processes and dynamics involving excited states one must introduce explicit time dependence into the electronic structure, as in real-time TDDFT and related time-dependent electronic methods.^32,33 A schematic overview of these different MD schemes and their treatment of electrons and nuclei is shown in Fig. 1.

Now let us turn to the quantum dynamics formalism limited here for simplicity to the time-evolution of wavefunctions. We will consider a wavefunction of the following form,

where r and R are the vectors of coordinates of light and heavy particles, such as the electrons and nuclei, respectively. We take {Φ_i(r, R)} to be the usual adiabatic eigenstate basis of the electronic Hamiltonian,

Ĥ_elΦ_i(r, R) = V_i(R)Φ_i(r, R), (2)

at a fixed nuclear geometry, R. We consider the case of these eigenstates being non-degenerate and real; functions V_i(R) are the Born–Oppenheimer electronic potential energy surfaces (PES). With that, the time-evolution of the R-dependent electronic expansion coefficients, i.e. of the nuclear wavefunctions, {ψ_i(R, t)}, satisfy the nuclear time-dependent Schrödinger equation (TDSE).³⁴

Given the mass- and time-scale separation of the nuclear and electronic motion, the most frequent scenario is that the nuclear wavefunction evolves on a single Born–Oppenheimer PES associated with the ground electronic eigenstate. Omitting the electronic state index, the nuclear wavefunction, ψ(R, t), solves the TDSE,

where Ĥ is the nuclear Hamiltonian operator comprised of the kinetic and potential energy operators, T̂ and V̂ respectively,

Ĥ = T̂ + V(R). (4)

Using here, for simplicity, the Cartesian coordinates for the nuclear positions, and labeling the mass for the nth dimension as m_n, the kinetic energy operator is given bywhere d is the number of dimensions. Generalizations of eqn (3) and (4) include the non-adiabatic TDSE,³⁴ with the time-evolution of ψ_i(R, t) coupled via the second and first derivatives of the electronic wavefunctions with respect to the nuclear positions, and time-dependent potentials, V̂ = V(R, t), due to the external electromagnetic fields.

Formally, the computational efforts of describing a fully coupled anharmoinc quantum system scale exponentially with the system size. Thus, full-dimensional quantum treatment of large-amplitude nuclear motion associated with chemical reactions, isomerizations and highly excited vibrational states, even on a single PES is limited to the systems of 4–5 atoms, or 10–12 degrees of freedom (DOFs). However, most often the nuclei behave as classical or nearly classical particles, and the NQEs are important only for selected DOFs describing light particles, such as protons, at low or moderate temperatures. Therefore, approximate approaches, such as quasi- and semiclassical or mixed quantum/classical nuclear dynamics, and reduced dimensionality computational models, are being developed and employed by researchers for conceptual and practical reasons. The dynamics methods incorporating the NQEs into the studies of large systems tend to be system-specific. Thus, feasibility of evaluating the electronic structure 'on-the-fly' during the propagation of a nuclear wavefunction, as opposed to precomputing the PES prior to the dynamics study, is a major practical factor in choosing an appropriate theoretical/computational approach for a given system.

In the remainder of this section we focus on methods that extend the classical/adiabatic picture. Sections 2.2 and 2.3 discuss conventional basis/grid-based and quantum-trajectory nuclear dynamics that incorporate nuclear quantum effects for selected degrees of freedom, while Section 2.4 describes real-time time-dependent DFT (rt-TDDFT) in a real-space multigrid implementation that provides a fully quantum description of the electronic dynamics on fixed or slowly evolving nuclear configurations. Together, these approaches complement classical and ab initio MD by enabling explicit treatment of tunnelling, isotope effects, and non-equilibrium electronic processes in polymeric and soft–matter systems.

2.2 Dynamics within the finite basis/discrete variable representation of wavefunctions

The most straightforward approach of solving the TDSE (3) is the finite basis representation (FBR). A wavefunction ψ(R, t) is expanded in a stationary basis of N_b functions, {ϕ_k}, k = [1, 2,…, N_b],For the sake of this discussion, we will assume this basis to be orthonormal. Then the Hamiltonian matrix H with the elements H_kl = 〈ϕ_k|Ĥ|ϕ_l〉 is constructed and diagonalized,

HF = FE, (7)

where the columns of the matrix F are the eigenvectors, and the elements E_k of the diagonal matrix E are the corresponding energy eigenvalues. Then, at any t the wavefunction within the FBR is given byAny system properties, such as expectation values, correlation functions or energy eigenstates can be computed from ψ(R, t) in a straightforward manner. Besides the obvious considerations of the basis size and the cost of finding the matrix eigenvectors, one limitation of the FBR is (i) expensive evaluation of N_b(N_b − 1)/2 elements of the potential energy matrix V (representing V in the Hamiltonian of eqn (4)), generally performed by numerical integration. Another limitation is that (ii) the FBR approach implies time-independent PES, V = V(R).

The discrete variable representation (DVR)^35–37 addresses point (i): the original basis {ϕ} in the coordinate space is transformed to make the position operator R diagonal. In the new basis the potential energy matrix V becomes diagonal to a very high accuracy (equivalent to evaluation of the integrals by quadrature associated with the FBR basis underlying the chosen DVR basis). The kinetic energy matrix T is no longer diagonal but it is sparse and its elements are computed analytically (see e.g. ref. 38). After construction of the Hamiltonian matrix in the DVR basis the calculation proceeds as outlined above for the FBR basis.

The second drawback, point (ii) above, is addressed by performing explicit time-evolution of wavefunctions using the split-operator/Fourier Transform (SOFT) method, or the Chebyshev expansion of the time-evolution operator,^39–43

U(τ) = exp(−iĤτ/ħ). (9)

Û(τ) advances a wavefunction from time t = 0 to time τ, and the Chebyshev expansion maps application of Û(τ) to the initial wavefunction ψ(R, 0) into a sequence of successive applications of Ĥ. The gist of the SOFT method is in the symmetric splitting of Ĥ in the exponent into the kinetic and potential energy parts, applied in the momentum and coordinate space, respectively, where T̂ and V̂ are the diagonal operators. The typical scheme,gives the time-propagation error proportional to τ³. The wavefunction is represented in the DVR basis and the diagonal representation of V̂ is invoked to perform the first and last step in the RHS of eqn (10). The wavefunction is Fourier transformed to and from the momentum space to apply the kinetic energy operator at the center of the RHS of eqn (10), which is often facilitated by the application of the Fast Fourier Transform algorithm⁴⁴ scaling as N_b log N_b. This approach does not require matrix operations and is readily applicable to the time-dependent potentials, V = V(R, t).

A state-of-the-art application of the FBR/DVR technique is a calculation of the rotation-bending energy levels of CH₅⁺ by Wang and Carrington⁴⁵ leading to a new assignment of the spectroscopic transitions. The study was enabled by the iterative matrix diagonalization following basis contraction to 50 000 functions and fixing the stretch coordinates. This application illustrates the need for the reduced dimensionality models and for the approximate inclusion of the NQEs into dynamics even at the cost of lower accuracy.

2.3 The quantum trajectory dynamics with approximate quantum potential

Numerical cost of classical dynamics – the time-evolution of trajectories according to the Newton's equations of motion (EOMs) – behind the standard molecular dynamics (MD) simulations scales linearly with the system size, making the trajectory framework appealing for the dynamics of large molecular systems within semiclassical or mixed quantum/classical trajectory-based methods. In this section we describe the approximate quantum-trajectory (QT) dynamics method⁴⁶ interfaced with the efficient electronic structure calculations on the fly⁴⁷ affording incorporation of the dominant NQEs into selected DOFs. Below is a summary of the Madelung–de Broglie–Bohm, or QT, formulation of the TDSE,⁴⁸ using the Cartesian coordinates and a single mass m for all DOFs for simplicity.

A complex time-dependent wavefunction is represented in terms of the real amplitude, A(R, t) ≥ 0, and phase, S(R, t),

Substituting eqn (11) into eqn (3), and switching to the Lagrangian frame-of-reference associated with a trajectory at the position R_t,leads to the following EOMs for the quantum trajectory described by the position and momentum (R_t, P_t):where the vector-function, P(R, t), denotes the gradient of the phase,

P(R, t) = ∇_RS(R, t). (15)

The quantum behavior of the QT comes from the non-local quantum potential, U, dependent on the second derivatives of the wavefunction modulus with respect to the nuclear DOFs,The wavefunction phase along the QT, S_t = S(R, t)|_R=Rt – or in the Lagrangian frame-of-reference defined by eqn (12) – evolves according to the quantum Hamilton-Jacobi equation,Note, that if A(R, t) ≠ 0, the quantum potential U, which is proportional to ħ²/m, formally vanishes in the classical limit of heavy particles as ħ → 0, and eqn (13), (14) and (17) become those of classical mechanics. The form of eqn (16) also suggests a simple way of incorporating the NQEs into the selected DOFs by omitting the ‘classical’ DOFs from the sum.

To complete the QT formulation of the TDSE, the trajectory eqn (13), (14) and (17) are complemented by the EOM on the probability density, ρ,

ρ(R, t) = A²(R, t). (18)

In the QT frame-of-reference the TDSE yields,eqn (19) describes evolution of ρ_t which conserves the probability within the volume element δR_t associated with the trajectory positioned at R_t. As shown for example in ref. 49, this probability, referred to as the trajectory weight w, is constant in time,

w = |ψ(R_t)|²δR_t = |ψ(R₀)|²δR₀. (20)

The weight conservation property allows for practical evaluation of the expectation values. Within the QT framework the initial wavefunction ψ(R, 0) is represented as an ensemble of N_tr trajectories interacting through U. Their momentum is defined according to eqn (15) and their weights are assigned according to the sampling scheme, typically from a random or pseudo-random sampling of ρ(R, 0). The expectation values of the position-dependent and certain other operators, such as the current density, Ô, are computed as sums over the QT ensemble,

Besides the analysis of the dynamics, eqn (21) is central to the construction of a cheap energy-conserving (for time-independent V) approximations to the quantum potential U. Note that in general the cost of exact calculation of U will scale exponentially with the system size, because it is the only non-classical non-local quantity within the QT framework.¹⁹

The approximate U is defined by the fitting of the nonclassical momentum, ∇_R(ln|ψ|) = ∇_R(ln A), viewed as an additional attribute of a QT. The fitting is performed using a small physics-inspired basis, f(R) = (1, R₁, R₂,…, R_d,…) of the size N_b, by minimizing I,

I = 〈ψ|(A⁻¹∇_RA − r̃)²|ψ〉, (22)

with respect to the expansion coefficients of a d-dimensional vector r̃. The elements of this vector are functions defined by the basis expansions,The nonclassical momentum, whose components are linear in R, corresponds to a Gausssian wavefunction, which is an analytic solution to the TDSE with a (time-dependent) parabolic potential.

The expansion coefficients of the kth spatial component are found from the Least Squares Fit with integration by parts: the optimal b_kn are determined by the first and second moments of the QT distribution computed according to eqn (21).⁵⁰ The optimal r̃ yields the following energy-conserving approximate quantum potential,

The corresponding quantum force, F_q = −∇_RU, is computed analytically. In this approach, calculation of ρ(R, t) is not needed to perform the dynamics, though its reconstruction is required to obtain the wavefunction ψ(R, t), or phase-dependent correlation functions.

Combination of the approximate QT dynamics with on-the-fly electronic structure (ES) computed at the density-functional tight-binding (DFTB) level is referred to as ‘QTES-DFTB’ dynamics.³¹ The DFTB is a semi-empirical electronic structure method, which allows a practical quantum-chemical description of bond breaking and reforming on the 100-picosecond time scale for molecular system consisting of a several hundred atoms.^51,52 As a proof-of-principle, the QTES-DFTB dynamics was applied to a model scattering of the hydrogen atom on a graphene flake.⁵³ The protonic wavefunction was represented by an ensemble of QTs shown in Fig. 2(a) colliding with a flexible graphene flake of 37 atoms. The approximate quantum force was included in the DOF corresponding to the normal collision of the proton. The NQEs and the associated H/D isotope dependence of adsorption were assessed by comparing the results of the QTES-DFTB dynamics performed with and without the quantum potential. The results suggest that NQEs can make graphene act as a quantum ‘sieve’ for the H/D separation thanks to a strong preferential absorption of D over H at the collision energy of 0.2 eV illustrated in Fig. 2(b).

Fig. 2 (a) The molecular model of a proton colliding with a graphene flake, with the QT ensemble (grey spheres) collectively representing the protonic wavefunction. (b) Adsorption probability of the hydrogen and deuterium atoms on C₃₇H₁₅ obtained from dynamics with the linearized quantum force (label ‘LQF’) in the collision coordinate added to the classical force, or from the purely classical mechanics (label ‘CM’). In the legend the results for the hydrogen/deuterium systems are labeled as H/D. The initial wave function of C₃₇H₁₅ was represented with multiple (11(14) for H(D)) trajectory ensembles with the graphene atoms moving classically (see Appendix B of ref. 31 for details). Adapted with permission from ref. 31. Copyright 2013 American Chemical Society.

2.4 Electron dynamics employing TDDFT with real-space multigrids (DFT-RMG)

The nuclear quantum approaches discussed above target selected nuclear DOFs, typically for relatively small fragments or reduced-dimensionality models.³⁰ In parallel, many questions in polymer and soft–matter photophysics require an explicit description of the electronic dynamics in large, heterogeneous systems, for example to access optical spectra, exciton localization, and ultrafast charge separation.⁵⁴ For such problems, the real-time time-dependent density functional theory (rt-TDDFT) implemented in the real-space multigrid (RMG) code provides a complementary, fully quantum description of the electrons.⁵⁵

RMG is a DFT-based electronic-structure code^56,57 designed for scalability on massively parallel supercomputers. In RMG, the Kohn–Sham equations^58,59 are discretized on three-dimensional real-space grids and take the form

where H denotes the Hamiltonian containing kinetic and potential-energy terms, Ψ_n is the nth single-electron Kohn–Sham orbital (n = 1,…,N), and S is the overlap matrix. The S matrix reduces to the identity for norm-conserving pseudopotentials⁵⁷ and is non-diagonal when ultrasoft pseudopotentials are used.⁶⁰ An adaptive high-order finite-difference discretization of the Laplacian in the kinetic-energy operator is employed, yielding plane-wave-like accuracy across a wide range of elements, as demonstrated in standard Delta-test benchmarks.^57,61 In practice, the Laplacian is represented by sparse banded matrices on three-dimensional real-space grids, and the finite-difference coefficients are optimized for each chemical composition to minimize the discretization error in atomic reference calculations. The effective potential V_eff in eqn (25) includes the electron–ion interaction (with ions represented by semi-local or nonlocal pseudopotentials), the Hartree electron–electron term, and an exchange–correlation contribution that can incorporate local, semilocal, hybrid, and dispersion-corrected functionals; further numerical details are given in ref. 57.

The rt-TDDFT in RMG provides an efficient framework for modeling electronically non-equilibrium processes in large molecular and soft–matter systems. Within this approach, the electronic dynamics can describe charge-transfer processes, electronically excited states, interaction with external time-dependent fields, optical response, and, when coupled to nuclear motion, non-adiabatic transitions. In contrast to linear-response formulations in the frequency domain, the rt-TDDFT is set up as an initial-value problem: once the initial electronic state and the time dependence of the Hamiltonian are specified, the subsequent evolution of the electronic density matrix is determined by the time-dependent Kohn–Sham equations.

In implementation of ref. 55 the electronic density matrix is propagated in a basis of the Kohn–Sham orbitals obtained from a preliminary ground-state DFT calculation. Let {ϕ_p} denote a set of one-electron orbitals defining the simulation (active) space, typically constructed by combining all occupied Kohn–Sham orbitals with a selected subset of low-lying virtual orbitals to define an excitation window. In this basis, the time evolution of the density matrix P(t) is governed by the Liouville–von Neumann equation

where H(P(t), t) is the time-dependent Hamiltonian matrix that includes the effective one-electron Hamiltonian and any explicit time-dependent external field, and [A, B] = AB − BA denotes the commutator.

At time t = 0 the density matrix P(0) encodes the chosen initial value electronic state. For the optical-response calculations, P(0) is typically taken to be the ground-state density matrix constructed from the occupied Kohn–Sham orbitals,

where C_pn are the orbital expansion coefficients in the active-space basis and the factor of two accounts for spin degeneracy. More generally, excited or non-equilibrium initial states can be prepared by populating selected virtual orbitals, constructing density matrix from linear combinations of Slater determinants, or localizing charge on selected fragments, as was done for charge-transfer states in fullerene collisions.^62,63 In all cases, the rt-TDDFT reduces the problem to propagating eqn (26) with a specified P(0) and time-dependent Hamiltonian.

From a numerical standpoint, time propagation in the rt-TDDFT is more demanding than the ground-state Born–Oppenheimer molecular dynamics because the time step Δt is constrained by the fastest electronic transitions rather than by the nuclear motion. A low-order propagation scheme generates small local errors in P(t) which can accumulate over many time steps and ultimately lead to loss of idempotency, violation of energy conservation, or even exponential divergence of the dynamics. Achieving long-time stability therefore requires a careful balance between accuracy and computational cost: the longer the propagation interval of interest, the more accurate the underlying propagator must be to control error growth.

In RMG, the formal solution of eqn (26) over a time step, from initial time t, to a final time t + Δt, is written in terms of a time-evolution operator

P(t + Δt) = U(t + Δt, t)P(t)U^†(t + Δt, t), (28)

where, in principle,and [scr T, script letter T] denotes the time-ordering operator. In practice, the time-ordered exponential is approximated using a rigorous Magnus expansion,

U(t + Δt, t) = exp[Ω(t, Δt)], (30)

where the operator Ω(t, Δt) = Ω₁ + Ω₂+ … is expressed as a series of integrals over nested commutators of the Hamiltonian matrix evaluated within the interval [t, t + Δt]Truncation of the Magnus series at any order provides a unitary propagator automatically.

In RMG, the exponential of Ω is evaluated using a commutator expansion that is formally equivalent to constructing the evolution operator by exact diagonalization of F, but at a significantly reduced cost. In addition, at each time step the Hamiltonian matrix H(P(t), t) is updated self-consistently with the evolving density, in close analogy to the ground-state SCF procedure. This self-consistency significantly improves the long-time stability of the propagation, enabling picosecond-scale simulations with time steps on the order of one atomic unit without noticeable drift in total energy or loss of norm.^55,64

The interaction of the electronic system with an external field is introduced in the length gauge by adding a dipolar coupling term to the Kohn–Sham Hamiltonian,

H(P(t), t) = H₀(P(t)) − f(t)[small mu, Greek, circumflex]·e. (33)

In eqn (33) H₀(P(t)) is the field-free Kohn–Sham matrix dependent on time only implicitly through the evolving density matrix P(t), [small mu, Greek, circumflex] = (μ_x, μ_y, μ_z) is the vector of dipole-operator matrices in the active-space orbital basis, e is a unit vector specifying the polarization direction, and f(t) is a scalar envelope describing the explicit time dependence of the applied field amplitude. For optical-response calculations, an impulsive “kick” perturbation, f(t) = δ(t), is typically employed to impart a small phase to the occupied Kohn–Sham orbitals. In the frequency domain, such a short pulse corresponds to a broad-band excitation, so a single propagation following the kick contains information about the entire linear absorption spectrum. Alternatively, monochromatic or narrow-band envelopes can be used to selectively excite specific transitions or to drive the system into tailored non-equilibrium states.

Once the density matrix has been propagated, time-dependent observables are obtained as expectation values of the corresponding operators. In particular, the time-dependent dipole moment

μ(t) = Tr[P(t)[small mu, Greek, circumflex]] (34)

is recorded during the propagation and its Fourier transform yields the frequency-dependent polarizability and absorption spectrum. So far, as described in ref. 55), only finite systems are considered, and the length-gauge dipole coupling is employed, which is well suited for molecules and clusters. Extension of the rt-TDDFT implementation in RMG to fully periodic polymeric systems and polymer–solid interfaces will require a velocity-gauge formulation, in which the external field is represented by a time-dependent vector potential that couples to the electronic current operator, and a consistent treatment of nonlocal pseudopotentials in the presence of electromagnetic fields.

In Section 3.4, we illustrate this rt-TDDFT/RMG framework using simulations of the optical response of large chlorophyll chromophores, which serve as prototypical conjugated units and dye moieties embedded in polymeric and soft–matter environments.

2.5 Summary of quantum-dynamical approaches

The methodologies outlined in this Section contribute to a flexible toolkit for treating the quantum effects in polymeric and soft–matter systems across multiple length and time scales. The DVR- and FBR-based nuclear quantum dynamics of Section 2.2 enable exact description of low-dimensional proton transfer, including tunneling, in model environments. The quantum-trajectory and QTES-DFTB schemes of Section 2.3 extend the reach of nuclear quantum dynamics to multidimensional systems and realistic polymer morphologies at a manageable computational cost. Finally, the rt-TDDFT/RMG framework of Section 2.4 provides a quantum dynamical description of electrons in large chromophores and polymer fragments at fixed nuclear geometries, providing access to charge and energy transport and optical response. In Section 3, these complementary approaches are combined and applied to specific case studies, including the hydroxide transport in anion-exchange membranes, the isotope effects on structure and charge transfer in conjugated polymers, and the electronic dynamics of chlorophyll chromophores relevant to polymeric and bio-inspired light-harvesting materials.

3 Applications

First, we present two examples of the NQE assessment based on the reduced dimensionality models, where the nuclear TDSE for a single proton solved using the DVR technique (Sections 3.1 and 3.2). Then, we describe a study where the QTES-DFTB was employed to evaluate delocalization of the protonic/deuteronic wavefunctions of multiple atoms, affecting the electron transfer in a donor/acceptor polymeric system (Section 3.3). Finally, in Section 3.4 we describe a large-scale electron dynamics in model chlorophyll chromophores.

3.1 The quantum dynamical effects during the hydroxide diffusion within the ion exchange membrane

Understanding the mechanisms of the hydroxide transport inside polyelectrolyte membranes and the factors influencing this process is highly desirable for the rational design of high-performing anion-exchange membranes (AEMs).^65–70 At the fundamental level the hydroxide transport occurs via the vehicular and structural diffusion. Formally, the latter process involves proton ‘hopping’ between the water molecule and the hydroxide, H₂O + OH⁻ → OH⁻ + H₂O within the aqueous environment in the presence of a stationary counterion, often terminating a sidechain of a polymeric membrane. In reality, however, the hydroxide motion – with or without hopping – involves rearrangement of 6–7 molecules coordinated to OH⁻, and is affected by the temperature, hydration level, cation identity and spacing, and the membrane composition and morphology. In case of hopping, the NQEs also play a role as they do in pure water.¹⁶

Recently, the present authors and co-workers⁷¹ analyzed the hydroxide transport within the environment of organometallic cobaltocenium-containing copolymers, focusing on the contribution of hopping, absent in the standard molecular dynamics (MD) studies due to the limitations of the classical force-fields. To allow for bond-breaking/bond formation the MD simulations were performed with the electronic structure computed on-the-fly using DFTB+.⁷² Similar on-the-fly electronic-structure MD can be carried out in widely used packages such as CP2K,⁷³ VASP,⁷⁴ and Quantum ESPRESSO,⁷⁵ where ground-state DFT forces are evaluated along classical nuclear trajectories in ab initio MD simulations of liquids, interfaces, and soft materials. Compared with such full-DFT AIMD approaches, DFTB+ offers a substantially lower computational cost while retaining an explicit electronic description, which is crucial for the large polymeric environments considered here. The proton tunneling was estimated from the proton dynamics on the time-dependent PES constructed along the one-dimensional reaction path for a hopping event of the MD trajectory, using the SOFT method of quantum dynamics.

Following the work of Zeldovich et al.^68,69 on the hydroxide diffusion in quaternary ammonium-based AEMs, we developed a simplified molecular model for the cobaltocenium-containing copolymers deemed promising for the AEM applications due to their mechanical, thermal, and chemical stability.^76–80 In our model consisting of 500–600 atoms per simulation cell, the hydrophilic domain is confined by two graphane sheets and includes one space-fixed cobaltocenium, one hydroxide and 10–40 water molecules. Thus, we circumvent the variability of the polymer configurations of the full atomistic description,⁸¹ while controlling the cation spacing, water density, temperature, cobaltocenium orientation and substitutions, all potentially affecting the OH⁻ transport.

To mimic the continuity of the hydrophilic domain of the AEM, the dynamics is performed using a periodic boundary condition imposed on the orthorhombic unit cell of size L_x × L_y × L_z for L_z = 40 Å. A representative simulation cell is shown in Fig. 3. Two hydrophobic graphane sheets parallel to the xy-plane and separated by 10 Å in the z-direction, followed by 30 Å of vacuum. The cell size defines the cation separation. Each cell contains one hydroxide, OH⁻, and one cobaltocenium, Co(Cp)₂⁺, with Co centered between the graphane layers, and 10–40 water molecules. This range corresponds to the number density varying from 7.9 to 31.7 molecules per 1000 ų, or from 23.7% to 95% of the bulk water density. Only atoms of water and hydroxide are treated as ‘movable’. The MD protocol is given in detail in ref. 71.

Fig. 3 The simulation cell of the hydroxide diffusion within a model polymeric environment with the cobaltocenium cation terminating the sidechain. Diffusion of the hydroxide proceeds via vehicular (solid blue arrow) and hopping (green dash arrow) mechanisms, the latter involving the hydrogen transfer (green dotted arrow). Adapted with permission from ref. 71. Copyright 2023 American Chemical Society.

The role of hopping is illustrated in Fig. 4 for the cell size of L_x = 12.26 Å and L_y = 13.36 Å at T = 300 K from 200 ps MD as a function of the number density of water. The largest effect of the hydroxide hopping on the diffusion coefficient is seen at the number density of water equal to about 60% of that for liquid water, when hopping increases the diffusion coefficient from 0.4 to 0.67 Ų ps⁻¹.

Fig. 4 The hydroxide diffusion coefficient D with and without hopping and the average number of hopping events under different hydration levels. Reproduced with permission from ref. 71. Copyright 2023 American Chemical Society.

To estimate the quantum behavior of the proton, we proceed as follows. The hopping event consists of a proton of H₂O transferring to the hydroxide oxygen which defines the reactive coordinate and the corresponding PES. Since the reactants and products are the same, the fully relaxed PES is a symmetric double well, and the proton transfer is reversible, unless the product configuration is stabilized by the changes in the molecular environment. Within a simple reaction path model, based on one of the hopping events from our MD simulations, we represented the effect of the environment by the time-dependence of the PES defined along the proton transfer coordinate (one spatial dimension). The oxygen atoms remain nearly stationary on the time scale of the proton hop, which for the chosen MD trajectory is about 100 fs. The time coordinate, parameterizing the environmental configuration effectively serves as the second dimension. The model limitations are: it is based on the classical MD trajectory; the quantum dynamics is performed only along the proton transfer coordinate, i.e., there are no quantum corrections to other protonic modes of motion.

The PES is constructed based on 120 fs of dynamics surrounding the proton hop at 70 fs of one of the MD trajectories. The changing environment was represented by 13 snapshots taken at 10 fs intervals, and the molecular representation was truncated to 20 water molecules. The potential energy curves were computed for the proton transfer between the donor and acceptor oxygens, labeled O^D and O^A respectively, along the reaction coordinate r, defined as

r = | HO^D| − |HO^A|. (35)

The energy curves are obtained by performing the partial energy minimization with respect to the positions of the hydrogen atoms attached to the donor and acceptor oxygens while scanning the position of the transferring proton between the donor and acceptor oxygens. The electronic structure calculations are performed at the B3LYP/6-31+G(d,p) level in gas phase for r in the range of r = [−1.5,1.5] a₀ using Q-Chem.⁸² The PES is constructed as the cubic spline fit; at the edges the PES curves are extrapolated as quadratic functions of r. Finally, the cubic spline is used to interpolate across the one-dimensional curves along the time-coordinate. The resulting PES, dependent on the reaction coordinate and time, is shown in Fig. 5(a). The reactant/product regions correspond to the positive/negative r-values.

Fig. 5 The proton transfer model: (a) the contour plot of the time-dependent potential energy surface; the reaction coordinate r and time are given on the vertical and horizontal axes, respectively (in a.u.) The initial wavepacket is centered at r₀ = 0.784a₀. (b) The reactant well populations, computed using exact quantum time-evolution, and classical dynamics of the trajectories sampling the Wigner distribution corresponding to the initial wavefunction. The vertical green dot-dashes mark the interval of the hop in the MD trajectory underlying the developed PES. Adapted with permission from ref. 71. Copyright 2023 American Chemical Society.

To estimate the quantum effects, we propagate a protonic wavefunction, initialized as the Gaussian wave packet,

using SOFT method on a grid, and compare the reactant well probabilities to the results of the classical evolution of 2000 classical trajectories, sampling the Wigner distribution corresponding to the initial wavepacket ψ(r, 0). The initial wavepacket position, r₀, and width, α, correspond to the vibrational ground state within the local harmonic approximation to the PES at the minimum of the reactant well. The initial momentum p₀ is directed towards the acceptor oxygen, and the corresponding translational energy is equivalent to 300 Kelvin. The numerical parameter values in atomic units (a.u.) are {α = 9.30, r₀ = −0.784, p₀ = −1.868} a.u. The reactant well probabilities are shown in Fig. 5(b). At the beginning, both classical and quantum simulations give similar reaction probability. Within 1000 a.u. of time, 25% of the reactant goes to the product side. Then, however, the proton transfer is reversed and at 1500 a.u. of time, 95% of the probability density is on the reactant side. This behavior correlates with the disappearing and reappearing of the PES reactant well induced by the environment. This quasi-oscillatory behavior persists over the following 2000 a.u. of time during which a shallow reactant well develops. After 3500 a.u. of time 60% of the quantum and 30% of the classical nuclei fall into the product well, yielding the ratio of the quantum to classical transfer probability of about 2. Overall, the proton hop becomes ‘irreversible’ after stabilization at a favorable configuration of the environment corresponding to the shortened distance between the donor and acceptor oxygens. The difference between the quantum and classical probabilities is attributed to the non-locality of the protonic wavefunction.

3.2 Isotope effect on stability of the polymeric poly(3-hexylthiophene)

This study was motivated by the experiments⁸³ performed on the polymeric poly(3-hexylthiophene) (P3HT) at four levels of deuteration – pristine (P-P3HT), main-chain deuterated (MD-P3HT), side-chain deuterated (SD-P3HT), and fully deuterated (FD-P3HT) – sketched in Fig. 6(a). According to the small angle neutron scattering (SANS) and wide-angle X-ray scattering (WAXS) measurements the deuteration of the thiophene ring (MD-P3HT and FD-P3HT), but not on the side-chain (SD-P3HT) reduced the crystallinity of P3HT. The WAXS data are shown in Fig. 6(b), where the deuteration number of 0, 1, 13 and 14, refers to P-P3HT, MD-P3HT, SD-P3HT and FD-P3HT, respectively. These results suggest that the difference in crystallinity is due to the quantum behavior of a single hydrogen/deuterium on the thiophene ring, and that the change of the dipole–dipole interactions, which largely determines stacking of the polymeric chains, may play a role.⁸⁴

Fig. 6 Crystallinity of protonated and deuterated P3HT. (a) Molecular structures of protonated and deuterated P3HT: pristine (P-P3HT), main chain (thiophene) deuterated (MD-P3HT), side-chain (hexyl) deuterated, fully deuterated (hexyl and thiophene). (b) X-ray (100) diffraction peak intensity of P3HT as a function of different numbers of deuterium atoms per P3HT monomer. The numbers of deuterium atoms per monomer for P-P3HT, MD-P3HT, SD-P3HT, and FD-P3HT are 0, 1, 13, and 14, respectively. The (100) peak intensity is proportional to the sample crystallinity. Reproduced with permission from ref. 83. Copyright 2017 American Chemical Society.

To understand the stability trend among the four P3HT isotopologues, the dependence of the zero-point energy (ZPE) and of the dipole moment on the isotope was examined computationally and theoretically.^83,85 The atomistic model, based on the crystallography data consists of four chains of P3HT (408 atoms total) shown in Fig. 7(a) with all but one chain omitted for clarity. The extent and orientation of all the sidechains are shown in Fig. 7(b). Each chain, shown in Fig. 7(c), has four hexylthiophene units (25 atoms) and is terminated with the hydrogen atoms. According to the electronic structure calculations 3-hexylthiophene exhibits an appreciable dipole moment of ∼1.1 Debye, largely localized on the thiophene ring and oriented in its plain. In crystalline P3HT the neighboring polymeric chains are arranged according to the parallel–displaced stacking of the thiophene rings, so that the dipole moments are in staggered anti-parallel orientation, which minimizes the dipole–dipole interactions (Fig. 7(a)).

Fig. 7 The model of the P3HT crystal consists of four P3HT chains of four thiophene hexyl units each. Relative orientation of the dipole moments, shown as blue arrows, leads to the anti-parallel stacking. Yellow/red spheres indicate the sulfur/hydrogen or deuterium on the thiophene ring; the remaining hydrogen atoms and carbons are shown in white and light blue spheres, respectively. (a) The outer layer (“MM region”) is represented by point charges from the DFT calculation. The inner layer (“QM region”) consists a single 3-hexylthiophene unit represented via DFT embedded in the external field of point charges in the MM region. Single proton or deuteron is treated as a quantum particle described via DVR; the DVR region is within the red box. (b) View along the polymer chains with the sidechains extended. (c) A close up of the chain with the quantum H/D. The DVR grid with the superimposed contours of the actual potential energy are shown in purple. Panels (a) and (b) are reproduced with permission from ref. 85. Copyright 2018 John Wiley and Sons.

First, let us focus on the ZPE effect on the stability of the crystalline model within the harmonic oscillator description of the bond vibrations, routinely computed by the electronic structure codes. The ZPE was calculated for the hydrogen and deuterium in crystalline P3HT and for a single isolated (or free-space) chain of P3HT employing three different electronic structure methods: LRC-ωPBEh and M06L paired with 6-31G(d) basis (as implemented in Q-Chem⁸²), and the DFTB method. The empirical dispersion correction D3 was used in all calculations. The vibrational frequencies of the thiophene CH/CD bond, obtained with all three methods, show that the crystalline environment flattens the PES where H/D connects to the thiophene ring, compared to the PES of an isolated P3HT chain. Thus, the force constants and the ZPEs are reduced for both H and D. The ZPE lowering is by about 50% more pronounced for the protonic species, as schematically depicted in Fig. 8, which contributes to the stabilization of the protonic crystallized species. The corresponding ZPEs are presented in Table 1. Within the harmonic approximation there is no isotope effect on the dipole moment, which is a constant.

Fig. 8 The ZPE level diagram: crystal field decreases the force constants, thus lowering the ZPE for both the hydrogen and the deuterium species leading to increased stability of the crystalline P3HT compared to an isolated P3HT chain. Reproduced with permission from ref. 83. Copyright 2017 American Chemical Society.

Table 1 The ZPE per deuterium or hydrogen obtained within the harmonic approximation to the PES for crystalline P3HT (ZPE_crystal) and for a single chain of P3HT (ZPE_free). All values are in mE_h. Adapted with permission from ref. 83. Copyright 2017 American Chemical Society

Method	Deuterium			Hydrogen
	ZPEcrystal	ZPEfree	ΔZPE	ZPEcrystal	ZPEfree	ΔZPE
LRC-ωPBEh-D3	8.184	8.243	−0.058	11.569	11.652	−0.083
M06L-D3	7.726	8.020	−0.294	10.922	11.337	−0.415
DFTB-D3	7.724	8.408	−0.683	10.919	11.885	−0.966

---

More accurate estimates of the NQEs associated with the isotope substitutions in P3HT are made using DVR^36,38 to obtain the energy eigenstates. The computational model is composed of three layers treated at different levels of theory (Fig. 7(a)). (i) The outer layer is the molecular mechanics (MM) region, where the interatomic interactions are modeled via the point charges. (ii) The inner layer is the quantum mechanics (QM) region where the electronic structure is described from the first principles via DFT to provide an adequately accurate PES for the proton/deuteron motion. (iii) The DVR region of the proton/deuteron quantum dynamics. The PES and the protonic/deuteronic wave function of the quantum nucleus are computed on a three-dimensional Cartesian grid. This approach includes the effects of the PES anharmonicity on the nuclear wave function, and gives the variance of the dipole moment on the grid.

(i) Within the MM region (Fig. 7(a) and (b)) all atoms are represented as point charges, defined by the Mulliken charges from the electronic structure calculation employing the long-range corrected (LRC) density functional with the valence basis, i.e. LRC-ωPBEh/6-31G(d). (ii) The intermediate-scale QM region, consisting of the 3-hexylthiophene unit (27 atoms) from the central part of the large-scale molecular model (Fig. 7(c)), is imposed on the MM region. The electronic structure calculations, performed using larger basis (LRC-ωPBEh-D3/6-31++G(d,p) method), incorporate the electrostatic interactions with the surrounding point charges from the MM layer. This is the model used to generate the PES and the dipole moment surfaces of the small-scale DVR region. (iii) The nuclear wave function for the quantum H/D is represented on a three-dimensional Cartesian grid of 31 × 41 × 41 = 52 111 points spaced by 0.05 Å in all directions (Fig. 7(c)). The axis of the grid are aligned with the vibrational modes of the quantum H/D: C–H stretch, C–H bend in the plane with the thiophene ring and the out-of-plane bend. The Hamiltonian constructed within the grid representation using five-point DVR is diagonalized to obtain the nuclear wavefunctions, used to evaluate the dipole moment and the transition dipole moments (see ref. 85 for full details).

The dependence of the dipole moment on the three mode coordinates is shown in Fig. 9. The in-plane bend leads to the largest changes in the magnitude of the dipole moment, while the out-of plane bend changes its orientation. The anisotropy of the dipole moment combined with the PES anharmonicity leads to non-zero transition dipole moments. Table 2 lists several dipole moment integrals, computed over the protonic and deuteronic wavefunctions. Overall, the total dipole moment vector of the P3HT unit consists of a “stationary” dipole moment of a thiophene ring independent of H/D vibrations, and of its instantaneous fluctuations associated with the motion of H/D quantified by their root-mean square (RMS) for the ground H/D vibrational states. As expected, the dipole moment fluctuations are larger for H than for D: RMS_H = [2.9,11.0,3.3] mDeb, RMS_D = [1.9,7.8,2.2] mDeb. Their role on the interchain dipole–dipole interaction, averaged over a delocalized wavefunction included (the first order perturbation theory⁸⁵), is estimated as 287.02 and 299.72 μE_h within the LRC-ωPBEh computational model described above. The NQE on the dipole leads to the energy of the protonic 3-hexylthiophene being lower by 12.7 μE_h than its deuteronic counterpart, which is comparable to the NQE on the protonic stabilization of the ZPE in the crystalline phase (Δ_H/DΔ_crystal/free ZPE) of 25 μE_h. These effects of the isotope substitution and isotopic purity on the polymer properties encourage further experimental and computational studies of the NQEs.

Fig. 9 Variation of the dipole magnitude and of the ground state wave function for H and D as a function of displacement from the equilibrium position along the normal mode direction. The nuclear wave function was calculated within the DVR approach on the ab initio PES computed at the LRC-ωPBEh-D3 theory level. Adapted with permission from ref. 83. Copyright 2017 American Chemical Society.

Table 2 The dipole moment (k = 0 in 〈0|x|k〉) and the transition dipole moments (k > 0) associated with the motion of H/D on the thiophene ring in P3HT. The vibrational energy, E(k), is listed in mE_h. The quantum numbers {l, m, n} are associated with the vibrational modes of C–H, i.e. stretch, out-of-plane bend and in-plane bend, respectively. Adapted with permission from ref. 85. Copyright 2018 John Wiley and Sons

k	lmn	Deuterium				Hydrogen
		E(k)	〈0|x|k〉	〈0|y|k〉	〈0|z|k〉	E(k)	〈0|x|k〉	〈0|y|k〉	〈0|z|k〉
0	000	8.10997	9.70 × 10^−1	−2.61 × 10^−2	−5.13 × 10^−1	11.43831	9.63 × 10^−1	−2.57 × 10^−2	−5.18 × 10^−1
1	001	10.54077	−1.62 × 10^−2	−8.79 × 10^−2	3.32 × 10^−3	14.86221	−1.95 × 10^−2	−1.04 × 10^−1	3.89 × 10^−3
2	010	11.95388	2.94 × 10^−2	8.85 × 10^−3	−4.44 × 10^−2	16.85546	3.60 × 10^−2	1.07 × 10^−2	−5.27 × 10^−2
9	001	17.83907	−2.41 × 10^−2	5.46 × 10^−4	−1.64 × 10^−2	25.07164	−2.95 × 10^−2	8.96 × 10^−4	−2.07 × 10^−2

---

3.3 The isotope effect on charge transfer in P3HT/PCBM

The QTES-DFTB approach was employed to gain insight into the unexpected dependence of the open circuit voltage, V_OC, in a blend of poly(3-hexylthiophene) (P3HT) and[6,6]-phenyl-C61-butyric acid methyl ester (PCBM) on the deuteration of the polymer.⁸⁶ On the fundamental level, for the organic donor/acceptor system, such as P3HT/PCBM, the efficiency of the current generation upon the UV irradiation is determined by the difference in electrochemical potentials, charge transfer integrals, and electron/hole mobilities defined by their electronic structure. The output voltage correlates with V_OC, and is thus associated with the difference of the energy levels between the highest occupied molecular orbital (HOMO) of an electron donor and the lowest unoccupied molecular orbital (LUMO) of an electron acceptor.⁸⁷ However, as reported in ref. 88 in P3HT/PCBM the H/D substitution on the sidechains reduced experimental V_OC by 0.24 V, while such substitution on the polymer backbone had very small effect. Compared to the pristine P3HT V_OC was lowered only 0.01 V in the latter case. These observations that could not be explained by the changes in the electronic structure motivated application of the QTES-DFTB approach to describe the protonic/deuteronic wavefunctions on the sidechains.

In this theoretical isotope effect study,⁸⁶ the P3HT isotopologues are the same as shown in Fig. 6(a), and the molecular model consists of 221 atoms of PCBM enveloped by the polymer (Fig. 10). The 39 highlighted protons (or deuterons) were treated as quantum particles, their wavefunctions represented by an ensemble of 4300 QTs. The heavy atoms of the sidechains moved according to the classical forces averaged over the quantum nuclei, while the remaining atoms were space-fixed. The quantum force was computed independently for each proton/deuteron from the fitting of the non-classical momentum in a quadratic basis to capture the anharmonicity of the PES. The average CH/CD bond-distances were obtained from 500 fs of QTES-DFTB dynamics (the charge transfer proceeds on the time-scale of 100 fs). Average hydrogen-carbon bond lengths for CH and CD at the end and in the middle of the sidechains are shown in Fig. 11.

Fig. 10 The molecular model of P3HT/PCBM consisting of 221 atoms: carbon/hydrogen/sulfur/oxygen = gray/white/yellow/red spheres, respectively. 39 hydrogens, highlighted in blue, are treated as quantum particles. Reproduced with permission from ref. 86. Copyright 2016 American Chemical Society.

Fig. 11 Average hydrogen-carbon bond length for the hexyl arms in P3HT/PCBM in the course of the QTES-DFTB dynamics as a function of time. Top panel: Internal H/D, located within –(CH₂)₅– of hexyl groups. Bottom panel: Terminal H/D, located within methyl ends of hexyl groups. Adapted with permission from ref. 86. Copyright 2016 American Chemical Society.

Then, the electronic energy levels associated with the charge transfer and charge separation, and the couplings were computed along the selected QTs using four (BLYP, B3LYP, CAM-B3LYP and LRC-ωPBEh) density functionals. As shown in Fig. 12 upon the photo-excitation the electrons moves from the HOMO localized on P3HT (the donor) to the LUMO+3 spanning both the donor and acceptor, which is coupled to LUMO+2 localized on PCBM (the acceptor) and also to LUMO+1 (not shown) and LUMO of the same character. Regardless of the DFT functional, the isotope-specific effect on the HOMO–LUMO gap and on the couplings between the electronic states was found negligible. For example, the TDDFT gaps, E_g, for the protonated(deuterated) species computed with the B3LYP, CAM-B3LYP and LRC-ωPBEh functionals were equal to E_g = 1.16(1.17), 2.22(2.23) and 2.346(2.353) eV, respectively. Thus, the variation in E_g did not explain the experimental observation for the V_OC lowering by 0.24 eV upon the side-chain deuteration.

Fig. 12 Frontier orbitals of P3HT-PCBM from CAM-B3LYP/6-31G(d). The initial excitation proceeds from the donor-localized HOMO to a charge transfer LUMO+3 coupled to the charge separated LUMO+2, LUMO+1 (not shown) and LUMO. Adapted with permission from ref. 86. Copyright 2016 American Chemical Society.

Finally, the effect of the isotope-specific ‘noise’ on the charge transfer was considered within the theory of Bittner and Silva⁸⁹ arguing that the noise may promote fast charge transfer via tunneling from a neutral state to a charge-separated state. As detailed in ref. 86 within the P3HT/PCBM model the charge transfer proceeds from LUMO+3 to LUMO+2, LUMO+1 and LUMO. The fluctuations of the energy gap between the charge transfer and charge separated states due to the nuclear wave functions were found larger for the protonic than for the deuteronic P3HT, and comparable to the coupling between the states. For example, within CAM-B3LYP/6-31G(d) method the LUMO+3/LUMO+2 gap was 74.9 ± 5.7 and and 78.2 ± 4.4 meV for the protonic and deuteronic species, while the coupling between these two states was 74.8 ± 0.21 and 78.2 ± 0.22 meV, respectively. For the B3LYP functional these quantities were 163.1 ± 67.8 and 137.6 ± 31.9 meV for the energy gaps, and 72.2 ± 0.19 and 73.7 ± 0.20 meV for the couplings, within the protonic and deuteronic P3HT, respectively. Based on the charge transfer probabilities within the four-state model, it was demonstrated that such isotope-dependence of the Hamiltonian elements may lead to appreciable changes on the charge transfer probability. Hence the isotope effect could account for the experimental trends by promoting the charge transfer in P3HT:PCBM and increasing the charge recombination on the donor in the deuterium-substituted P3HT:PCBM.

3.4 Chlorophylls as model chromophores for electron dynamics and optical response simulations

In this Section we complement the nuclear quantum dynamics of Section 3.1–3.3 by the quantum-dynamical description of electrons in chlorophylls a and b, as representative macrocylcic chromophores, and model their optical response using the rt-TDDFT/RMG method.⁵⁵

Chlorophylls are bioorganic pigments found in green plants, algae, and photosynthetic microorganisms such as cyanobacteria and marine phytoplankton. They play a central role in photosynthesis by harvesting sunlight and converting it into electronic excitation followed by photoinduced charge separation, which ultimately drives CO₂ reduction and carbohydrate synthesis. Their strong visible-light absorption, rich photophysics, and metal–macrocycle redox chemistry make chlorophylls attractive as sustainable, non-toxic components for photocatalysis, energy conversion, optoelectronics, sensing, and biomedical applications.^90–94

Chlorophylls consist of a rigid, partially hydrogenated porphyrin macrocycle (a chlorin), which forms an extended aromatic π-electron system with a Mg²⁺ ion chelated at its center, and a flexible hydrocarbon tail known as phytol. Two nearly degenerate highest-occupied π orbitals and two nearly degenerate lowest-unoccupied π* orbitals generate the strong Soret (or B) band in the blue–UV region (typically ∼300–450 nm) and the lower-energy Q bands in the red/near-IR range (∼500–700 nm). The exact positions, splittings, and intensities of these Soret and Q bands are controlled by the metal center and asymmetry introduced by substituents in the chlorin macrocycle. Among several natural variants, the most abundant are chlorophyll a and chlorophyll b, differing by a methyl versus formyl substituent and shown in Fig. 13 along with their UV-vis absorption spectra.

Fig. 13 The rt-TDDFT/RMG simulations of chlorophylls a and b. (a) Molecular structure of chlorophyll b used in the calculations, with the formyl (–CHO) group highlighted (orange dotted oval); in chlorophyll a this position is occupied by a methyl (–CH₃) group. The color scheme for the atoms: H = white, C = light blue, O = red, N = dark blue, Mg = green. (b) Total energy for both chlorophylls during the rt-TDDFT propagation, illustrating excellent energy conservation. (c) Comparison of photoabsorption spectra from the rt-TDDFT/RMG and the linear-response TDDFT (RPA/6-311+G(d)) calculations. The RPA spectra are broadened with Gaussians to mimic experimental line shapes. (d) Experimental photoabsorption spectra of chlorophylls a and b. Reproduced with permission from ref. 55. Copyright 2025 American Chemical Society.

While chlorophylls are highly efficient light and energy harvesters in biological systems, their direct use as electrode materials in supercapacitors and related electrochemical devices is hampered by the poor electrical conductivity of monomeric chlorophyll molecules. To address this limitation, electropolymerization of suitably functionalized chlorins was used to prepare conjugated poly(chlorophyll) films.^91,92 Electropolymerized chlorophyll electrodes, including Ni-based poly(chlorophyll) systems, exhibit enhanced capacitance and improved charge-transport properties compared to their monomeric counterparts, while retaining broad visible absorption and redox activity. Hybrid devices that integrate a chlorophyll-based photoactive layer with a supercapacitor element show stable photocharging–discharging behaviour under illumination, highlighting the potential of chlorophyll-derived polymers as non-toxic, degradable, and renewable components for combined solar-energy conversion and storage.⁹¹ In parallel, the photoredox and photosensitizing capabilities of chlorophylls and their derivatives motivate their use as visible-light photoredox catalysts and as photosensitizers for singlet-oxygen generation in antimicrobial photodynamic therapy and in chlorophyll-loaded polymer matrices, fibres, and coatings with antimicrobial activity.^90,93,94 Together, these developments position chlorophyll-based and chlorophyll-containing polymers as promising platforms for environmentally friendly energy, catalytic, and biomedical materials.

In the rt-TDDFT/RMG study of ref. 55, isolated chlorophyll a (Chl-a) and chlorophyll b (Chl-b) are considered as representatives of macrocyclic chromophores that can be embedded in or converted into polymeric architectures. Chlorophyll geometries were first optimized at a conventional DFT level using a Gaussian basis, then placed in large orthorhombic simulation cells with sufficient vacuum to eliminate interactions with periodic images. Within RMG (see disucssion in Section 2.4), the Kohn–Sham orbitals are represented on real-space grids and the PBE functional and norm-conserving pseudopotentials are employed. The active space for the rt-TDDFT includes all occupied Kohn–Sham orbitals and a set of low-lying virtual orbitals spanning the low-energy π → π* excitation window relevant for the Q- and Soret-band manifolds.

Starting from the converged ground-state, described by the density matrix P(0), a weak impulsive kick along a selected polarization direction was applied and density matrix P(t) was propagated in time according to the Liouville–von Neumann equation using the Magnus-expansion-based propagator of Section 2.4. The kick amplitude was chosen so that the response remains in the linear regime. The time-dependent dipole moment, µ(t) = Tr[P(t)[small mu, Greek, circumflex]], was recorded and the absorption spectrum was obtained from its Fourier transform. The excellent long-time numerical stability of the propagation, with no noticeable drift in the total energy over the duration of the simulation, provides a stringent test of the rt-TDDFT/RMG algorithm for these large, anisotropic chromophores.

The computed spectra for Chl-a and Chl-b are shown in Fig. 13 along with the experimental ones. They exhibit well-separated Q- and Soret-band regions with peak positions and relative intensities consistent with frequency-domain linear-response (lr) TDDFT benchmarks within the random-phase approximation using 6-311+G(d) Gaussian basis (RPA/6-311+G(d)).⁵⁷ The rt-TDDFT peak positions agree closely with the RPA results, with only small differences in relative intensities. Experimental spectra for the chlorophylls in diethyl ether⁹⁵ show Soret bands with maxima near 428 nm (Chl-a) and 453 nm (Chl-b) and Q bands near 661 nm (Chl-a) and 642 nm (Chl-b). The TDDFT spectra are obtained in vacuum for a single geometry and therefore neglect thermal motion and solvent broadening, but the level of agreement is comparable to previous CAM-B3LYP TDDFT studies.⁹⁶ The simulated spectra reproduce the Q-band region with peak maxima close to 640 nm and the stronger Q-band peak for Chl-a relative to Chl-b. In the blue region, the simulations resolve multiple overlapping features of the experimental Soret bands.

The rt-TDDFT simulations of the chlorophylls provide a convenient platform to explore how local electrostatic environments and chemical modification of the macrocycle affect the optical response of polymerized chlorophyll materials and chlorophyll–polymer hybrids. Within the RMG framework, one can readily introduce explicit nearby fragments to mimic the effect of a polymer matrix, a supporting electrode, or adjacent chromophores. Changes in the splitting and intensity distribution of the Q- and Soret-band features upon such perturbations can be related to shifts in local site energies and couplings, providing microscopic input to exciton and charge-transfer models for extended chlorophyll assemblies. For example, in the context of electropolymerized PolyChl films and chlorophyll-based photoelectrodes,⁹¹ one may probe how polymerization-induced changes in the macrocycle environment or chain conformation influence the absorption cross section, the polarization dependence of the optical response, and the distribution of charge following photoexcitation.

More broadly, the chlorophyll calculations demonstrate that the rt-TDDFT implemented in RMG can treat electronically complex chromophores at a size and level of detail comparable to those encountered in polymerized chlorophyll systems, chlorophyll-doped polymer networks, and chlorophyll-loaded fibres.^90,92,93 While the present examples focus on isolated molecules in vacuum, the same methodology can be combined with ensembles of nuclear configurations generated from classical or quantum molecular dynamics in a polymeric environment. The nuclear quantum effects and thermal disorder in polymer matrices can be mapped onto distributions of chlorophyll excitation energies, transition dipoles, and charge-transfer pathways, bridging the nuclear quantum phenomena discussed in Section 3.1–3.3 with the electronic response of polymerized chlorophyll materials and other bioinspired polymer–chromophore architectures.

4 Summary and outlook

To summarize, in this article we discussed several complementary theoretical and modeling approaches of incorporating the quantum behavior affecting structure, dynamics and function in polymeric and soft–matter systems. The presented studies included hydrogen-bonded networks, ion and charge transport, and photoactive chromophores. For selected nuclear degrees of freedom we employed grid-based nuclear quantum dynamics using discrete variable and Fourier bases to describe proton transfer, tunneling, and isotope effects in reduced-dimensionality models. The nuclear quantum dynamics was extended to larger systems by employing the quantum trajectory and quantum–thermal bath frameworks combined with on-the-fly electronic structure calculations, which makes description of the high-dimensional polymeric environments computationally feasible. For the electronic degrees of freedom we presented a simulation of the electronic dynamics and optical response in large chromophores, i.e. chlorophylls, using the real-time time-dependent density functional theory implemented in the real-space multigrid (RMG) code. These methods were illustrated on several case studies, including the proton and hydroxide transport in hydrated polymer membranes, charge transfer in conjugated polymers, and the optical spectra of chlorophyll chromophores, as relevant to polymerized chlorophyll materials and chlorophyll–polymer hybrids. These examples highlight practical routes of including the quantum effects in simulations of soft functional materials. Further development of theoretical approaches and computational tools, in particular those bypassing the eigenstate representation for both nuclei and electrons and incorporating hierarchical description of 'quantumness' of the nuclei, is highly desirable for facilitation of the fundamental studies and practical development of the polymeric and soft materials.

Author contributions

Conceptualization: SG, JJ and VR; methodology: SG, JJ, VR; funding and computational resource acquisition: SG, JJ; writing original draft: SG, JJ, VR; review & editing: SG, JJ, VR.

Conflicts of interest

There are no conflicts to declare.

Data availability

No new data were presented in this work.

Acknowledgements

This material is based upon work supported by the National Science Foundation of USA under Awards CHE-2308922, CHE-1056188, CHE-1565985, and OIA-1655740, and by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences Separation Science program under Awards DE-SC0020272 and DE-SC0025229. JJ acknowledges support by the Laboratory Directed Research and Development Funding, Non-Equilibrium And emergent Transients (NEAT) Initiative, Oak Ridge National Laboratory. Computations were performed in part using allocation DMR110037 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by the grants #2138259, #2138286, #2138307, #2137603, and #2138296 awarded by the National Science Foundation of U.S.A. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory (project #MAT739), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This research was conducted at Oak Ridge National Laboratory’s Center for Nanophase Materials Sciences (CNMS), which is a DOE Office of Science User Facility.

Notes and references

1. I. M. Ansari, E. R. Heller, G. Trenins and J. O. Richardson, Heavy-Atom Tunnelling in Singlet Oxygen Deactivation Predicted by Instanton Theory with Branch-Point Singularities, Nat. Commun., 2024, 15, 4335,  DOI:10.1038/s41467-024-48463-2.
2. J. F. Rowen, F. Beyer, T. Schleif and W. Sander, Isomer-Specific Heavy-Atom Tunneling in the Ring Expansion of Fluorenylazirines, J. Org. Chem., 2023, 88(13), 7893–7900,  DOI:10.1021/acs.joc.3c00484.
3. A. Kundu, Y. Song and G. Galli, Influence of nuclear quantum effects on the electronic properties of amorphous carbon, Proc. Natl. Acad. Sci. U. S. A., 2022, 119(31), e2203083119,  DOI:10.1073/pnas.2203083119.
4. Y. Liu, S. Shen, O. V. Prezhdo, R. Long and W.-H. Fang., Nuclear Quantum Effects Accelerate Hot Carrier Relaxation but Slow Down Recombination in Metal Halide Perovskites, J. Am. Chem. Soc., 2025, 147(13), 11543–11554,  DOI:10.1021/jacs.5c02139 PMID: 40106363.
5. D. S. Kim, O. Hellman, J. Herriman, H. L. Smith, J. Y. Y. Lin, N. Shulumba, J. L. Niedziela, C. W. Li, D. L. Abernathy and B. Fultz, Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon, Proc. Natl. Acad. Sci. U. S. A., 2018, 115(9), 1992–1997,  DOI:10.1073/pnas.1707745115.
6. C. P. Herrero, R. Ramírez and G. Herrero-Saboya, Nuclear quantum effects in structural and elastic properties of cubic silicon carbide, Phys. Rev. B, 2024, 109, 104112,  DOI:10.1103/PhysRevB.109.104112.
7. K.-Y. Chiang, J. Hunger, M. Bonn and Y. Nagata, Experimental quantification of nuclear quantum effects on the hydrogen bond of liquid water, Sci. Adv., 2025, 11(14), eadv7218,  DOI:10.1126/sciadv.adv7218.
8. M. L. Berrens, A. Kundu, M. F. Calegari Andrade, T. Anh Pham, G. Galli and D. Donadio, Nuclear Quantum Effects on the Electronic Structure of Water and Ice, J. Phys. Chem. Lett., 2024, 15(26), 6818–6825,  DOI:10.1021/acs.jpclett.4c01315. PMID: 38916450.
9. I. Savchenko, B. Gu, T. Heine, J. Jakowski and S. Garashchuk, Nuclear quantum effects on adsorption of H₂ and isotopologues on metal ions, Chem. Phys. Lett., 2017, 670, 64–70,  DOI:10.1016/j.cplett.2016.12.069.
10. L. Gem Gao, R. Ming Zhang, X. Xu and D. G. Truhlar, Quantum Effects on H₂ Diffusion in Zeolite RHO:Inverse Kinetic Isotope Effect for Sieving, J. Am. Chem. Soc., 2019, 141(34), 13635–13642,  DOI:10.1021/jacs.9b06506.
11. A. C. Thakur and R. C. Remsing, Nuclear quantum effects in the acetylene:ammonia plastic co-crystal, J. Chem. Phys., 2024, 160(2), 024502,  DOI:10.1063/5.0179161.
12. H. E. Sauceda, V. Vassilev-Galindo, S. Chmiela, K.-R. Mueller and A. Tkatchenko, Dynamical strengthening of covalent and non-covalent molecular interactions by nuclear quantum effects at finite temperature, Nat. Commun., 2021, 12, 442,  DOI:10.1038/s41467-020-20212-1.
13. G. Giubertoni, M. Bonn and S. Woutersen, D₂O as an Imperfect Replacement for H₂O: Problem or Opportunity for Protein Research?, J. Phys. Chem. B, 2023, 127(38), 8086–8094,  DOI:10.1021/acs.jpcb.3c04385.
14. S. van de Bund and G. J. Ackland, Competition between superconductivity and molecularization in the quantum nuclear behavior of lanthanum hydride, Phys. Rev. B, 2023, 108, 184102,  DOI:10.1103/PhysRevB.108.184102.
15. L. Li, J. Jakowski, C. Do and K. Hong, Deuteration and Polymers: Rich History with Great Potential, Macromolecules, 2021, 54(8), 3555–3584,  DOI:10.1021/acs.macromol.0c02284.
16. T. E. Markland and M. Ceriotti, Nuclear quantum effects enter the mainstream, Nat. Rev. Chem., 2018, 2, 0109,  DOI:10.1038/s41570-017-0109.
17. B. E. Ugur and M. A. Webb, Nuclear quantum effects in molecular liquids across chemical space, Nat. Commun., 2025, 16, 5786,  DOI:10.1038/s41467-025-60850-x.
18. D. Ray and M. Parrinello, Kinetics from Metadynamics: Principles, Applications, and Outlook, J. Chem. Theory Comput., 2023, 19(17), 5649–5670,  DOI:10.1021/acs.jctc.3c00660.
19. V. A. Rassolov and S. Garashchuk, Computational complexity in quantum chemistry, Chem. Phys. Lett., 2008, 464(4), 262–264,  DOI:10.1016/j.cplett.2008.09.026.
20. A. Chin, J. Keeling, D. Segal and H. Wang, Algorithms and software for open quantum system dynamics, J. Chem. Phys., 2025, 163(5), 050401,  DOI:10.1063/5.0289390.
21. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017, ISBN 9780198803195 DOI:10.1093/oso/9780198803195.001.0001.
22. D. Frenkel and B. Smit, Understanding Molecular Simulation, 2nd edn, Academic Press, 2002 DOI:10.1016/B978-012267351-1/50000-6.
23. O. T. Unke, S. Chmiel, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Machine Learning Force Fields, Chem. Rev., 2021, 121(16), 10142–10186,  DOI:10.1021/acs.chemrev.0c01111.
24. C. Ribaldone and S. Casassa, Born–Oppenheimer Molecular Dynamics with a Linear Combination of Atomic Orbitals and Hybrid Functionals for Condensed Matter Simulations Made Possible. Theory and Performance for the Microcanonical and Canonical Ensembles, J. Chem. Theory Comput., 2024, 20(9), 3954–3975,  DOI:10.1021/acs.jctc.3c01231.
25. T. D. Kuhne, Second generation Car-Parrinello molecular dynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 391–406.
26. D. Marx and J. Hutter, Ab Initio Molecular Dynamics:Basic Theory and Advanced Methods, Cambridge University Press, 2009.
27. D. Ray, N. Ansari, V. Rizzi, M. Invernizzi and M. Parrinello, Rare Event Kinetics from Adaptive Bias Enhanced Sampling, J. Chem. Theory Comput., 2022, 18(11), 6500–6509.
28. D. Ray and M. Parrinello, Kinetics from Metadynamics:Principles, Applications, and Outlook, J. Chem. Theory Comput., 2023, 19(17), 5649–5670,  DOI:10.1021/acs.jctc.3c00660 PMID: 37585703.
29. K. Zhu, E. Trizio, J. Zhang, R. Hu, L. Jiang, T. Hou and L. Bonati, Enhanced Sampling in the Age of Machine Learning: Algorithms and Applications, Chem. Rev., 2026, 126(1), 671–713.
30. S. Garashchuk, J. Huang, B. G. Sumpter and J. Jakowski, From classical to quantum dynamics of atomic and ionic species interacting with graphene and its analogue, in Properties and Functionalization of Graphene, of Theoretical and Computational Chemistry, ed. T. Dinadayalane and F. Hagelberg, Elsevier, 2022, vol. 21, pp. 61–86 DOI:10.1016/B978-0-12-819514-7.00001-4.
31. S. Garashchuk, J. Jakowski, L. Wang and B. G. Sumpter, Quantum Trajectory-Electronic Structure Approach for Exploring Nuclear Effects in the Dynamics of Nanomaterials, J. Chem. Theory Comput., 2013, 9, 5221–5235,  DOI:10.1021/ct4006147.
32. T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg and S. Tretiak, Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials, Chem. Rev., 2020, 120(4), 2215–2287,  DOI:10.1021/acs.chemrev.9b00447 . PMID: 32040312.
33. X. Li, N. Govind, C. Isborn, III, A. Eugene DePrince and K. Lopata, Real-Time Time-Dependent Electronic Structure Theory, Chem. Rev., 2020, 120, 9951–9993,  DOI:10.1021/acs.chemrev.0c00223.
34. J. C. Tully, Nonadiabatic dynamcs, in Modern methods for multidimensional dynamics computations in chemistry, ed. D. Thompson, World Scientific, 1998, pp. 34–72.
35. R. W. Heather and J. C. Light, Discrete Variable Theory of Triatomic Photodissociation, J. Chem. Phys., 1983, 79, 147.
36. J. C. Light, I. P. Hamilton and J. V. Lill, Generalized Discrete Variable Approximation in Quantum Mechanics, J. Chem. Phys., 1985, 82, 1400.
37. J. V. Lill, G. A. Parker and J. C. Light, The Discrete Variable-Finite Basis Approach to Quantum Scattering, J. Chem. Phys., 1986, 85, 900.
38. D. T. Colbert and W. H. Miller, A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method, J. Chem. Phys., 1992, 96(3), 1982–1991,  DOI:10.1063/1.462100.
39. D. Kosloff and R. Kosloff, A Fouirier method solution for the time-dependent Schrodinger equation as a tool in molecular dynamics, J. Comput. Phys., 1983, 52(1), 35–53,  DOI:10.1016/0021-9991(83)90015-3.
40. C. Leforestier, et al., A comparison of different popagation schemes for the time-dependent Schrodinger equation, J. Comput. Phys., 1991, 94(1), 59–80,  DOI:10.1016/0021-9991(91)90137-A.
41. R. Kosloff, The Fourier Method, in Numerical grid methods and their applications to Schrodinger equation, ed. C. Cerjan, 1993, vol. 412, pp. 175–194, of NATO Advanced Science Institutes Series C, NATO, SCI COMM, 1993. NATO Advanced Research Workshop on Grid Methods in Atomic and Molecular Quantum Calculations, Corte, Fracne, Sep 27-Octs 03, 1992.
42. R. Kosloff, Propagation methods for quantum molecular dynamics, Ann. Rev. Phys. Chem., 1994, 45, 145–178,  DOI:10.1146/annurev.pc.45.100194.001045.
43. S. R. Billeter and W. F. VanGunsteren, A comparison of different numerical propagation schemes for solving the timedependent Schrodinger equation in the position representation in one dimension for mixed quantum- and molecular dynamics simulations, Mol. Simul., 1995, 15(5), 301–322,  DOI:10.1080/08927029508022343.
44. M. D. Feit, J. A. Fleck and A. Steiger, Solution of the Schrödinger equation by a spectral method, J. Comput. Phys., 1982, 47, 412.
45. X.-G. Wang and T. Carrington Jr., Calculated rotation-bending energy levels of CH + 5 and a comparison with experiment, J. Chem. Phys., 2016, 144(20), 204304,  DOI:10.1063/1.4948549.
46. S. Garashchuk and V. A. Rassolov, Semiclassical Dynamics with Quantum Trajectories: Formulation and Comparison with the Semiclassical Initial Value Representation Propagator, J. Chem. Phys., 2003, 118, 2482–2490.
47. S. Garashchuk and M. V. Volkov, The energy-conserving dynamics of quantum-classical systems based on quantum trajectories, Mol. Phys., 2012, 110, 985–993,  DOI:10.1080/00268976.2012.675449.
48. D. Bohm, A suggested interpretation of the quantum theory in terms of “hidden” variables, I and II, Phys. Rev., 1952, 85, 167.
49. S. Garashchuk and V. A. Rassolov, Energy conserving approximations to the quantum potential: Dynamics with linearized quantum force, J. Chem. Phys., 2004, 120, 1181–1190.
50. S. Garashchuk and V. A. Rassolov, Quantum dynamics with Bohmian trajectories: Energy conserving approximation to the quantum potential, Chem. Phys. Lett., 2003, 376, 358–363.
51. M. Elstner and G. Seifert, Density Functional Tight Binding, Philos. Trans. R. Soc. Lond., B, Biol. Sci., 2014, 372(2011), 20120483.
52. M. Gaus, Q. Cui and M. Elstner, DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB), J. Chem. Theory Comput., 2011, 7(4), 931–948.
53. Quantum Trajectory Electronic Structure (code), 2016, https://github.com/garashch/qtes-code.
54. D. Lingerfelt, P. Ganesh, B. G. Sumpter and J. Jakowski, From ground to excited electronic state dynamics of electron and ion irradiated graphene nanomaterials, in Properties and Functionalization of Graphene, of Theoretical and Computational Chemistry, ed. T. Dinadayalane and F. Hagelberg, Elsevier, 2022, vol. 21, pp. 87–107 DOI:10.1016/B978-0-12-819514-7.00003-8.
55. J. Jakowski, W. Lu, E. Briggs, D. Lingerfelt, B. G. Sumpter, P. Ganesh and J. Bernholc, Simulation of 24 000 Electrons Dynamics: Real-Time Time-Dependent Density Functional Theory (TDDFT) with the Real-Space Multigrids (RMG), J. Chem. Theory Comput., 2025, 21, 1322–1339.
56. E. L. Briggs, D. J. Sullivan and J. Bernholc, Real-space multigrid-based approach to large-scale electronic structure calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54(20), 14362–14375.
57. E. L. Briggs, W. Lu and J. Bernholc, Adaptive finite differencing in high accuracy electronic structure calculations, npj Comput. Mater., 2024, 10, 17,  DOI:10.1038/s41524-024-01203-y.
58. P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev., 1964, 136, B864–B871.
59. W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev., 1965, 140, A1133–A1138.
60. D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 7892,  DOI:10.1103/PhysRevB.41.7892.
61. K. Lejaeghere, et al., Reproducibility in density functional theory calculations of solids, Science, 2016, 351(6280), aad3000,  DOI:10.1126/science.aad3000.
62. J. Jakowski, S. Irle, B. Sumpter and K. Morokuma, Modeling charge transfer in fullerene collisions via real-time electron dynamics, J. Phys. Chem. Lett., 2012, 3, 1536–1542.
63. J. Jakowski, S. Irle and K. Morokuma, Time-dependent quantum dynamical simulations of C2 condensation under extreme conditions, Phys. Chem. Chem. Phys., 2012, 14, 6273–6279.
64. J. Jakowski and K. Morokuma, Liouville-von Neumann molecular dynamics, J. Chem. Phys., 2009, 130, 224106.
65. W. Zhang, D. Dong, D. Bedrov and A. C. T. Van Duin, Hydroxide Transport and Chemical Degradation in Anion Exchange Membranes: A Combined Reactive and Non-Reactive Molecular Simulation Study, J. Mater. Chem. A, 2019, 7(10), 5442–5452.
66. G. E. Lindberg, C. Knight, R. Jorn, J. F. Dama and G. A. Voth, Multiscale Simulation of Hydroxide Solvation and Transport in Anion Exchange Membranes, ECS Trans., 2011, 41(1), 1785–1793.
67. G. E. Lindberg, C. Knight, L. E. Felberg and G. A. Voth, Molecular Dynamics Simulations of Hydroxide Solvation and Transport in Anionic Exchange Membranes, ECS Trans., 2013, 50(2), 2053.
68. T. Zelovich, L. Vogt-Maranto, M. A. Hickner, S. J. Paddison, C. Bae, D. R. Dekel and M. E. Tuckerman, Hydroxide Ion Diffusion in Anion-Exchange Membranes at Low Hydration: Insights from Ab Initio Molecular Dynamics, Chem. Mater., 2019, 31(15), 5778–5787.
69. T. Zelovich, Z. Long, M. Hickner, S. J. Paddison, C. Bae and M. E. Tuckerman, Ab Initio Molecular Dynamics Study of Hydroxide Diffusion Mechanisms in Nanoconfined Structural Mimics of Anion Exchange Membranes, J. Phys. Chem. C, 2019, 123(8), 4638–4653.
70. C. Chen, Y.-L. Steve Tse, G. E. Lindberg, C. Knight and G. A. Voth, Hydroxide Solvation and Transport in Anion Exchange Membranes, J. Am. Chem. Soc., 2016, 138(3), 991–1000.
71. S. Wickramasinghe, A. Hoehn, S. T. Wetthasinghe, H. Lin, Q. Wang, J. Jakowski, V. Rassolov, C. Tang and S. Garashchuk, Theoretical Examination of the Hydroxide Transport in Cobaltocenium-Containing Polyelectrolytes, J. Phys. Chem. B, 2023, 127(47), 10129–10141,  DOI:10.1021/acs.jpcb.3c04118.
72. B. Hourahine, et al., DFTB+, a Software Package for Efficient Approximate Density Functional Theory Based Atomistic Simulations, J. Chem. Phys., 2020, 152(12), 124101.
73. T. Kühne, et al., CP2K: an electronic structure and molecular dynamics software package – Quickstep: Efficient and accurate electronic structure calculations, J. Chem. Phys., 2020, 150, 194103,  DOI:10.1063/5.0007045.
74. J. Hafner, Ab-initio simulations of materials using VASP:Density-functional theory and beyond, J. Comput. Chem., 2008, 29(13), 2044–2078,  DOI:10.1002/jcc.21057.
75. P. Giannozzi, et al., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys. Condens. Matter, 2009, 21, 395502.
76. T. Zhu, S. Xu, A. Rahman, E. Dogdibegovic, P. Yang, P. Pageni, M. Pabel Kabir, X.-D. Zhou and C. Tang, Cationic Metallo-Polyelectrolytes for Robust Alkaline Anion-Exchange Membranes, Angew. Chem., 2018, 130(9), 2412–2416.
77. T. Zhu, Y. Sha, J. Yan, P. Pageni, M. Anisur Rahman, Y. Yan and C. Tang, Metallo-Polyelectrolytes as a Class of Ionic Macromolecules for Functional Materials, Nat. Commun., 2018, 9(1), 1–15.
78. T. Zhu, Y. Sha, H. Adabi Firouzjaie, X. Peng, Y. Cha, D. M. M. Mevan Dissanayake, M. D. Smith, A. K. Vannucci, W. E. Mustain and C. Tang, Rational Synthesis of Metallo-Cations Toward Redox-and Alkaline-Stable Metallo-Polyelectrolytes, J. Am. Chem. Soc., 2019, 142(2), 1083–1089.
79. W. Yang, S. Liu, J. Yan, F. Zhong, N. Jia, Y. Yan and Q. Zhang, Metallo-Polyelectrolyte-Based Robust Anion Exchange Membranes via Acetalation of a Commodity Polymer, Macromolecules, 2021, 54(19), 9145–9154.
80. L. Ren, C. G. Hardy, S. Tang, D. B. Doxie, N. Hamidi and C. Tang, Preparation of Side-Chain 18-E Cobaltocenium-Containing Acrylate Monomers and Polymers, Macromolecules, 2010, 43(22), 9304–9310.
81. D. Dong, W. Zhang, A. C. T. Van Duin and D. Bedrov, Grotthuss Versus Vehicular Transport of Hydroxide in Anion-Exchange Membranes: Insight from Combined Reactive and Nonreactive Molecular Simulations, J. Phys. Chem. Lett., 2018, 9(4), 825–829.
82. Y. Shao, et al., Advances in molecular quantum chemistry contained in the Q-Chem 4 program package, Mol. Phys., 2015, 113, 184–215,  DOI:10.1080/00268976.2014.952696.
83. J. Jakowski, J. Huang, S. Garashchuk, Y. Luo, K. Hong, J. Keum and B. G. Sumpter, Deuteration as a Means to Tune Crystallinity of Conducting Polymers, J. Phys. Chem. Lett., 2017, 8, 4333–4340,  DOI:10.1021/acs.jpclett.7b01803.
84. D. V. Talapin, E. V. Shevchenko, C. B. Murray, A. V. Titov and P. Král, Dipole-Dipole Interactions in Nanoparticle Superlattices, Nano Lett., 2007, 7(5), 1213–1219,  DOI:10.1021/nl070058c . PMID: 17397231.
85. J. Jakowski, J. Huang, B. G. Sumpter and S. Garashchuk, Theoretical assessment of the nuclear quantum effects on polymer crystallinity via perturbation theory and dynamics, Int. J. Quantum Chem., 2018, 118, e25712,  DOI:10.1002/qua.25712.
86. L. Wang, J. Jakowski, S. Garashchuk and B. G. Sumpter, Understanding How Isotopes Affect Charge Transfer in P3HT/PCBM: A Quantum Trajectory-Electronic Structure Study with Nonlinear Quantum Corrections, J. Chem. Theory Comput., 2016, 12, 4487–4500,  DOI:10.1021/acs.jctc.6b00126.
87. R. Kroon, M. Lenes, J. C. Hummelen, P. W. M. Blom and B. de Boer, Small Bandgap Polymers for Organic Solar Cells (Polymer Material Development in the Last 5 Years), Polym. Rev., 2008, 48(3), 531–582,  DOI:10.1080/15583720802231833.
88. M. Shao, et al., The isotopic effects of deuteration on optoelectronic properties of conducting polymers, Nat. Commun., 2014, 5(1), 3180.
89. E. R. Bittner and C. Silva, Noise-induced quantum coherence drives photo-carrier generation dynamics at polymeric semiconductor heterojunctions, Nat. Commun., 2014, 5, 3119,  DOI:10.1038/ncomms4119.
90. J.-Q. Cai, X.-M. Liu, Z.-J. Gao, L.-L. Li and H. Wang, Chlorophylls derivatives: Photophysical properties, assemblies, nanostructures and biomedical applications, Mater. Today, 2021, 45, 77–92,  DOI:10.1016/j.mattod.2020.11.001.
91. C. Zhang, S. Duan, M. Zhou, Z. Liu, H. Ren, S. Ichi Sasaki and X.-F. Wang, Electropolymerized chlorophyll derivative biopolymers for supercapacitors, Chem. Eng. J., 2022, 450, 138000,  DOI:10.1016/j.cej.2022.138000.
92. S. Shanmugam, J. Xu and C. Boyer, Utilizing the electron transfer mechanism of chlorophyll a under light for controlled radical polymerization, Chem. Sci., 2015, 6, 1341,  10.1039/c4sc03342f.
93. P. M. Tyubaeva, et al., Electrospinning of Poly-3-Hydroxybutyrate Fibers Loaded with Chlorophyll for Antibacterial Purposes, Polymers, 2024, 16, 3221,  DOI:10.3390/polym16223221.
94. J. M. R. Narayanam and C. R. J. Stephenson, Visible light photoredox catalysis: applications in organic synthesis, Chem. Soc. Rev., 2011, 40, 102–113,  10.1039/B913880N.
95. H. K. Lichtenthaler and C. Buschmann, Chlorophylls and Carotenoids: Measurement and Characterization by UV-VIS Spectroscopy, Curr. Protocols Food Anal. Chem., 2001, 1(1), F4.3.1–F4.3.8,  DOI:10.1002/0471142913.faf0403s01.
96. B. F. Milne, Y. Toker, A. Rubio and S. Brøndsted Nielsen, Unraveling the Intrinsic Color of Chlorophyll, Angew. Chem., 2015, 127, 2198–2201,  DOI:10.1002/ange.201410899.

---

Footnote

† This manuscript has been authored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The publisher acknowledges the US government license to provide public access under the DOE Public Access Plan (https://energy.gov/downloads/doe-public-access-plan).

---