Volume 133, Issue 24, 28 December 2010

We present a plane wave basis set implementation for the calculation of electronic coupling matrix elements of electron transferreactions within the framework of constrained density functional theory (CDFT). Following the work of Wu and Van Voorhis [J. Chem. Phys. 125, 164105 (2006)], the diabatic wavefunctions are approximated by the Kohn–Sham determinants obtained from CDFT calculations, and the coupling matrix element calculated by an efficient integration scheme. Our results for intermolecular electron transfer in small systems agree very well with highlevel ab initio calculations based on generalized Mulliken–Hush theory, and with previous local basis set CDFT calculations. The effect of thermal fluctuations on the coupling matrix element is demonstrated for intramolecular electron transfer in the tetrathiafulvalenediquinone (QTTFQ^{−}) anion. Sampling the electronic coupling along density functional based molecular dynamics trajectories, we find that thermal fluctuations, in particular the slow bending motion of the molecule, can lead to changes in the instantaneous electron transfer rate by more than an order of magnitude. The thermal average, is significantly higher than the value obtained for the minimum energy structure, . While CDFT in combination with generalized gradient approximation (GGA) functionals describes the intermolecular electron transfer in the studied systems well, exact exchange is required for QTTFQ^{−} in order to obtain coupling matrix elements in agreement with experiment (3.9 mH). The implementation presented opens up the possibility to compute electronic coupling matrix elements for extended systems where donor, acceptor, and the environment are treated at the quantum mechanical (QM) level.
 ARTICLES

 Theoretical Methods and Algorithms

Sampling rare events in nonequilibrium and nonstationary systems
View Description Hide DescriptionAlthough many computational methods for rare event sampling exist, this type of calculation is not usually practical for general nonequilibrium conditions, with macroscopically irreversible dynamics and away from both stationary and metastable states. A novel method for calculating the timeseries of the probability of a rare event is presented which is designed for these conditions. The method is validated for the cases of the Glauber–Ising model under timevarying shear flow, the Kawasaki–Ising model after a quench into the region between nucleation dominated and spinodal decomposition dominated phase change dynamics, and the parallel open asymmetric exclusion process. The method requires a subdivision of the phase space of the system: it is benchmarked and found to scale well for increasingly fine subdivisions, meaning that it can be applied without detailed foreknowledge of the physically important reaction pathways.

A simple protocol for the probability weights of the simulated tempering algorithm: Applications to firstorder phase transitions
View Description Hide DescriptionThe simulated tempering (ST) is an important method to deal with systems whose phase spaces are hard to sample ergodically. However, it uses accepting probabilities weights, which often demand involving and time consuming calculations. Here it is shown that such weights are quite accurately obtained from the largest eigenvalue of the transfer matrix—a quantity straightforward to compute from direct Monte Carlo simulations—thus simplifying the algorithm implementation. As tests, different systems are considered, namely, Ising, Blume–Capel, Blume–Emery–Griffiths, and Bell–Lavis liquid water models. In particular, we address firstorder phase transition at low temperatures, a regime notoriously difficulty to simulate because the large freeenergy barriers. The good results found (when compared with other well established approaches) suggest that the ST can be a valuable tool to address strong firstorder phase transitions, a possibility still not well explored in the literature.

Nonlocal van der Waals density functional: The simpler the better
View Description Hide DescriptionWe devise a nonlocal correlation energy functional that describes the entire range of dispersion interactions in a seamless fashion using only the electron density as input. The new functional is considerably simpler than its predecessors of a similar type. The functional has a tractable and robust analytic form that lends itself to efficient selfconsistent implementation. When paired with an appropriate exchange functional, our nonlocal correlation model yields accurate interaction energies of weaklybound complexes, not only near the energy minima but also far from equilibrium. Our model exhibits an outstanding precision at predicting equilibrium intermonomer separations in van der Waals complexes. It also gives accurate covalent bond lengths and atomization energies. Hence the functional proposed in this work is a computationally inexpensive electronic structure tool of broad applicability.

Effective local potentials for excited states
View Description Hide DescriptionThe constrained variational Hartree–Fock method for excited states of the same symmetry as the ground state [Chem. Phys. Lett. 287, 189 (1998)] is combined with the effective local potential (ELP) method [J. Chem. Phys. 125, 081104 (2006)] to generate Kohn–Shamtype exactexchange potentials for singly excited states of manyelectron systems. Illustrative examples include the three lowest states of the Li and Na atoms and the three lowest states of He and Be. For the systems studied, excitedstate ELPs differ from the corresponding groundstate potentials in two respects: They are less negative and have small additional “bumps” in the outer electron region. The technique is general and can be used to approximate excitedstate exchangecorrelation potentials for other orbitaldependent functionals.

Electronic coupling matrix elements from charge constrained density functional theory calculations using a plane wave basis set
View Description Hide DescriptionWe present a plane wave basis set implementation for the calculation of electronic coupling matrix elements of electron transferreactions within the framework of constrained density functional theory (CDFT). Following the work of Wu and Van Voorhis [J. Chem. Phys. 125, 164105 (2006)], the diabatic wavefunctions are approximated by the Kohn–Sham determinants obtained from CDFT calculations, and the coupling matrix element calculated by an efficient integration scheme. Our results for intermolecular electron transfer in small systems agree very well with highlevel ab initio calculations based on generalized Mulliken–Hush theory, and with previous local basis set CDFT calculations. The effect of thermal fluctuations on the coupling matrix element is demonstrated for intramolecular electron transfer in the tetrathiafulvalenediquinone (QTTFQ^{−}) anion. Sampling the electronic coupling along density functional based molecular dynamics trajectories, we find that thermal fluctuations, in particular the slow bending motion of the molecule, can lead to changes in the instantaneous electron transfer rate by more than an order of magnitude. The thermal average, is significantly higher than the value obtained for the minimum energy structure, . While CDFT in combination with generalized gradient approximation (GGA) functionals describes the intermolecular electron transfer in the studied systems well, exact exchange is required for QTTFQ^{−} in order to obtain coupling matrix elements in agreement with experiment (3.9 mH). The implementation presented opens up the possibility to compute electronic coupling matrix elements for extended systems where donor, acceptor, and the environment are treated at the quantum mechanical (QM) level.

Multidimensional optical spectroscopy of a single molecule in a currentcarrying state
View Description Hide DescriptionThe nonlinear optical signals from an open system consisting of a molecule connected to metallic leads, in response to a sequence of impulsive pulses, are calculated using a superoperator formalism. Two detection schemes are considered: coherent stimulated emission and incoherent fluorescence. The two provide similar but not identical information. The necessary superoperator correlation functions are evaluated either by converting them to ordinary (Hilbert space) operators which are then expanded in manybody states, or by using Wick's theorem for superoperators to factorize them into nonequilibrium two point Green's functions. As an example we discuss a stimulated Raman process that shows resonances involving two different charge states of the molecule in the same signal.

Densityfunctional expansion methods: Evaluation of LDA, GGA, and metaGGA functionals and different integral approximations
View Description Hide DescriptionWe extend the Kohn–Sham potential energy expansion (VE) to include variations of the kinetic energy density and use the VE formulation with a 631G* basis to perform a “Jacob's ladder” comparison of small molecule properties using density functionals classified as being either LDA, GGA, or metaGGA. We show that the VE reproduces standard Kohn–Sham DFT results well if all integrals are performed without further approximation, and there is no substantial improvement in using metaGGA functionals relative to GGA functionals. The advantages of using GGA versus LDA functionals becomes apparent when modeling hydrogen bonds. We furthermore examine the effect of using integral approximations to compute the zerothorder energy and firstorder matrix elements, and the results suggest that the origin of the shortrange repulsive potential within selfconsistent charge densityfunctional tightbinding methods mainly arises from the approximations made to the firstorder matrix elements.

Comparison of oneparticle basis set extrapolation to explicitly correlated methods for the calculation of accurate quartic force fields, vibrational frequencies, and spectroscopic constants: Application to H_{2}O, N_{2}H^{+}, NO_{2} ^{+}, and C_{2}H_{2}
View Description Hide DescriptionOneparticle basis set extrapolation is compared with one of the new R12 methods for computing highly accurate quartic force fields (QFFs) and spectroscopic data, including molecular structures, rotational constants, and vibrational frequencies for the H_{2}O, N_{2}H^{+}, NO_{2} ^{+}, and C_{2}H_{2} molecules. In general, agreement between the spectroscopic data computed from the best R12 and basis set extrapolation methods is very good with the exception of a few parameters for N_{2}H^{+} where it is concluded that basis set extrapolation is still preferred. The differences for H_{2}O and NO_{2} ^{+} are small and it is concluded that the QFFs from both approaches are more or less equivalent in accuracy. For C_{2}H_{2}, however, a known oneparticle basis set deficiency for C–C multiple bonds significantly degrades the quality of results obtained from basis set extrapolation and in this case the R12 approach is clearly preferred over oneparticle basis set extrapolation. The R12 approach used in the present study was modified in order to obtain high precision electronic energies, which are needed when computing a QFF. We also investigated including corecorrelation explicitly in the R12 calculations, but conclude that current approaches are lacking. Hence corecorrelation is computed as a correction using conventional methods. Considering the results for all four molecules, it is concluded that R12 methods will soon replace basis set extrapolation approaches for high accuracy electronic structure applications such as computing QFFs and spectroscopic data for comparison to highresolution laboratory or astronomical observations, provided one uses a robust R12 method as we have done here. The specific R12 method used in the present study, CCSD(T)_{R12}, incorporated a reformulation of one intermediate matrix in order to attain machine precision in the electronic energies. Final QFFs for N_{2}H^{+} and NO_{2} ^{+} were computed, including basis set extrapolation, corecorrelation, scalar relativity, and higherorder correlation and then used to compute highly accurate spectroscopic data for all isotopologues. Agreement with highresolution experiment for ^{14}N_{2}H^{+} and ^{14}N_{2}D^{+} was excellent, but for ^{14}N^{16}O_{2} ^{+} agreement for the two stretching fundamentals is outside the expected residual uncertainty in the theoretical values, and it is concluded that there is an error in the experimental quantities. It is hoped that the highly accurate spectroscopic data presented for the minor isotopologues of N_{2}H^{+} and NO_{2} ^{+} will be useful in the interpretation of future laboratory or astronomical observations.

Aggregation work at polydisperse micellization: Ideal solution and “dressed micelle” models comparing to molecular dynamics simulations
View Description Hide DescriptionGeneral thermodynamic relations for the work of polydisperse micelle formation in the model of ideal solution of molecular aggregates in nonionic surfactantsolution and the model of “dressed micelles” in ionic solution have been considered. In particular, the dependence of the aggregation work on the total concentration of nonionic surfactant has been analyzed. The analogous dependence for the work of formation of ionic aggregates has been examined with regard to existence of two variables of a state of an ionic aggregate, the aggregation numbers of surface active ions and counterions. To verify the thermodynamicmodels, the molecular dynamics simulations of micellization in nonionic and ionic surfactantsolutions at two total surfactant concentrations have been performed. It was shown that for nonionic surfactants, even at relatively high total surfactant concentrations, the shape and behavior of the work of polydisperse micelle formation found within the model of the ideal solution at different total surfactant concentrations agrees fairly well with the numerical experiment. For ionic surfactantsolutions, the numerical results indicate a strong screening of ionic aggregates by the bound counterions. This fact as well as independence of the coefficient in the law of mass action for ionic aggregates on total surfactant concentration and predictable behavior of the “waterfall” lines of surfaces of the aggregation work upholds the model of “dressed” ionic aggregates.

Local CC2 response method for triplet states based on Laplace transform: Excitation energies and firstorder properties
View Description Hide DescriptionA new multistate local CC2 response method for calculating excitation energies and firstorder properties of excited triplet states in extended molecular systems is presented. The Laplace transform technique is employed to partition the left/right local CC2 eigenvalue problems as well as the linear equations determining the Lagrange multipliers needed for the properties. The doubles part in the equations can then be inverted onthefly and only effective equations for the singles part must be solved iteratively. The local approximation presented here is adaptive and statespecific. The densityfitting method is utilized to approximate the electronrepulsion integrals. The accuracy of the new method is tested by comparison to canonical reference values for a set of 12 test molecules and 62 excited triplet states. As an illustrative application example, the lowest four triplet states of 3(5(5(4(bis(4(hexyloxy)phenyl)amino)phenyl)thiophene2yl)thiophene2yl)2cyanoacrylic acid, an organic sensitizer for solarcell applications, are computed in the present work. No triplet chargetransfer states are detected among these states. This situation contrasts with the singlet states of this molecule, where the lowest singlet state has been recently found to correspond to an excited state with a pronounced chargetransfer character having a large transition strength.

A smooth, nonsingular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach
View Description Hide DescriptionPolarizable continuum models (PCMs) are a widely used family of implicit solvent models based on reactionfield theory and boundaryelement discretization of the solute/continuum interface. An often overlooked aspect of these theories is that discretization of the interface typically does not afford a continuous potential energy surface for the solute. In addition, we show that discretization can lead to numerical singularities and violations of exact variational conditions. To fix these problems, we introduce the switching/Gaussian (SWIG) method, a discretization scheme that overcomes several longstanding problems with PCMs. Our approach generalizes a procedure introduced by York and Karplus [J. Phys. Chem. A 103, 11060 (1999)], extending it beyond the conductorlike screening model. Comparison to other purportedly smooth PCM implementations reveals certain artifacts in these alternative approaches, which are avoided using the SWIG methodology. The versatility of our approach is demonstrated via geometry optimizations, vibrational frequency calculations, and molecular dynamics simulations, for solutes described using quantum mechanics and molecular mechanics.

An adaptive coupledcluster theory: @CC approach
View Description Hide DescriptionA formulation of an adaptive coupledcluster theory is presented. The method automatically “adjusts” to any state of an electronic system and converges to the full CI limit, thus being capable of describing both single and multireference phenomena. Adaptivity is accomplished through a guided selection of a compact set of cluster amplitudes as required for a proper description of the electronic system under consideration. The approach suggested is of “blackbox” type. A special importanceselection function (discriminatory function) is explicitly introduced for the guided selection of variables involved in the theoretical model. The method is tested on molecules which exhibit strong multireference character in the region of chemical bond elongation. An unambiguous comparison with formally exact full CI solutions shows that the method is capable of providing mHartee accuracy using a rather compact set of cluster amplitudes.

The fluxflux correlation function for anharmonic barriers
View Description Hide DescriptionThe fluxflux correlation function formalism is a standard and widely used approach for the computation of reaction rates. In this paper we introduce a method to compute the classical and quantum fluxflux correlation functions for anharmonic barriers essentially analytically through the use of the classical and quantum normal forms. In the quantum case we show that for a general f degreeoffreedom system having an index one saddle the quantum normal form reduces the computation of the fluxflux correlation function to that of an effective onedimensional anharmonic barrier. The example of the computation of the quantum fluxflux correlation function for a fourth order anharmonic barrier is worked out in detail, and we present an analytical expression for the quantum mechanical microcanonical fluxflux correlation function. We then give a discussion of the shorttime and harmonic limits.

Escape from the potential well: Competition between long jumps and long waiting times
View Description Hide DescriptionWithin a concept of the fractional diffusionequation and subordination, the paper examines the influence of a competition between long waiting times and long jumps on the escape from the potential well. Applying analytical arguments and numerical methods, we demonstrate that the presence of long waiting times distributed according to a powerlaw distribution with a diverging mean leads to very general asymptotic properties of the survival probability. The observed survival probability asymptotically decays like a power law whose form is not affected by the value of the exponent characterizing the power law jump length distribution. It is demonstrated that this behavior is typical of and generic for systems exhibiting long waiting times. We also show that the survival probability has a universal character not only asymptotically, but also at small times. Finally, it is indicated which properties of the first passage time density are sensitive to the exact value of the exponent characterizing the jump length distribution.

Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: A comparison of simulation techniques
View Description Hide DescriptionOver the last number of years several simulation methods have been introduced to study rare events such as nucleation. In this paper we examine the crystal nucleation rate of hard spheres using three such numerical techniques: molecular dynamics, forward flux sampling, and a Bennett–Chandlertype theory where the nucleation barrier is determined using umbrella sampling simulations. The resulting nucleation rates are compared with the experimental rates of Harland and van Megen [Phys. Rev. E 55, 3054 (1997)], Sinn et al. [Prog. Colloid Polym. Sci. 118, 266 (2001)], Schätzel and Ackerson [Phys. Rev. E 48, 3766 (1993)], and the predicted rates for monodisperse and 5% polydisperse hard spheres of Auer and Frenkel [Nature 409, 1020 (2001)]. When the rates are examined in units of the longtime diffusion coefficient, we find agreement between all the theoretically predicted nucleation rates, however, the experimental results display a markedly different behavior for low supersaturation. Additionally, we examined the precritical nuclei arising in the molecular dynamics, forward flux sampling, and umbrella sampling simulations. The structure of the nuclei appears independent of the simulation method, and in all cases, the nuclei contains on average significantly more facecenteredcubic ordered particles than hexagonalclosepacked ordered particles.

Understanding freeenergy perturbation calculations through a model of harmonic oscillators: Theory and implications to improve the sampling efficiency by molecular simulation
View Description Hide DescriptionFreeenergyperturbation calculation is frequently used to calculate freeenergy differences because it is easy to implement and the computation is fast. However, the calculation is subject to large inaccuracies in some circumstances due to the insufficient sampling of the relevant tails of the energydifference distributions. Here we expand this knowledge of insufficient sampling into a twodimensional (2D) energy space using a model of harmonic oscillators. We show analytically the relation between the energies of the sampling system and those of the desired target energy spaces, which provide the basis to understand the difficulties in freeenergyperturbation calculations. We clarify the reasons of the inaccurate calculation in the different harmonic cases that stem from the spatial separations of the reference and the target energy pairs located in the twodimensional energy space. The potentialenergy space introduced into this 2D energyspace model provides additional clues to improve the sampling efficiency. Based on this understanding, we propose two ways to calculate the freeenergy differences using the two schemes of the distribution method. We show that the distribution method implemented in the appropriate energy space—the energydifference space and the potentialenergy space, respectively—can improve the calculation of free energies in different circumstances. This analysis implies that the sampling can be improved if it is directed toward the appropriate region in the potentialenergy space, which is easily implemented in various types of freeenergy calculations. To test this, we calculate the freeenergy surface of alanine dipeptide in gas phase and in aqueous phase, respectively. We demonstrate that the freeenergy surface calculation is improved when the biased sampling of the potential energy is integrated into the sampling scheme.

A cutoff phenomenon in accelerated stochastic simulations of chemical kinetics via flow averaging (FLAVORSSA)
View Description Hide DescriptionWe present a simple algorithm for the simulation of stiff, discretespace, continuoustime Markov processes. The algorithm is based on the concept of flow averaging for the integration of stiff ordinary and stochastic differential equations and ultimately leads to a straightforward variation of the the wellknown stochastic simulation algorithm (SSA). The speedup that can be achieved by the present algorithm [flow averaging integrator SSA (FLAVORSSA)] over the classical SSA comes naturally at the expense of its accuracy. The error of the proposed method exhibits a cutoff phenomenon as a function of its speedup, allowing for optimal tuning. Two numerical examples from chemical kinetics are provided to illustrate the efficiency of the method.
 Gas Phase Dynamics and Structure: Spectroscopy, Molecular Interactions, Scattering, and Photochemistry

Ground and excited electronic states of azobenzene: A quantum Monte Carlo study
View Description Hide DescriptionLarge–scale quantum Monte Carlo (QMC) calculations of ground and excited singlet states of both conformers of azobenzene are presented. Remarkable accuracy is achieved by combining medium accuracy quantum chemistry methods with QMC. The results not only reproduce measured values with chemical accuracy but the accuracy is sufficient to identify part of experimental results which appear to be biased. Novel analysis of nodal surface structure yields new insights and control over their convergence, providing boost to the chemical accuracy electronic structure methods of large molecular systems.

Metastable anions of dinitrobenzene: Resonances for electron attachment and kinetic energy release
View Description Hide DescriptionAttachment of free, lowenergy electrons to dinitrobenzene (DNB) in the gas phase leads to DNB^{−} as well as several fragment anions. DNB^{−}, (DNBH)^{−}, (DNBNO)^{−}, (DNB2NO)^{−}, and (DNBNO_{2})^{−} are found to undergo metastable (unimolecular) dissociation. A rich pattern of resonances in the yield of these metastable reactions versus electron energy is observed; some resonances are highly isomerspecific. Most metastable reactions are accompanied by large average kinetic energy releases (KER) that range from 0.5 to 1.32 eV, typical of complex rearrangement reactions, but (1,3DNBH)^{−} features a resonance with a KER of only 0.06 eV for loss of NO. (1,3DNBNO)^{−} offers a rare example of a sequential metastable reaction, namely, loss of NO followed by loss of CO to yield C_{5}H_{4}O^{−} with a large KER of 1.32 eV. The G4(MP2) method is applied to compute adiabatic electron affinities and reaction energies for several of the observed metastable channels.

Microwave spectra, structure, and dynamics of the weakly bound complex, N_{2} CO_{2}
View Description Hide DescriptionThe Fourier transformmicrowave spectra of the various isotopologs of the weakly bound complex of carbon dioxide with the most abundant molecule in the atmosphere, nitrogen, have been measured. The structure of the complex has been determined and evidence for the inversion of the N_{2} is presented. The molecule is Tshaped, with the OCO forming the cross of the T, a structure consistent with that deduced from a previous rotationally resolved infrared experiment. A significant wideamplitude bending motion of the N_{2} is deduced from the values of the (nearly identical) nuclear quadrupole coupling constants of the nitrogen nuclei. The spectroscopic results are compared with highquality ab initio calculations. We examine the consequences of the N_{2} CO_{2} complex formation in the atmosphere upon the greenhouse warming potential of carbon dioxide.