Search

We describe how to simulate charge diffusion in organic semiconductors using a recently introduced mixed quantum-classical method, the mapping approach to surface hopping (MASH). In contrast to standard fewest-switches surface hopping, this method propagates the classical degrees of freedom deterministically on the most populated adiabatic electronic state. This correctly preserves the equilibrium distribution of a quantum charge coupled to classical phonons, allowing one to time-average along trajectories to improve the statistical convergence of the calculation. We illustrate the method with an application to a standard model for the charge transport in the direction of maximum mobility in crystalline rubrene. Because of its consistency with the equilibrium distribution, the present method gives a time-dependent diffusion coefficient that plateaus correctly to a long-time limiting value. The resulting mobility is somewhat higher than that of the relaxation time approximation, which uses a phenomenological relaxation parameter to obtain a non-zero diffusion coefficient from a calculation with static phonon disorder. However, it is very similar to the mobility obtained from Ehrenfest dynamics, at least in the parameter regimes we have investigated here. This is somewhat surprising because Ehrenfest dynamics overheats the electronic subsystem and is therefore inconsistent with the equilibrium distribution., Comment: 13 pages, 5 figures

We compare the recently introduced multi-state mapping approach to surface hopping (MASH) with the F\"orster and Redfield theories of excitation energy transfer. Whereas F\"orster theory relies on weak coupling between chromophores, and Redfield theory assumes the electronic excitations to be weakly coupled to fast chromophore vibrations, MASH is free from any perturbative or Markovian approximations. We illustrate this with an example application to the rate of energy transfer in a Frenkel-exciton dimer, showing that MASH interpolates correctly between the opposing regimes in which the F\"orster and Redfield results are reliable. We then compare the three methods for a realistic model of the Fenna-Matthews-Olson complex with a structured vibrational spectral density and static disorder in the excitation energies. In this case there are no exact results for comparison so we use MASH to assess the validity of F\"orster and Redfield theories. We find that F\"orster theory is the more accurate of the two on the picosecond timescale, as has been shown previously for a simpler model of this particular light-harvesting complex. We also explore various ways to sample the initial electronic state in MASH and find that they all give very similar results for exciton dynamics., Comment: 10 pages, 3 figures

We describe how the fast quasi-centroid molecular dynamics (f-QCMD) method can be applied to condensed phase systems by approximating the quasi-centroid potential of mean force as a sum of inter- and intra-molecular corrections to the classical interaction potential. The corrections are found by using a regularised iterative Boltzmann inversion procedure to recover the inter- and intra-molecular quasi-centroid distribution functions obtained from a path integral molecular dynamics simulation. The resulting methodology is found to give good agreement with a previously published QCMD dipole absorption spectrum for liquid water, and satisfactory agreement for ice. It also gives good agreement with spectra from a recent implementation of CMD that uses a pre-computed elevated temperature potential of mean force. Modern centroid molecular dynamics methods therefore appear to be reaching a consensus regarding the impact of nuclear quantum effects on the vibrational spectra of water and ice., Comment: 10 pages, 4 figures, plus supplementary information and toc figure

We describe a multiple electronic state adaptation of the mapping approach to surface hopping introduced recently by Mannouch and Richardson (J. Chem. Phys. 158, 104111 (2023)). This adaptation treats populations and coherences on an equal footing and is guaranteed to give populations in any electronic basis that tend to the correct quantum-classical equilibrium values in the long-time limit (assuming ergodicity). We demonstrate its accuracy by comparison with exact benchmark results for three- and seven-state models of the Fenna-Matthews-Olson complex, obtaining electronic populations and coherences that are significantly more accurate than those of fewest switches surface hopping and at least as good as those of any other semiclassical method we are aware of. Since these results were obtained by adapting the scheme of Mannouch and Richardson, we go on to compare our results with theirs for a variety of problems with two electronic states. We find that their method is sometimes more accurate, and especially so in the Marcus inverted regime. However, in other situations the accuracies are comparable, and since our scheme can be used with multiple electronic states it can be applied to a wider variety of systems., Comment: 15 pages, 9 figures

We assess the cavity molecular dynamics method for the calculation of vibrational polariton spectra, using liquid water as a specific example. We begin by disputing a recent suggestion that nuclear quantum effects may lead to a broadening of polariton bands, finding instead that they merely result in anharmonic red shifts in the polariton frequencies. We go on to show that our simulated cavity spectra can be reproduced to graphical accuracy with a harmonic model that uses just the cavity-free spectrum and the geometry of the cavity as input. We end by showing that this harmonic model can be combined with the experimental cavity-free spectrum to give results in good agreement with optical cavity measurements. Since the input to our harmonic model is equivalent to the input to the transfer matrix method of applied optics, we conclude that cavity molecular dynamics cannot provide any more insight into the effect of vibrational strong coupling on the absorption spectrum than this transfer matrix method, which is already widely used by experimentalists to corroborate their cavity results., Comment: 10 pages, 5 figures

Machine learning frameworks based on correlations of interatomic positions begin with a discretized description of the density of other atoms in the neighbourhood of each atom in the system. Symmetry considerations support the use of spherical harmonics to expand the angular dependence of this density, but there is as yet no clear rationale to choose one radial basis over another. Here we investigate the basis that results from the solution of the Laplacian eigenvalue problem within a sphere around the atom of interest. We show that this generates the smoothest possible basis of a given size within the sphere, and that a tensor product of Laplacian eigenstates also provides the smoothest possible basis for expanding any higher-order correlation of the atomic density within the appropriate hypersphere. We consider several unsupervised metrics of the quality of a basis for a given dataset, and show that the Laplacian eigenstate basis has a performance that is much better than some widely used basis sets and is competitive with data-driven bases that numerically optimize each metric. In supervised machine learning tests, we find that the optimal function smoothness of the Laplacian eigenstates leads to comparable or better performance than can be obtained from a data-driven basis of a similar size that has been optimized to describe the atom-density correlation for the specific dataset. We conclude that the smoothness of the basis functions is a key and hitherto largely overlooked aspect of successful atomic density representations.

We describe a fast implementation of the quasi-centroid molecular dynamics (QCMD) method in which the quasi-centroid potential of mean force is approximated as a separable correction to the classical interaction potential. This correction is obtained by first calculating quasi-centroid radial and angular distribution functions in a short path integral molecular dynamics simulation, and then using iterative Boltzmann inversion to obtain an effective classical potential that reproduces these distribution functions in a classical NVT simulation. We illustrate this approach with example applications to the vibrational spectra of gas phase molecules, obtaining excellent agreement with QCMD reference calculations for water and ammonia and good agreement with the quantum mechanical vibrational spectrum of methane., Comment: 7 pages, 3 figures, and supplementary material

We show that the centroid molecular dynamics (CMD) method provides a realistic way to calculate the thermal diffusivity $a=\lambda/\rho c_{\rm V}$ of a quantum mechanical liquid such as para-hydrogen. Once $a$ has been calculated, the thermal conductivity can be obtained from $\lambda=\rho c_{\rm V}a$, where $\rho$ is the density of the liquid and $c_{\rm V}$ is the constant-volume heat capacity. The use of this formula requires an accurate quantum mechanical heat capacity $c_{\rm V}$, which can be obtained from a path integral molecular dynamics simulation. The thermal diffusivity can be calculated either from the decay of the equilibrium density fluctuations in the liquid or by using the Green-Kubo relation to calculate the CMD approximation to $\lambda$ and then dividing this by the corresponding approximation to $\rho c_{\rm V}$. We show that both approaches give the same results for liquid para-hydrogen and that these results are in good agreement with experimental measurements of the thermal conductivity over a wide temperature range. In particular, they correctly predict a decrease in the thermal conductivity at low temperatures -- an effect that stems from the decrease in the quantum mechanical heat capacity and has eluded previous para-hydrogen simulations. We also show that the method gives equally good agreement with experimental measurements for the thermal conductivity of normal liquid helium., Comment: 13 pages, 9 figures, 3 tables

We show that the stochastic Schr\"odinger equation (SSE) provides an ideal way to simulate the quantum mechanical spin dynamics of radical pairs. Electron spin relaxation effects arising from fluctuations in the spin Hamiltonian are straightforward to include in this approach, and their treatment can be combined with a highly efficient stochastic evaluation of the trace over nuclear spin states that is required to compute experimental observables. These features are illustrated in example applications to a flavin-tryptophan radical pair of interest in avian magnetoreception, and to a problem involving spin-selective radical pair recombination along a molecular wire. In the first of these examples, the SSE is shown to be both more efficient and more widely applicable than a recent stochastic implementation of the Lindblad equation, which only provides a valid treatment of relaxation in the extreme-narrowing limit. In the second, the exact SSE results are used to assess the accuracy of a recently-proposed combination of Nakajima-Zwanzig theory for the spin relaxation and Schulten-Wolynes theory for the spin dynamics, which is applicable to radical pairs with many more nuclear spins. An appendix analyses the efficiency of trace sampling in some detail, highlighting the particular advantages of sampling with SU(N) coherent states., Comment: 14 pages, 4 figures

We present a simple method for the calculation of reaction rates in the Fermi golden-rule limit, which accurately captures the effects of tunnelling and zero-point energy. The method is based on a modification of the recently proposed golden-rule quantum transition state theory (GR-QTST) of Thapa, Fang and Richardson. While GR-QTST is not size consistent, leading to the possibility of unbounded errors in the rate, our modified method has no such issue and so can be reliably applied to condensed phase systems. Both methods involve path-integral sampling in a constrained ensemble; the two methods differ, however, in the choice of constraint functional. We demonstrate numerically that our modified method is as accurate as GR-QTST for the one-dimensional model considered by Thapa and coworkers. We then study a multi-dimensional spin-boson model, for which our method accurately predicts the true quantum rate, while GR-QTST breaks down with an increasing number of boson modes in the discretisation of the spectral density. Our method is able to accurately predict reaction rates in the Marcus inverted regime, without the need for the analytic continuation required by Wolynes theory., Comment: 13 pages, 4 figures

We revisit the well-known aqueous ferrous-ferric electron transfer reaction in order to address recent suggestions that nuclear tunnelling can lead to significant deviation from the linear response assumption inherent in the Marcus picture of electron transfer. A recent study of this reaction by Richardson and coworkers has found a large difference between their new path-integral method, GR-QTST, and the saddle point approximation of Wolynes (Wolynes theory). They suggested that this difference could be attributed to the existence of multiple tunnelling pathways, leading Wolynes theory to significantly overestimate the rate. This was used to argue that the linear response assumptions of Marcus theory may break down for liquid systems when tunnelling is important. If true, this would imply that the commonly used method for studying such systems, where the problem is mapped onto a spin-boson model, is invalid. However, we have recently shown that size inconsistency in GR-QTST can lead to poor predictions of the rate in systems with many degrees of freedom. We have also suggested an improved method, the path-integral linear golden-rule (LGR) approximation, which fixes this problem. Here we demonstrate that the GR-QTST results for ferrous-ferric electron transfer are indeed dominated by its size consistency error. Furthermore, by comparing the LGR and Wolynes theory results, we confirm the established picture of nuclear tunnelling in this system. Finally, by comparing our path-integral results to those obtained by mapping onto the spin-boson model, we reassess the importance of anharmonic effects and the accuracy of this commonly used mapping approach., Comment: 15 pages, 5 figures

We present a general quantum instanton approach to calculating reaction rates for systems with two electronic states and arbitrary values of the electronic coupling. This new approach, which we call the non-adiabatic quantum instanton (NAQI) approximation, reduces to Wolynes theory in the golden rule limit and to a recently proposed projected quantum instanton (PQI) method in the adiabatic limit. As in both of these earlier theories, the NAQI approach is based on making a saddle point approximation to the time integral of a reactive flux autocorrelation function, although with a generalised definition of the projection operator onto the product states. We illustrate the accuracy of the approach by comparison with exact rates for one dimensional scattering problems and discuss its applicability to more complex reactions., Comment: 9 pages, 6 figures

Radical pair recombination reactions are known to be sensitive to extremely weak magnetic fields, and can therefore be said to function as molecular magnetoreceptors. The classic example is a carotenoid-porphyrin-fullerene (C+PF-) radical pair that has been shown to provide a "proof-of-principle" for the operation of a chemical compass [K. Maeda et al., Nature 453, 387 (2008)]. Previous simulations of this radical pair have employed semiclassical approximations, which are routinely applicable to its 47 coupled electronic and nuclear spins. However, calculating the exact quantum mechanical spin dynamics presents a significant challenge, and has not been possible before now. Here we use a recently developed method to perform numerically converged simulations of the C+PF- quantum mechanical spin dynamics, including all coupled spins. Comparison of these quantum mechanical simulations with various semiclassical approximations reveals that, while it is not perfect, the best semiclassical approximation does capture essentially all of the relevant physics in this problem., Comment: 8 pages, 4 figures

We analyse the golden rule limit of the recently proposed isomorphic ring polymer (iso-RP) method. This method aims to combine an exact expression for the quantum mechanical partition function of a system with multiple electronic states with a pre-existing mixed quantum-classical (MQC) dynamics approximation, such as fewest switches surface hopping. Since the choice of the MQC method adds a degree of flexibility, we simplify the analysis by assuming that the dynamics used correctly reproduces the exact golden rule rate for a non-adiabatic (e.g., electron transfer) reaction in the high temperature limit. Having made this assumption, we obtain an expression for the iso-RP rate in the golden rule limit that is valid at any temperature. We then compare this rate with the exact rate for a series of simple spin-boson models. We find that the iso-RP method does not correctly predict how nuclear quantum effects affect the reaction rate in the golden rule limit. Most notably, it does not capture the quantum asymmetry in a conventional (Marcus) plot of the logarithm of the reaction rate against the thermodynamic driving force, and it also significantly overestimates the correct quantum mechanical golden rule rate for activationless electron transfer reactions. These results are analysed and their implications discussed for the applicability of the iso-RP method to more general non-adiabatic reactions., Comment: 8 pages, 2 figures

Relaxation processes can have a large effect on the spin selective electron transfer reactions of radical pairs. These processes are often treated using phenomenological relaxation superoperators or with some model for the microscopic relaxation mechanism treated within Bloch-Redfield-Wangsness theory. Here, we demonstrate that an alternative perturbative relaxation theory, based on the Nakajima-Zwanzig equation, has certain advantages over Redfield theory. In particular, the Nakajima-Zwanzig equation does not suffer from the severe positivity problem of Redfield theory in the static disorder limit. Combining the Nakajima-Zwanzig approach consistently with the Schulten-Wolynes semiclassical method, we obtain an efficient method for modeling the spin dynamics of radical pairs containing many hyperfine-coupled nuclear spins. This is then used to investigate the spin-dependent electron transfer reactions and intersystem crossing of dimethyljulolidine-naphthalenediimide (DMJ-NDI) radical ion pairs. By comparing our simulations with experimental data, we find evidence for a field-independent contribution to the triplet quantum yields of these reactions which cannot be explained by electron spin relaxation alone.

We present a simple interpolation formula for the rate of an electron transfer reaction as a function of the electronic coupling strength. The formula only requires the calculation of Fermi Golden Rule and Born-Oppenheimer rates and so can be combined with any methods that are able to calculate these rates. We first demonstrate the accuracy of the formula by applying it to a one dimensional scattering problem for which the exact quantum mechanical, Fermi Golden Rule, and Born-Oppenheimer rates are readily calculated. We then describe how the formula can be combined with the Wolynes theory approximation to the Golden Rule rate, and the ring polymer molecular dynamics (RPMD) approximation to the Born-Oppenheimer rate, and used to capture the effects of nuclear tunnelling, zero point energy, and solvent friction on condensed phase electron transfer reactions. Comparison with exact hierarchical equations of motion (HEOM) results for a demanding set of spin-boson models shows that the interpolation formula has an error comparable to that of RPMD rate theory in the adiabatic limit, and that of Wolynes theory in non-adiabatic limit, and is therefore as accurate as any method could possibly be that attempts to generalise these methods to arbitrary electronic coupling strengths., Comment: 17 pages, 5 figures

Magnetic field effects on radical pair reactions arise due to the interplay of coherent electron spin dynamics and spin relaxation effects, a rigorous treatment of which requires the solution of the Liouville-von Neumann equation. However, it is often found that simple incoherent kinetic models of the radical pair singlet-triplet intersystem crossing provide an acceptable description of experimental measurements. In this paper we outline the theoretical basis for this incoherent kinetic description, elucidating its connection to exact quantum mechanics. We show in particular how the finite lifetime of the radical pair spin states, as well as any additional spin-state dephasing, leads to incoherent intersystem crossing. We arrive at simple expressions for the radical pair spin state interconversion rates to which the functional form proposed recently by Steiner et al. [J. Phys. Chem. C 122, 11701 (2018)] can be regarded as an approximation. We also test the kinetic master equation against exact quantum dynamical simulations for a model radical pair and for a series of $\text{PTZ}^{\bullet+}\text{-Ph}_\text{n}\text{-PDI}^{\bullet-}$ molecular wires.

Progress in the atomic-scale modelling of matter over the past decade has been tremendous. This progress has been brought about by improvements in methods for evaluating interatomic forces that work by either solving the electronic structure problem explicitly, or by computing accurate approximations of the solution and by the development of techniques that use the Born-Oppenheimer (BO) forces to move the atoms on the BO potential energy surface. As a consequence of these developments it is now possible to identify stable or metastable states, to sample configurations consistent with the appropriate thermodynamic ensemble, and to estimate the kinetics of reactions and phase transitions. All too often, however, progress is slowed down by the bottleneck associated with implementing new optimization algorithms and/or sampling techniques into the many existing electronic-structure and empirical-potential codes. To address this problem, we are thus releasing a new version of the i-PI software. This piece of software is an easily extensible framework for implementing advanced atomistic simulation techniques using interatomic potentials and forces calculated by an external driver code. While the original version of the code was developed with a focus on path integral molecular dynamics techniques, this second release of i-PI not only includes several new advanced path integral methods, but also offers other classes of algorithms. In other words, i-PI is moving towards becoming a universal force engine that is both modular and tightly coupled to the driver codes that evaluate the potential energy surface and its derivatives.

Radical pair recombination reactions are normally described using a quantum mechanical master equation for the electronic and nuclear spin density operator. The electron spin state selective (singlet and triplet) recombination processes are described with a Haberkorn reaction term in this master equation. Here we consider a general spin state selective electron transfer reaction of a radical pair and use Nakajima-Zwanzig theory to derive the master equation for the spin density operator, thereby elucidating the relationship between non-adiabatic reaction rate theory and the Haberkorn reaction term. A second order perturbation theory treatment of the diabatic coupling naturally results in the Haberkorn master equation with an additional reactive scalar electron spin coupling term. This term has been neglected in previous spin chemistry calculations, but we show that it will often be quite significant. We also show that beyond second order in perturbation theory, i.e., beyond the Fermi golden rule limit, an additional reactive singlet-triplet dephasing term appears in the master equation. A closed form expression for the reactive scalar electron spin coupling in terms of the Marcus theory parameters that determine the singlet and triplet recombination rates is presented. By performing simulations of radical pair reactions with the exact Hierarchical Equations of Motion (HEOM) method, we demonstrate that our master equations provide a very accurate description of radical pairs undergoing spin-selective non-adiabatic electron transfer reactions. The existence of a reactive electron spin coupling may well have implications for biologically relevant radical pair reactions such as those which have been suggested to play a role in avian magnetoreception.

Radical pair recombination reactions are known to be sensitive to the application of both low and high magnetic fields. The application of a weak magnetic field reduces the singlet yield of a singlet-born radical pair, whereas the application of a strong magnetic field increases the singlet yield. The high field effect arises from energy conservation: when the magnetic field is stronger than the sum of the hyperfine fields in the two radicals, ${\rm S}\to {\rm T}_{\pm}$ transitions become energetically forbidden, thereby reducing the number of pathways for singlet to triplet interconversion. The low field effect arises from symmetry breaking: the application of a weak magnetic field lifts degeneracies among the zero field eigenstates and increases the number of pathways for singlet to triplet interconversion. However, the details of this effect are more subtle, and have not previously been properly explained. Here we present a complete analysis of the low field effect in a radical pair containing a single proton, and in a radical pair in which one of the radicals contains a large number of hyperfine-coupled nuclear spins. We find that the new transitions that occur when the field is switched on are between ${\rm S}$ and ${\rm T}_0$ in both cases, and not between ${\rm S}$ and ${\rm T}_{\pm}$ as has previously been claimed. We then illustrate this result by using it in conjunction with semiclassical spin dynamics simulations to account for the observation of a biphasic--triphasic--biphasic transition with increasing magnetic field strength in the magnetic field effect on the time-dependent survival probability of a photoexcited carotenoid-porphyrin-fullerene radical pair., Comment: Article: 8 pages and 4 figures. Supplementary material: 6 pages and 1 figure

We describe a simple quantum mechanical method that can be used to obtain accurate numerical results over long time scales for the spin correlation tensor of an electron spin that is hyperfine coupled to a large number of nuclear spins. This method does not suffer from the statistical errors that accompany a Monte Carlo sampling of the exact eigenstates of the central spin Hamiltonian obtained from the algebraic Bethe ansatz, or from the growth of the truncation error with time in the time-dependent density matrix renormalization group (t-DMRG) approach. As a result, it can be applied to larger central spin problems than the algebraic Bethe ansatz, and for longer times than the t-DMRG algorithm. It is therefore an ideal method to use to solve central spin problems, and we expect that it will also prove useful for a variety of related problems that arise in a number of different research fields., Comment: Article: 5 pages, 5 figures Supplementary material: 15 pages, 6 figures

The Wolynes theory of electronically nonadiabatic reaction rates [P. G. Wolynes, J. Chem. Phys. 87, 6559 (1987)] is based on a saddle point approximation to the time integral of a reactive flux autocorrelation function in the nonadiabatic (golden rule) limit. The dominant saddle point is on the imaginary time axis at $t_{\rm sp}=i\lambda_{\rm sp}\hbar$, and provided $\lambda_{\rm sp}$ lies in the range $-\beta/2\le\lambda_{\rm sp}\le\beta/2$, it is straightforward to evaluate the rate constant using information obtained from an imaginary time path integral calculation. However, if $\lambda_{\rm sp}$ lies outside this range, as it does in the Marcus inverted regime, the path integral diverges. This has led to claims in the literature that Wolynes theory cannot describe the correct behaviour in the inverted regime. Here we show how the imaginary time correlation function obtained from a path integral calculation can be analytically continued to $\lambda_{\rm sp}<-\beta/2$, and the continuation used to evaluate the rate in the inverted regime. Comparisons with exact golden rule results for a spin-boson model and a more demanding (asymmetric and anharmonic) model of electronic predissociation show that the theory it is just as accurate in the inverted regime as it is in the normal regime., Comment: 9 pages, 8 figures

We have used an efficient new quantum mechanical method for radical pair recombination reactions to study the spin-dependent charge recombination along PTZ$^{\bullet+}$--Ph$_n$--PDI$^{\bullet-}$ molecular wires. By comparing our results to the experimental data of E. Weiss {\em et al.} [J. Am. Chem. Soc. {\bf 126}, 5577 (2004)], we are able to extract the spin-dependent (singlet and triplet) charge recombination rate constants for wires with $n=2-5$. These spin-dependent rate constants have not been extracted previously from the experimental data because they require fitting its magnetic field-dependence to the results of quantum spin dynamics simulations. We find that the triplet recombination rate constant decreases exponentially with the length of the wire, consistent with the superexchange mechanism of charge recombination. However, the singlet recombination rate constant is nearly independent of the length of the wire, suggesting that the singlet pathway is dominated by an incoherent hopping mechanism. A simple qualitative explanation for the different behaviours of the two spin-selective charge recombination pathways is provided in terms of Marcus theory. We also find evidence for a magnetic field-independent background contribution to the triplet yield of the charge recombination reaction, and suggest several possible explanations for it. Since none of these explanations is especially compelling given the available experimental evidence, and since the result appears to apply more generally to other molecular wires, we hope that this aspect of our study will stimulate further experimental work., Comment: 12 pages, 10 figures

Long-range ion induced water-water correlations were recently observed in femtosecond elastic second harmonic scattering experiments of electrolyte solutions. To further the qualitative understanding of these correlations, we derive an analytical expression that quantifies ion induced dipole-dipole correlations in a non-interacting gas of dipoles. This model is a logical extension of Debye-H\"uckel theory that can be used to qualitatively understand how the combined electric field of the ions induces correlations in the orientational distributions of the water molecules in an aqueous solution. The model agrees with results from molecular dynamics simulations and provides an important starting point for further theoretical work.

The standard quantum mechanical expressions for the singlet and triplet survival probabilities and product yields of a radical pair recombination reaction involve a trace over the states in a combined electronic and nuclear spin Hilbert space. If this trace is evaluated deterministically, by performing a separate time-dependent wavepacket calculation for each initial state in the Hilbert space, the computational effort scales as $O(Z^2\log Z)$, where $Z$ is the total number of nuclear spin states. Here we show that the trace can also be evaluated stochastically, by exploiting the properties of spin coherent states. This results in a computational effort of $O(MZ\log Z)$, where $M$ is the number of Monte Carlo samples needed for convergence. Example calculations on a strongly-coupled radical pair with $Z>10^6$ show that the singlet yield can be converged to graphical accuracy using just $M=200$ samples, resulting in a speed up by a factor of $>5000$ over a standard deterministic calculation. We expect that this factor will greatly facilitate future quantum mechanical simulations of a wide variety of radical pairs of interest in chemistry and biology., Comment: 6 pages, 5 figures

We present a new method for calculating the product yield of a radical pair recombination reaction in the presence of a weak time-dependent magnetic field. This method successfully circumvents the computational difficulties presented by a direct solution of the Liouville-von Neumann equation for a long-lived radical pair containing many hyperfine-coupled nuclear spins. Using a modified formulation of Floquet theory, treating the time-dependent magnetic field as a perturbation, and exploiting the slow radical pair recombination, we show that one can obtain a good approximation to the product yield by considering only nearly-degenerate sub-spaces of the Floquet space. Within a significant parameter range, the resulting method is found to give product yields in good agreement with exact quantum mechanical results for a variety of simple model radical pairs. Moreover it is considerably more efficient than the exact calculation, and it can be applied to radical pairs containing significantly more nuclear spins. This promises to open the door to realistic theoretical investigations of the effect of radiofrequency electromagnetic radiation on the photochemically induced radical pair recombination reactions in the avian retina which are believed to be responsible for the magnetic compass sense of migratory birds., Comment: 11 pages, 7 figures

The magnetoelectroluminescence of conjugated organic polymer films is widely accepted to arise from a polaron pair mechanism, but their magnetoconductance is less well understood. Here we derive a new relationship between the experimentally measurable magnetoelectroluminescence and magnetoconductance and the theoretically calculable singlet yield of the polaron pair recombination reaction. This relationship is expected to be valid regardless of the mechanism of the magnetoconductance, provided the mobilities of the free polarons are independent of the applied magnetic field (i.e., provided one discounts the possibility of spin-dependent transport). We also discuss the semiclassical calculation of the singlet yield of the polaron pair recombination reaction for materials such as poly(2,5-dioctyloxy-paraphenylene vinylene) (DOO-PPV), the hyperfine fields in the polarons of which can be extracted from light-induced electron spin resonance measurements. The resulting theory is shown to give good agreement with experimental data for both normal (H-) and deuterated (D-) DOO-PPV over a wide range of magnetic field strengths once singlet-triplet dephasing is taken into account. Without this effect, which has not been included in any previous simulation of magnetoelectroluminescence, it is not possible to reproduce the experimental data for both isotopologues in a consistent fashion. Our results also indicate that the magnetoconductance of DOO-PPV cannot be solely due to the effect of the magnetic field on the dissociation of polaron pairs.

We employ classical and ring polymer molecular dynamics simulations to study the effect of nuclear quantum fluctuations on the structure and the water exchange dynamics of aqueous solutions of lithium and fluoride ions. While we obtain reasonably good agreement with experimental data for solutions of lithium by augmenting the Coulombic interactions between the ion and the water molecules with a standard Lennard-Jones ion-oxygen potential, the same is not true for solutions of fluoride, for which we find that a potential with a softer repulsive wall gives much better agreement. A small degree of destabilization of the first hydration shell is found in quantum simulations of both ions when compared with classical simulations, with the shell becoming less sharply defined and the mean residence time of the water molecules in the shell decreasing. In line with these modest differences, we find that the mechanisms of the exchange processes are unaffected by quantization, so a classical description of these reactions gives qualitatively correct and quantitatively reasonable results. We also find that the quantum effects in solutions of lithium are larger than in solutions of fluoride. This is partly due to the stronger interaction of lithium with water molecules, partly due to the lighter mass of lithium, and partly due to competing quantum effects in the hydration of fluoride, which are absent in the hydration of lithium., Comment: 12 pages, 8 figures

Light nuclei at room temperature and below exhibit a kinetic energy which significantly deviates from the predictions of classical statistical mechanics. This quantum kinetic energy is responsible for a wide variety of isotope effects of interest in fields ranging from chemistry to climatology. It also furnishes the second moment of the nuclear momentum distribution, which contains subtle information about the chemical environment and has recently become accessible to deep inelastic neutron scattering experiments. Here we show how, by combining imaginary time path integral dynamics with a carefully designed generalized Langevin equation, it is possible to dramatically reduce the expense of computing the quantum kinetic energy. We also introduce a transient anisotropic Gaussian approximation to the nuclear momentum distribution which can be calculated with negligible additional effort. As an example, we evaluate the structural properties, the quantum kinetic energy, and the nuclear momentum distribution for a first-principles simulation of liquid water.

We present a practical semiclassical method for computing the electron spin dynamics of a radical in which the electron spin is hyperfine coupled to a large number of nuclear spins. This can be used to calculate the singlet and triplet survival probabilities and quantum yields of radical recombination reactions in the presence of magnetic fields. Our method differs from the early semiclassical theory of Schulten and Wolynes [J. Chem. Phys. 68, 3292 (1978)] in allowing each individual nuclear spin to precess around the electron spin, rather than assuming that the hyperfine coupling-weighted sum of nuclear spin vectors is fixed in space. The downside of removing this assumption is that one can no longer obtain a simple closed-form expression for the electron spin correlation tensor: our method requires a numerical calculation. However, the computational effort increases only linearly with the number of nuclear spins, rather than exponentially as in an exact quantum mechanical calculation. The method is therefore applicable to arbitrarily large radicals. Moreover, it approaches quantitative agreement with quantum mechanics as the number of nuclear spins increases and the environment of the electron spin becomes more complex, owing to the rapid quantum decoherence in complex systems. Unlike the Schulten-Wolynes theory, the present semiclassical theory predicts the correct long-time behaviour of the electron spin correlation tensor, and it therefore correctly captures the low magnetic field effect in the singlet yield of a radical recombination reaction with a slow recombination rate.

We describe how the semiclassical theory of radical pair recombination reactions recently introduced by two of us [D. E. Manolopoulos and P. J. Hore, J. Chem. Phys. 139, 124106 (2013)] can be generalised to allow for different singlet and triplet recombination rates. This is a non-trivial generalisation because when the recombination rates are different the recombination process is dynamically coupled to the coherent electron spin dynamics of the radical pair. Furthermore, because the recombination operator is a two-electron operator, it is no longer sufficient simply to consider the two electrons as classical vectors: one has to consider the complete set of 16 two-electron spin operators as independent classical variables. The resulting semiclassical theory is first validated by comparison with exact quantum mechanical results for a model radical pair containing 12 nuclear spins. It is then used to shed light on the spin dynamics of a carotenoid-porphyrin-fullerene (CPF) triad containing considerably more nuclear spins which has recently been used to establish a 'proof of principle' for the operation of a chemical compass [K. Maeda et al., Nature 453, 387 (2008)]. We find in particular that the intriguing biphasic behaviour that has been observed in the effect of an Earth-strength magnetic field on the time-dependent survival probability of the photo-excited C+PF- radical pair arises from a delicate balance between its asymmetric recombination and the relaxation of the electron spin in the carotenoid radical.

Two of the most successful methods that are presently available for simulating the quantum dynamics of condensed phase systems are centroid molecular dynamics (CMD) and ring polymer molecular dynamics (RPMD). Despite their conceptual differences, practical implementations of these methods differ in just two respects: the choice of the Parrinello-Rahman mass matrix and whether or not a thermostat is applied to the internal modes of the ring polymer during the dynamics. Here we explore a method which is halfway between the two approximations: we keep the path integral bead masses equal to the physical particle masses but attach a Langevin thermostat to the internal modes of the ring polymer during the dynamics. We justify this by showing analytically that the inclusion of an internal mode thermostat does not affect any of the desirable features of RPMD: thermostatted RPMD (TRPMD) is equally valid with respect to everything that has actually been proven about the method as RPMD itself. In particular, because of the choice of bead masses, the resulting method is still optimum in the short-time limit, and the transition state approximation to its reaction rate theory remains closely related to the semiclassical instanton approximation in the deep quantum tunneling regime. In effect, there is a continuous family of methods with these properties, parameterised by the strength of the Langevin friction. Here we explore numerically how the approximation to quantum dynamics depends on this friction, with a particular emphasis on vibrational spectroscopy. We find that a broad range of frictions approaching optimal damping give similar results, and that these results are immune to both the resonance problem of RPMD and the curvature problem of CMD.

Recent developments in path integral methodology have significantly reduced the computational expense of including quantum mechanical effects in the nuclear motion in ab initio molecular dynamics simulations. However, the implementation of these developments requires a considerable programming effort, which has hindered their adoption. Here we describe i-PI, an interface written in Python that has been designed to minimise the effort required to bring state-of-the-art path integral techniques to an electronic structure program. While it is best suited to first principles calculations and path integral molecular dynamics, i-PI can also be used to perform classical molecular dynamics simulations, and can just as easily be interfaced with an empirical forcefield code. To give just one example of the many potential applications of the interface, we use it in conjunction with the CP2K electronic structure package to showcase the importance of nuclear quantum effects in high pressure water.

Even at room temperature, quantum mechanics plays a major role in determining the quantitative behaviour of light nuclei, changing significantly the values of physical properties such as the heat capacity. However, other observables appear to be only weakly affected by nuclear quantum effects (NQEs): for instance, the melting temperatures of light and heavy water differ by less than 4 K. Recent theoretical work has attributed this to a competition between intra and inter molecular NQEs, which can be separated by computing the anisotropy of the quantum kinetic energy tensor. The principal values of this tensor change in opposite directions when ice melts, leading to a very small net quantum mechanical effect on the melting point. This paper presents the first direct experimental observation of this phenomenon, achieved by measuring the deuterium momentum distributions n(p) in heavy water and ice using Deep Inelastic Neutron Scattering (DINS), and resolving their anisotropy. Results from the experiments, supplemented by a theoretical analysis, show that the anisotropy of the quantum kinetic energy tensor can also be captured for heavier atoms such as oxygen.

The quantum nature of nuclei plays an important role in the accurate modelling of light atoms such as hydrogen, but it is often neglected in simulations due to the high computational overhead involved. It has recently been shown that zero-point energy effects can be included comparatively cheaply in simulations of harmonic and quasi-harmonic systems by augmenting classical molecular dynamics with a generalized Langevin equation (GLE). Here we describe how a similar approach can be used to accelerate the convergence of path integral (PI) molecular dynamics to the exact quantum mechanical result in more strongly anharmonic systems exhibiting both zero point energy and tunnelling effects. The resulting PI-GLE method is illustrated with applications to a double-well tunnelling problem and to liquid water.

Computing averages over a target probability density by statistical re-weighting of a set of samples with a different distribution is a strategy which is commonly adopted in fields as diverse as atomistic simulation and finance. Here we present a very general analysis of the accuracy and efficiency of this approach, highlighting some of its weaknesses. We then give an example of how our results can be used, specifically to assess the feasibility of high-order path integral methods. We demonstrate that the most promising of these techniques -- which is based on re-weighted sampling -- is bound to fail as the size of the system is increased, because of the exponential growth of the statistical uncertainty in the re-weighted average.

Numerous studies have identified large quantum mechanical effects in the dynamics of liquid water. In this paper, we suggest that these effects may have been overestimated due to the use of rigid water models and flexible models in which the intramolecular interactions were described using simple harmonic functions. To demonstrate this, we introduce a new simple point charge model for liquid water, q-TIP4P/F, in which the O--H stretches are described by Morse-type functions. We have parameterized this model to give the correct liquid structure, diffusion coefficient, and infra-red absorption frequencies in quantum (path integral-based) simulations. By comparing classical and quantum simulations of the liquid, we find that quantum mechanical fluctuations increase the rates of translational diffusion and orientational relaxation in our model by a factor of around 1.15. This effect is much smaller than that observed in all previous simulations of simple empirical water models, which have found a quantum effect of at least 1.4 regardless of the quantum simulation method or the water model employed. The small quantum effect in our model is a result of two competing phenomena. Intermolecular zero point energy and tunneling effects destabilize the hydrogen bonding network, leading to a less viscous liquid with a larger diffusion coefficient. However this is offset by intramolecular zero point motion, which changes the average water monomer geometry resulting in a larger dipole moment, stronger intermolecular interactions, and slower diffusion. We end by suggesting, on the basis of simulations of other potential energy models, that the small quantum effect we find in the diffusion coefficient is associated with the ability of our model to produce a single broad O-H stretching band in the infra-red absorption spectrum., Comment: Published in JCP

The path integral molecular dynamics (PIMD) method provides a convenient way to compute the quantum mechanical structural and thermodynamic properties of condensed phase systems at the expense of introducing an additional set of high-frequency normal modes on top of the physical vibrations of the system. Efficiently sampling such a wide range of frequencies provides a considerable thermostatting challenge. Here we introduce a simple stochastic path integral Langevin equation (PILE) thermostat which exploits an analytic knowledge of the free path integral normal mode frequencies. We also apply a recently-developed colored-noise thermostat based on a generalized Langevin equation (GLE), which automatically achieves a similar, frequency-optimized sampling. The sampling efficiencies of these thermostats are compared with that of the more conventional Nos\'e-Hoover chain (NHC) thermostat for a number of physically relevant properties of the liquid water and hydrogen-in-palladium systems. In nearly every case, the new PILE thermostat is found to perform just as well as the NHC thermostat while allowing for a computationally more efficient implementation. The GLE thermostat also proves to be very robust delivering a near-optimum sampling efficiency in all of the cases considered. We suspect that these simple stochastic thermostats will therefore find useful application in many future PIMD simulations., Comment: Accepted for publication on JCP