Transition metal oxides have attracted much attention as photo(electrochemical)-catalysts but practical applications are typically hampered by their low and anisotropic charge mobility. A deep understanding of excess charge carrier transport in these materials requires a dynamical treatment of nuclear motion that goes well beyond standard approaches. Here we introduce DeepPolaron, a machine learning framework boosting the accessible time scale of first principles molecular dynamics of adiabatic polaron transport by three orders of magnitude at a virtually negligible loss in accuracy. We apply our method to excess electron and hole transport in titanium dioxide rutile and anatase. We find that the excess electron in rutile relaxes to a polaron predominantly localized on a single Ti atom with hopping occurring only along the [001] direction, associated with an activation energy of 39 meV and a room temperature mobility of 4.4 × 10^-2 cm^2/Vs in good agreement with experiment. In contrast the hole polaron in anatase is localized on a single O atom, and due to poor O 2p orbital overlap with first nearest neighbors charge transport occurs primarily to second nearest neighbors, with a large activation energy of 139 meV resulting in a small room temperature mobility of 1.4 × 10^-3 cm^2/Vs. This work provides a finite temperature first-principles characterization of small polaron transport in rutile and anatase, with a methodology that is directly transferable to other small polaron forming materials and interfacial charge-transfer processes.

High-mobility organic molecular crystals such as rubrene are important materials for organic electronics, yet a quantitatively predictive description of their charge transport properties remains challenging. Direct mixed quantum-classical nonadiabatic molecular dynamics simulations provide a promising route by explicitly propagating the charge carrier wavefunction, without assuming a specific transport mechanism. However, previous large-scale simulations of apolar molecular crystals have commonly neglected dynamic electrostatic disorder, since evaluation of electrostatic interactions is computationally demanding and the approximation appears plausible for apolar systems. Here, we use the damped shifted-force (DSF) real-space electrostatic summation method, combined with an efficient addition-subtraction scheme, to include dynamic electrostatic disorder in fragment orbital-based surface hopping (FOB-SH) simulations of room-temperature hole transport in rubrene. We find that electrostatic interactions increase the reorganization energy for (hypothetical) nearest neighbour hopping by 29 and 39 meV along the a and b-directions, respectively, relative to a baseline value of 152 meV obtained without electrostatics. In FOB-SH simulations, electrostatic interactions lead to increased site energy disorder, reducing the spatial extent of the hole wavefunction, as measured by a decrease in the inverse participation ratio from 13 to 9, and lowering the predicted mobility along the high-mobility direction from 35 to 21 cm^2 V^-1s^-1, in close agreement with experiment.

Recent experiments have shown contradictory effects of strong light-matter coupling on exciton-exciton annihilation (EEA) in organic molecular systems. In this work, we perform numerical simulations of polariton dynamics and reveal the role of strong coupling in changing the EEA rate. The results of our simulations suggest that strong coupling allows to partially overcome disorder in systems with poor exciton mobility via delocalisation of excitons owing to the interaction with the common cavity mode. This leads to an enhanced connectivity between excitons and, consequently, to an increase in the EEA rate. Conversely, in systems with high exciton mobility, in which disorder has a much smaller effect on excitation energy transfer, excitons can interact strongly even without coupling to the cavity photons at the exciton densities at which EEA typically occurs. In this case, the EEA rate can be even lower than in bare molecules due to the existence of a competing decay channel associated with photon leakage through the cavity mirrors. We also find that in the weak coupling regime, the EEA rate appears to be suppressed due to this decay channel regardless of the exciton transport properties. Our simulations resolve the experimental controversy on the effect of strong coupling on EEA and provide guidance for minimising the EEA rate towards a more feasible realisation of Bose-Einstein condensation of polaritons.

Cable bacteria are multicellular bacteria capable of centimeter-scale conduction through a regular fiber network embedded in their cell envelope. The conductivity of these fibers is extremely high for biological materials, and rivals that of the best synthetic conductive polymers, but the underlying electron transport mechanism remains elusive. Recent microscopic and spectroscopic evidence indicates that each fiber embeds a bundle of intertwined nanoribbons as the conductive conduit. Each nanoribbon consists of a one-dimensional nickel-organic framework, built from stacked nickel bis(1,2-dithiolene) oligomers (NiBiD units) as molecular building blocks. Here we performed DFT calculations of nanoribbon model structures, in order to characterize their electronic properties, examine potential stacking configurations and verify whether these structures can support efficient conductance. Our simulations indicate that nanoribbons are comprised of tightly stacked AA or AB-type packings of NiBiD units. In the most energetically stable structure (AB-type) some Ni centers are predicted to be 5-fold coordinated due to formation of an inter-layer Ni-S coordination bond. In several energetically low-lying structures, the electronic coupling between neighboring molecules exceeds the critical threshold for charge delocalization permitting efficient charge transport beyond small polaron hopping. Our results hence reveal that nanoribbons based on NiBiD units exhibit favorable charge transfer properties that may explain the unusually high conductivities measured in the fibers of cable bacteria.

Multiheme cytochromes have emerged as a promising new class of bioelectronic materials for potential use in soft and biocompatible electronics. To aid these developments it is essential to understand how conductance decays with increasing size of multiheme proteins and which parts of the proteins govern their high conductance. Here we present a systematic investigation of six native multiheme cytochromes forming electronic junctions of increasing size. We find that the computed conductances are best described by a protein-specific intercept model with a single exponential distance decay constant of β = 1.8 nm-1 shared across all six proteins and an increasing intercept with protein size. Moreover, our calculations suggest that conductance for larger multiheme cytochromes tends to be less sensitive to thermal protein fluctuations than for smaller proteins. Both observations can be explained by the higher density of protein electronic states with increasing protein size. Fe ions appear to be unimportant for electronic conduction and the contribution of heme-cofactors is relatively small compared to their essential role for electron transfer in native biological environment. Our calculations suggest that it is mainly the protein frame and the contacts between amino acids and electrodes that govern electronic conductance in the series of multiheme cytochromes investigated.

A molecular understanding of the solvation and dynamics of ions under static electric fields is crucial for modeling a wide range of natural and technological processes. Yet, traditional simulation methods suffer from a trade-off that has to be made between accuracy and statistical convergence. To bridge this gap, herein, we extend our recently introduced perturbed neural network potential molecular dynamics (PNNP MD) approach to investigate the solvation structures and ionic transport mechanisms of electrified alkali cationic solutions. We obtain ionic conductivities for Li+, Na+, and Cs+ from the field dependence of the ionic current density in good agreement with experiment. Surprisingly, the migration mechanism is found to be strikingly different for the three ions, despite their similar ionic conductivities. While Li+ conducts predominantly through vehicular migration of a stable 4-fold coordinated ion at all field strengths, Cs+ conducts strictly through a structural diffusion mechanism, where 9-12 transient first shell water coordination bonds are continuously broken and reformed. Notably, aqueous Na+ emerges as a "Goldilocks" ion: its ion-water interactions are strong enough to maintain distinct 5-6-fold coordination shells at zero field (unlike Cs+) yet labile enough to be strongly perturbed by electric fields (unlike Li+). As a consequence, we observe an electric-field-induced transition from vehicular to structural ionic transport for Na+ that is accompanied by a marked increase in ionic current density. Our results imply that the conductance mechanism of ions with moderate ion-solvent interactions can be effectively tuned by external electric fields.

We have recently introduced excitonic state-based surface hopping, a powerful non-adiabatic molecular dynamics method capable of simulating charge photogeneration in nanoscale molecular systems. Here, we assess the impact of the initial conditions, treatment of thermal fluctuations, and model dimensionality when simulating exciton dissociation at an organic donor-acceptor interface. We find that realistic modeling of the initial photoexcitation into the partially delocalized excitonic band states is essential to capture the short-term relaxation dynamics that are often observed in experiment. Yet, the population dynamics on longer timescales remain largely insensitive to the initialization procedure. Two-dimensional systems are found to exhibit a significantly increased density of hybrid Frenkel exciton-charge transfer states at the same excitation energy when compared to one-dimensional systems, resulting in accelerated exciton decay. Moreover, the increase in the ratio of charge transfer to Frenkel exciton states alongside dimensionality entropically favors charge generation. Finally, we find that the impact of neglecting the thermal fluctuations of electronic couplings depends on the specific parameter regime of the interface. In the incoherent hopping regime, their neglect leads to only a very small decrease in the exciton decay rate. In the transient delocalization regime, the delocalization of charge transfer states away from (or close to) the interface is overestimated (or underestimated), resulting in a markedly different exciton decay profile when thermal fluctuations are neglected.

One of the distinguishing aspects of CP2K is its seamless integration of diverse structural and transition-state optimization techniques with advanced sampling approaches including Monte Carlo, molecular dynamics, and metadynamics, enabling the efficient exploration of complex potential- and free-energy landscapes, including rare events. These capabilities are combined with a broad hierarchy of energy and force evaluation methods, ranging from classical and machine-learned interaction potentials and mixed quantum-classical multiscale and semiempirical schemes, to highly accurate quantum-mechanical electronic-structure approaches. At the heart of the latter lies the Gaussian and plane-wave framework, along with its augmented all-electron generalization, which have been described in detail in our previous code review [T. D. Kühne et al., J. Chem. Phys. 152, 194103 (2020)]. Building on this foundation, the present work revisits the methods within CP2K that turn electronic structure into dynamics, transport, and spectroscopic response. Particular emphasis is placed on the coupling between static response calculations and nuclear motion: spectra may be evaluated at optimized structures, averaged over thermally sampled configurations, obtained from time-correlation functions along ab-initio or path integral molecular trajectories, or followed in real time together with electronic and nuclear dynamics. The same modular structure also enables equilibrium and biased transport simulations, from Kubo-type linear response to open-boundary approaches under external potentials, highlighting CP2K's unique capability to unify quantum chemistry with quantum and statistical mechanics within a versatile, holistic simulation environment.

Detailed balance, the correct thermalization of electronic state populations, is an essential and elusive property of quantum–classical non-adiabatic dynamics methods. While some methods can reproduce detailed balance through physically well-motivated algorithmic adaptations, or by construction of a conserved Hamiltonian function, the physical mechanism leading to detailed balance is not understood from first principles. Coupled trajectory mixed quantum–classical (CTMQC) dynamics may provide some insight into the question, as it can be derived from first principles in the exact factorization theorem of full quantum mechanics. Although we find that the current conventional flavor of CTMQC, which conserves energy across the ensemble of trajectories (known as CTMQC-E), fails to reproduce detailed balance as in Ehrenfest dynamics, we show that a similar variant, where total energy is conserved on each trajectory independently, provides a major improvement over Ehrenfest with respect to detailed balance. Moreover, we show that the theory achieves convergence of the mean electronic potential energy with the number of energy levels that successively increase in energy. This new variant is shown to, by simulations on the Tully models and double arch model, retain a good description electronic populations and coherence compared to exact quantum dynamics. We explain the thermalization mechanism through the additional terms that distinguish CTMQC from Ehrenfest dynamics. We show that the improvement can be explained via geometric contributions to the nuclear force, resulting from the quantum momentum, which act to oppose motion when electrons decohere upward in energy and act to enhance motion otherwise, somewhat emulating the mechanism of frustrated hops. These results have considerable implications for the applicability of CTMQC to condensed phase simulations.

Direct intermolecular interactions, being short-range and, hence, strongly dependent on local disorder, are commonly neglected in simulations of organic exciton-polaritons, and molecules are only indirectly linked through the cavity electromagnetic field, irrespective of intermolecular separations. Whereas accounting for direct intermolecular interactions has presumably little effect on the properties of organic exciton-polaritons in many experimental systems, disregarding these interactions might no longer be valid in systems with a certain degree of structural ordering, such as organic crystals, conjugated polymers, or molecular aggregates. In these systems, intermolecular couplings are comparable to the experimentally achieved collective light-matter coupling strengths and, therefore, may modify the properties of polaritons compared to those in weakly interacting molecular systems. To test this, we incorporate nearest-neighbor excitonic couplings into the multi-mode Tavis-Cummings Hamiltonian and perform numerical simulations of polariton delocalization and transport. The simulation results suggest that negative, or J-aggregate-type, coupling alters the energetics of lower polariton states such that these states become more robust against static excitation energy disorder. This results in better transport properties, in particular in less velocity renormalization and the transition from ballistic transport to diffusion occurring at larger exciton fractions than in the absence of excitonic coupling. In contrast, positive, or H-aggregate-type, coupling makes the lower polariton states more sensitive to disorder and deteriorates their propagation. These findings emphasize the importance of intermolecular interactions as a control parameter for achieving long-range excitation energy transfer in optical microcavities.

Nonadiabatic molecular dynamics simulation of charge and exciton transport in molecular materials and biological systems are often carried out in a (quasi-)diabatic or site basis. Such simulations require the calculation of the electrostatic site energy of all possible charge or excited states of the system at each molecular dynamics step, which quickly becomes computationally prohibitive when Ewald summation is used. By combining the damped shifted force real space electrostatic summation method with a suitable addition-subtraction scheme, we show that the calculation of electrostatic energy and forces for Nmol site energies can be carried out at a small and system size independent overhead compared to the calculation for a single site energy. This advance enables us to include full electrostatic interactions in nonadiabatic molecular dynamics simulations for charge and exciton transport. Applying our computational scheme to hole transport in crystalline anthracene, we find that upon inclusion of electrostatic site energy fluctuations (also sometimes termed diagonal electrostatic disorder) the inverse participation ratio measuring hole delocalization decreases from ∼5 to ∼4 concomitant with a decrease in the hole mobility by about 9% along the b-crystallographic direction and by 30% along the a-direction. Accounting for electrostatics improves the agreement with experimental time-of-flight mobilities and mobility anisotropy, but it does not alter the charge transport mechanism, transient delocalization. Our work confirms that omission of electrostatic site energy disorder is a reasonable approximation for acenes, yet electrostatics is required to obtain near-quantitative agreement with experiment, even for apolar systems.

The field of organic photovoltaics has witnessed a renaissance in recent years owing to the development of non-fullerene acceptor materials reaching record power conversion efficiencies of > 20%. New computational models are needed to rationalise the regimes of photophysics reached in these materials. Here we report on a novel implementation of eXcitonic state-based Surface Hopping, a powerful non-adiabatic molecular dynamics method for the simulation of photo-induced charge generation in truly nanoscale donor-acceptor interfaces. We observe a transition from an inefficient cold to an efficient hot exciton dissociation mechanism as the electronic coupling between the molecules or the dielectric constant of the materials is increased. The hot pathway is observed to occur by Frenkel excitons converting into transiently delocalised hybrid exciton-charge transfer states that subsequently form charge separated states. This avoids the formation of kinetically trapped interfacial charge transfer states that are prone to non-radiative recombination.

Atomic-scale understanding of important geochemical processes including sorption, dissolution, nucleation, and crystal growth is difficult to obtain from experimental measurements alone and would benefit from strong continuous progress in molecular simulation. To this end, we present a reactive neural network potential-based molecular dynamics approach to simulate the interaction of aqueous ions on mineral surfaces in contact with liquid water, taking Fe(II) on hematite(001) as a model system. We show that a single neural network potential predicts rate constants for water exchange for aqueous Fe(II) and for the exergonic chemisorption of aqueous Fe(II) on hematite(001) in good agreement with experimental observations. The neural network potential developed herein allows one to converge free energy profiles and transmission coefficients at density functional theory-level accuracy outperforming state-of-the-art classical force field potentials. This suggests that machine learning potential molecular dynamics should become the method of choice for atomistic studies of geochemical processes.

Field-induced electron spin resonance provides valuable insights into the interplay between spin and charge dynamics in organic semiconductors. We apply this technique to ion-gel-gated capacitors and conventional field-effect transistors to study the temperature-dependent carrier dynamics of high-mobility rubrene single-crystals. Unlike previous measurements on other molecular and polymer semiconductors, we observe remarkably long spin relaxation times-on the order of microseconds-persisting from room temperature down to 15 K. Such long relaxation times are caused by the rapid transient-localization motion of charge carriers, which induces efficient motional narrowing. Additionally, by leveraging the high injection efficiency of ion-gel-gated devices, we observe spin lifetimes shortening at high carrier concentrations. This is attributed to emerging spin-spin dipolar interactions and can be modelled using an approach adapted from fluid-phase nuclear magnetic resonance. Our work demonstrates that field-induced electron spin resonance provides a powerful probe of the transient-localization physics of high-mobility molecular crystals.

The protonation state of molecules and surfaces is pivotal in various disciplines, including (electro-)catalysis, geochemistry, biochemistry, and pharmaceutics. Accurately and efficiently determining acidity constants is critical yet challenging, particularly when explicitly considering the electronic structure, thermal fluctuations, anharmonic vibrations, and solvation effects. In this research, we employ thermodynamic integration accelerated by committee Neural Network potentials, training a single machine learning model that accurately describes the relevant protonated, deprotonated, and intermediate states. We investigate two deprotonation reactions at the BiVO4 (010)-water interface, a promising candidate for efficient photocatalytic water splitting. Our results illustrate the convergence of the required ensemble averages over simulation time and of the final acidity constant as a function of the Kirkwood coupling parameter. We demonstrate that simulation times on the order of nanoseconds are required for statistical convergence. This time scale is currently unachievable with explicit ab-initio molecular dynamics simulations at the hybrid DFT level of theory. In contrast, our machine learning workflow only requires a few hundred DFT single point calculations for training and testing. Exploiting the extended time scales accessible, we furthermore asses the effect of commonly applied bias potentials. Thus, our study significantly advances calculating free energy differences with ab-initio accuracy.

Abstract The interaction of condensed phase systems with external electric fields is of major importance in a myriad of processes in nature and technology, ranging from the field-directed motion of cells (galvanotaxis), to geochemistry and the formation of ice phases on planets, to field-directed chemical catalysis and energy storage and conversion systems including supercapacitors, batteries and solar cells. Molecular simulation in the presence of electric fields would give important atomistic insight into these processes but applications of the most accurate methods such as ab-initio molecular dynamics (AIMD) are limited in scope by their computational expense. Here we introduce Perturbed Neural Network Potential Molecular Dynamics (PNNP MD) to push back the accessible time and length scales of such simulations. We demonstrate that important dielectric properties of liquid water including the field-induced relaxation dynamics, the dielectric constant and the field-dependent IR spectrum can be machine learned up to surprisingly high field strengths of about 0.2 V Å−1 without loss in accuracy when compared to ab-initio molecular dynamics. This is remarkable because, in contrast to most previous approaches, the two neural networks on which PNNP MD is based are exclusively trained on molecular configurations sampled from zero-field MD simulations, demonstrating that the networks not only interpolate but also reliably extrapolate the field response. PNNP MD is based on rigorous theory yet it is simple, general, modular, and systematically improvable allowing us to obtain atomistic insight into the interaction of a wide range of condensed phase systems with external electric fields.

Multiheme cytochromes (MHCs) are bacterial electron-transfer proteins. We show from optical spectra and calculations that some of these cytochromes probably contain occupied and unoccupied bands formed from heme π and π* orbitals that span the protein. In the fully oxidised proteins, the unoccupied π*-bands are energetically above the redox-active frontier orbitals, which according to NMR data and calculations, are formed of Fe3+ t2g and porphyrin π-orbitals. These orbitals on different hemes are electronically coupled according to EPR data and calculations, but only weakly so. We suggest a role for the heme bands in the electronic conductivity of single MHCs in bioelectronic junctions that is distinct from the role of the redox-active Fe3+ t2g and porphyrin π-orbitals in physiological electron transfer.

Thermoelectric materials convert a temperature gradient into a voltage. This phenomenon is relatively well understood for inorganic materials but much less so for organic semiconductors (OSs). These materials present a challenge because the strong thermal fluctuations of electronic coupling between the molecules result in partially delocalized charge carriers that cannot be treated with traditional theories for thermoelectricity. Here, we develop a quantum dynamical simulation approach revealing in atomistic detail how the charge carrier wave function moves along a temperature gradient in an organic molecular crystal. We find that the wave function propagates from hot to cold in agreement with the experiment, and we obtain a Seebeck coefficient in good agreement with experimental measurements that are also reported in this work. Detailed analysis reveals that gradients in thermal electronic disorder play an important role in determining the magnitude of the Seebeck coefficient, opening unexplored avenues for the design of OSs with improved Seebeck coefficients.