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 2 /D 2 isotope effect (on the order of 20) and unexpectedly fast (compared to classical models) decay rate of the singlet oxygen in water; 1 the ring expansion of fluorenylazirines in argon matrix affected by the remote substituents; 2 critical role of the quantum vibronic coupling on the electronic properties and conductivity mechanism in amorphous carbon, 3 and on the hot carrier dynamics in metal halide perovskites. 4 The anharmonic effects on the zero-point energy are seen in the structural and elastic properties of silicon 5 and in cubic silicon carbide. 6 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][10][11][12][13] to superconductivity. 14 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 2 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' 16 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 Schro ¨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, 19 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. 20 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 realtime electronic-structure simulations. Section 3 presents their 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 electronnuclear correlation and the development of the geminal theory of electronic structure.
This journal is © The Royal Society of Chemistry 2026 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
In simulations of soft materials, a hierarchy of moleculardynamics (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][22][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.
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 groundstate potential-energy surface computed on-the-fly from electronic-structure methods, most often density-functional theory (DFT). [24][25][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][28][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 timedependent 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 quantumtrajectory 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 timedependent 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 {F i (r, R)} to be the usual adiabatic eigenstate basis of the electronic Hamiltonian,
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, {c i (R, t)}, satisfy the nuclear time-dependent Schro ¨dinger equation (TDSE). 34 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, c(R, t), solves the TDSE, ĤcðR; tÞ ¼ i h @ @t cðR; tÞ;
where H ˆis the nuclear Hamiltonian operator comprised of the kinetic and potential energy operators, T ˆand V ˆrespectively,
Using here, for simplicity, the Cartesian coordinates for the nuclear positions, and labeling the mass for the nth dimension This journal is © The Royal Society of Chemistry 2026 Chem. Commun., 2026, 62, 7454-7472 | 7457
as m n , the kinetic energy operator is given by
where d is the number of dimensions. Generalizations of eqn ( 3) and ( 4) include the non-adiabatic TDSE, 34 with the time-evolution of c 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 systemspecific. 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 realtime 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.
The most straightforward approach of solving the TDSE (3) is the finite basis representation (FBR). A wavefunction c(R, t) is expanded in a stationary basis of
For the sake of this discussion, we will assume this basis to be orthonormal. Then the Hamiltonian matrix H with the elements H kl = hf k |H ˆ|f l i is constructed and diagonalized,
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 by
Any system properties, such as expectation values, correlation functions or energy eigenstates can be computed from c(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][36][37] addresses point (i): the original basis {f} 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][40][41][42][43]
U ˆ(t) advances a wavefunction from time t = 0 to time t, and the Chebyshev expansion maps application of U ˆ(t) to the initial wavefunction c(R, 0) into a sequence of successive applications of H ˆ. The gist of the SOFT method is in the symmetric splitting of H ˆ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 t 3 . 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 Feature Article ChemComm Open Access Article. Published on 24 March 2026. Downloaded on 4/29/2026 2:11:15 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. This journal is © The Royal Society of Chemistry 2026
Transform algorithm 44 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 5 + by
Wang and Carrington 45 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.
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 46 interfaced with the efficient electronic structure calculations on the fly 47 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, 48 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) Z 0, and phase, S(R, t), cðR; tÞ ¼ AðR; tÞ exp i h 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,
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, UðR; tÞ ¼ À h 2 2mAðR; tÞ
The wavefunction phase along the QT, S t = S(R, t)| R=R t -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) a 0, the quantum potential U, which is proportional to h 2 /m, formally vanishes in the classical limit of heavy particles as h -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,
In the QT frame-of-reference the TDSE yields,
eqn ( 19) describes evolution of r t which conserves the probability within the volume element dR 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,
The weight conservation property allows for practical evaluation of the expectation values. Within the QT framework the initial wavefunction c(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(R, 0). The expectation values of the position-dependent and certain other operators, such as the current density, O ˆ, 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. 19 The approximate U is defined by the fitting of the nonclassical momentum, r R (ln|c|) = r R (ln A), viewed as an additional attribute of a QT. The fitting is performed using a small This journal is © The Royal Society of Chemistry 2026
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). 50 The optimal r ˜yields the following energy-conserving approximate quantum potential,
The corresponding quantum force, F q = Àr R U, is computed analytically. In this approach, calculation of r(R, t) is not needed to perform the dynamics, though its reconstruction is required to obtain the wavefunction c(R, t), or phasedependent correlation functions. Combination of the approximate QT dynamics with on-thefly electronic structure (ES) computed at the density-functional tight-binding (DFTB) level is referred to as 'QTES-DFTB' dynamics. 31 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. 53 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).
The nuclear quantum approaches discussed above target selected nuclear DOFs, typically for relatively small fragments or reduced-dimensionality models. 30 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. 54 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. 55 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 Open Access Article. Published on 24 March 2026. Downloaded on 4/29/2026 2:11:15 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.
where H denotes the Hamiltonian containing kinetic and potential-energy terms, C 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 57 and is non-diagonal when ultrasoft pseudopotentials are used. 60 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 timedependent 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 {f 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 activespace 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 Dt 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 + Dt, is written in terms of a time-evolution operator
where, in principle,
and T denotes the time-ordering operator. In practice, the time-ordered exponential is approximated using a rigorous Magnus expansion,
where the operator X(t, Dt) = X 1 + X 2 + . . . is expressed as a series of integrals over nested commutators of the Hamiltonian matrix evaluated within the interval [t, t + Dt]
Truncation of the Magnus series at any order provides a unitary propagator automatically.
In RMG, the exponential of X 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 picosecondscale simulations with time steps on the order of one atomic unit without noticeable drift in total energy or loss of norm. 55,64 This journal is © The Royal Society of Chemistry 2026
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,
In eqn (33) H 0 (P(t)) is the field-free Kohn-Sham matrix dependent on time only implicitly through the evolving density matrix P(t), m = (m x , m y , m 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) = d(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 broadband 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, timedependent observables are obtained as expectation values of the corresponding operators. In particular, the time-dependent dipole moment
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 velocitygauge 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.
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 quantumtrajectory 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 bioinspired light-harvesting materials.
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.
Understanding the mechanisms of the hydroxide transport inside polyelectrolyte membranes and the factors influencing this process is highly desirable for the rational design of highperforming anion-exchange membranes (AEMs). [65][66][67][68][69][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 2 O + OH À -OH À + H 2 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 hoppinginvolves 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. 16 Recently, the present authors and co-workers 71 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+. 72 Similar on-the-fly electronic-structure MD can be carried out in widely used packages such as CP2K, 73 VASP, 74 and Quantum ESPRESSO, 75 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][77][78][79][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, 81 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) 2 + , 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 Å 3 , 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.
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 Å 2 ps À1 .
To estimate the quantum behavior of the proton, we proceed as follows. The hopping event consists of a proton of H 2 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 timedependence 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 This journal is © The Royal Society of Chemistry 2026 Chem. Commun., 2026, 62, 7454-7472 | 7463
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 0 using Q-Chem. 82 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 onedimensional 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.
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 c(r, 0). The initial wavepacket position, r 0 , and width, a, 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 0 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 {a = 9.30, r 0 = À0.784, p 0 = À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.
This study was motivated by the experiments 83 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. 84 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 B1.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)).
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-oPBEh and M06L paired with 6-31G(d) basis (as implemented in Q-Chem 82 ), 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 This journal is © The Royal Society of Chemistry 2026 Chem. Commun., 2026, 62, 7454-7472 | 7465
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.
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 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 85 ), is estimated as 287.02 and 299.72 mE h within the LRC-oPBEh computational model described above. The NQE on the dipole leads to the energy of the protonic 3-hexylthiophene being lower by 12.7 mE h than its deuteronic counterpart, which is comparable to the NQE on the protonic stabilization of the ZPE in the crystalline phase (D H/ D D crystal/free ZPE) of 25 mE h . These effects of the isotope substitution and isotopic purity on the polymer properties encourage further experimental and computational studies of the NQEs.
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-C61butyric acid methyl ester (PCBM) on the deuteration of the polymer. 86 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. 87 However, as reported in ref. 88 This journal is © The Royal Society of Chemistry 2026 Chem. Commun., 2026, 62, 7454-7472 | 7467
Finally, the effect of the isotope-specific 'noise' on the charge transfer was considered within the theory of Bittner and Silva 89 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 AE 5.7 and and 78.2 AE 4.4 meV for the protonic and deuteronic species, while the coupling between these two states was 74.8 AE 0.21 and 78.2 AE 0.22 meV, respectively. For the B3LYP functional these quantities were 163.1 AE 67.8 and 137.6 AE 31.9 meV for the energy gaps, and 72.2 AE 0.19 and 73.7 AE 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.
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. 55 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 2 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][91][92][93][94] Chlorophylls consist of a rigid, partially hydrogenated porphyrin macrocycle (a chlorin), which forms an extended aromatic p-electron system with a Mg 2+ ion chelated at its center, and a flexible hydrocarbon tail known as phytol. Two nearly degenerate highest-occupied p orbitals and two nearly degenerate lowest-unoccupied p* orbitals generate the strong Soret (or B) band in the blue-UV region (typically B300-450 nm) and the lower-energy Q bands in the red/near-IR range (B500-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.
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 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. 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, m(t) = Tr[P(t)m], 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 wellseparated Q-and Soret-band regions with peak positions and relative intensities consistent with frequency-domain linear-response (lr) TDDFT benchmarks within the randomphase approximation using 6-311+G(d) Gaussian basis (RPA/6-311+G(d)). 57 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 95 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. 96 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, 91 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 chargetransfer 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.
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 quantumthermal bath frameworks combined with on-the-fly electronic structure calculations, which makes description of the highdimensional 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.
This journal is © The Royal Society of Chemistry 2026