Volume 140, Issue 18, 14 May 2014

Standard density functional approximations often give questionable results for oddelectron radical complexes, with the error typically attributed to selfinteraction. In density corrected density functional theory (DCDFT), certain classes of density functional theory calculations are significantly improved by using densities more accurate than the selfconsistent densities. We discuss how to identify such cases, and how DCDFT applies more generally. To illustrate, we calculate potential energy surfaces of HO·Cl^{−} and HO·H2O complexes using various common approximate functionals, with and without this density correction. Commonly used approximations yield wrongly shaped surfaces and/or incorrect minima when calculated self consistently, while yielding almost identical shapes and minima when density corrected. This improvement is retained even in the presence of implicit solvent.
 ARTICLES

 Theoretical Methods and Algorithms

An efficient algorithm for multipole energies and derivatives based on spherical harmonics and extensions to particle mesh Ewald
View Description Hide DescriptionNextgeneration molecular force fields deliver accurate descriptions of noncovalent interactions by employing more elaborate functional forms than their predecessors. Much work has been dedicated to improving the description of the electrostatic potential (ESP) generated by these force fields. A common approach to improving the ESP is by augmenting the point charges on each center with higherorder multipole moments. The resulting anisotropy greatly improves the directionality of the noncovalent bonding, with a concomitant increase in computational cost. In this work, we develop an efficient strategy for enumerating multipole interactions, by casting an efficient spherical harmonic based approach within a particle mesh Ewald (PME) framework. Although the derivation involves lengthy algebra, the final expressions are relatively compact, yielding an approach that can efficiently handle both finite and periodic systems without imposing any approximations beyond PME. Forces and torques are readily obtained, making our method well suited to modern molecular dynamics simulations.

Pathway structure determination in complex stochastic networks with nonexponential dwell times
View Description Hide DescriptionAnalysis of complex networks has been widely used as a powerful tool for investigating various physical, chemical, and biological processes. To understand the emergent properties of these complex systems, one of the most basic issues is to determine the structure and topology of the underlying networks. Recently, a new theoretical approach based on firstpassage analysis has been developed for investigating the relationship between structure and dynamic properties for network systems with exponential dwell time distributions. However, many real phenomena involve transitions with nonexponential waiting times. We extend the firstpassage method to uncover the structure of distinct pathways in complex networks with nonexponential dwell time distributions. It is found that the analysis of early time dynamics provides explicit information on the length of the pathways associated to their dynamic properties. It reveals a universal relationship that we have condensed in one general equation, which relates the number of intermediate states on the shortest path to the early time behavior of the firstpassage distributions. Our theoretical predictions are confirmed by extensive Monte Carlo simulations.

The RotnePragerYamakawa approximation for periodic systems in a shear flow
View Description Hide DescriptionRotnePragerYamakawa approximation is a commonly used approach to model hydrodynamic interactions between particles suspended in fluid. It takes into account all the longrange contributions to the hydrodynamic tensors, with the corrections decaying at least as fast as the inverse fourth power of the interparticle distances, and results in a positive definite mobility matrix, which is fundamental in Brownian dynamics simulations. In this communication, we show how to construct the RotnePragerYamakawa approximation for the bulk system under shear flow, which is modeled using the LeesEdwards boundary conditions.

General theory of multistage geminate reactions of isolated pairs of reactants. I. Kinetic equations
View Description Hide DescriptionGeneral matrix approach to the consideration of multistage geminate reactions of isolated pairs of reactants depending on reactant mobility is formulated on the basis of the concept of “effective” particles. Various elementary reactions (stages of multistage reaction including physicochemical processes of internal quantum state changes) proceeding with the participation of isolated pairs of reactants (or isolated reactants) are taken into account. Investigation has been made in terms of kinetic approach implying the derivation of general (matrix) kinetic equations for local and mean probabilities of finding any of the reaction species in the sample under study (or for local and mean concentrations). The recipes for the calculation of kinetic coefficients of the equations for mean quantities in terms of relative coordinates of reactants have been formulated in the general case of inhomogeneous reacting systems. Important specific case of homogeneous reacting systems is considered.

Evaluating interaction energies of weakly bonded systems using the BuckinghamHirshfeld method
View Description Hide DescriptionWe present the finalized BuckinghamHirshfeld method (BHDDFT) for the evaluation of interaction energies of nonbonded dimers with Density Functional Theory (DFT). In the method, dispersion energies are evaluated from static multipole polarizabilities, obtained onthefly from Coupled Perturbed KohnSham calculations and partitioned into diatomic contributions using the iterative Hirshfeld partitioning method. The dispersion energy expression is distributed over four atoms and has therefore a higher delocalized character compared to the standard pairwise expressions. Additionally, full multipolar polarizability tensors are used as opposed to effective polarizabilities, allowing to retain the anisotropic character at no additional computational cost. A density dependent damping function for the BLYP, PBE, BP86, B3LYP, and PBE0 functionals has been implemented, containing two global parameters which were fitted to interaction energies and geometries of a selected number of dimers using a bivariate RMS fit. The method is benchmarked against the S22 and S66 data sets for equilibrium geometries and the S22x5 and S66x8 data sets for interaction energies around the equilibrium geometry. Best results are achieved using the B3LYP functional with mean average deviation values of 0.30 and 0.24 kcal/mol for the S22 and S66 data sets, respectively. This situates the BHDDFT method among the best performing dispersion inclusive DFT methods. Effect of counterpoise correction on DFT energies is discussed.

Improving long time behavior of Poisson bracket mapping equation: A nonHamiltonian approach
View Description Hide DescriptionUnderstanding nonadiabatic dynamics in complex systems is a challenging subject. A series of semiclassical approaches have been proposed to tackle the problem in various settings. The Poisson bracket mapping equation (PBME) utilizes a partial Wigner transform and a mapping representation for its formulation, and has been developed to describe nonadiabatic processes in an efficient manner. Operationally, it is expressed as a set of Hamilton's equations of motion, similar to more conventional classical molecular dynamics. However, this original Hamiltonian PBME sometimes suffers from a large deviation in accuracy especially in the long time limit. Here, we propose a nonHamiltonian variant of PBME to improve its behavior especially in that limit. As a benchmark, we simulate spinboson and photosynthetic model systems and find that it consistently outperforms the original PBME and its Ehrenfest style variant. We explain the source of this improvement by decomposing the components of the mapping Hamiltonian and by assessing the energy flow between the system and the bath. We discuss strengths and weaknesses of our scheme with a viewpoint of offering future prospects.

Total photoionization crosssections of excited electronic states by the algebraic diagrammatic constructionStieltjesLanczos method
View Description Hide DescriptionHere, we extend the ab initio method for molecular photoionization crosssections introduced in Gokhberg et al. [J. Chem. Phys.130, 064104 (2009)] and benchmarked in Ruberti et al. [J. Chem. Phys.139, 144107 (2013)] to the calculation of total photoionization crosssections of molecules in electronically excited states. The method is based on the ab initio description of molecular electronic states within the manyelectron Green's function approach, known as algebraic diagrammatic construction (ADC), and on the application of StieltjesChebyshev moment theory to Lanczos pseudospectra of the ADC electronic Hamiltonian. The intermediate state representation of the dipole operator in the ADC basis is used to compute the transition moments between the excited states of the molecule. We compare the results obtained using different levels of the manybody theory, i.e., ADC(1), ADC(2), and ADC(2)x for the first two excited states of CO, N2, and H2O both at the ground state and the excited state equilibrium or saddle point geometries. We find that the single excitation ADC(1) method is not adequate even at the qualitative level and that the inclusion of double electronic excitations for description of excited state photoionization is essential. Moreover, we show that the use of the extended ADC(2)x method leads to a substantial systematic difference from the strictly secondorder ADC(2). Our calculations demonstrate that a theoretical modelling of photoionization of excited states requires an intrinsically double excitation theory with respect to the ground state and cannot be achieved by the standard single excitation methods with the ground state as a reference.

Compact twoelectron wave function for bond dissociation and Van der Waals interactions: A natural amplitude assessment
View Description Hide DescriptionElectron correlations in molecules can be divided in short range dynamical correlations, long range Van der Waals type interactions, and near degeneracy static correlations. In this work, we analyze for a onedimensional model of a twoelectron system how these three types of correlations can be incorporated in a simple wave function of restricted functional form consisting of an orbital product multiplied by a single correlation function f (r 12) depending on the interelectronic distance r 12. Since the three types of correlations mentioned lead to different signatures in terms of the natural orbital (NO) amplitudes in twoelectron systems, we make an analysis of the wave function in terms of the NO amplitudes for a model system of a diatomic molecule. In our numerical implementation, we fully optimize the orbitals and the correlation function on a spatial grid without restrictions on their functional form. Due to this particular form of the wave function, we can prove that none of the amplitudes vanishes and moreover that it displays a distinct sign pattern and a series of avoided crossings as a function of the bond distance in agreement with the exact solution. This shows that the wave function ansatz correctly incorporates the long range Van der Waals interactions. We further show that the approximate wave function gives an excellent binding curve and is able to describe static correlations. We show that in order to do this the correlation function f (r 12) needs to diverge for large r 12 at large internuclear distances while for shorter bond distances it increases as a function of r 12 to a maximum value after which it decays exponentially. We further give a physical interpretation of this behavior.

Stochastic quasisteady state approximations for asymptotic solutions of the chemical master equation
View Description Hide DescriptionIn this paper, we propose two methods to carry out the quasisteady state approximation in stochastic models of enzyme catalytic regulation, based on WKB asymptotics of the chemical master equation or of the corresponding partial differential equation for the generating function. The first of the methods we propose involves the development of multiscale generalisation of a WKB approximation of the solution of the master equation, where the separation of time scales is made explicit which allows us to apply the quasisteady state approximation in a straightforward manner. To the lowest order, the multiscale WKB method provides a quasisteady state, Gaussian approximation of the probability distribution. The second method is based on the HamiltonJacobi representation of the stochastic process where, as predicted by large deviation theory, the solution of the partial differential equation for the corresponding characteristic function is given in terms of an effective action functional. The optimal transition paths between two states are then given by those paths that maximise the effective action. Such paths are the solutions of the Hamilton equations for the Hamiltonian associated to the effective action functional. The quasisteady state approximation is applied to the Hamilton equations thus providing an approximation to the optimal transition paths and the transition time between two states. Using this approximation we predict that, unlike the meanfield quasisteady approximation result, the rate of enzyme catalysis depends explicitly on the initial number of enzyme molecules. The accuracy and validity of our approximated results as well as that of our predictions regarding the behaviour of the stochastic enzyme catalytic models are verified by direct simulation of the stochastic model using Gillespie stochastic simulation algorithm.

An ab initio approach to freeenergy reconstruction using logarithmic mean force dynamics
View Description Hide DescriptionWe present an ab initio approach for evaluating a free energy profile along a reaction coordinate by combining logarithmic mean force dynamics (LogMFD) and firstprinciples molecular dynamics. The mean force, which is the derivative of the free energy with respect to the reaction coordinate, is estimated using density functional theory (DFT) in the present approach, which is expected to provide an accurate free energy profile along the reaction coordinate. We apply this new method, firstprinciples LogMFD (FPLogMFD), to a glycine dipeptide molecule and reconstruct one and twodimensional free energy profiles in the framework of DFT. The resultant free energy profile is compared with that obtained by the thermodynamic integration method and by the previous LogMFD calculation using an empirical forcefield, showing that FPLogMFD is a promising method to calculate free energy without empirical forcefields.

Transformation of potential energy surfaces for estimating isotopic shifts in anharmonic vibrational frequency calculations
View Description Hide DescriptionA transformation of potential energy surfaces (PES) being represented by multimode expansions is introduced, which allows for the calculation of anharmonic vibrational spectra of any isotopologue from a single PES. This simplifies the analysis of infrared spectra due to significant CPUtime savings. An investigation of remaining deviations due to truncations and the socalled multilevel approximation is provided. The importance of vibrationalrotational couplings for small molecules is discussed in detail. In addition, an analysis is proposed, which provides information about the quality of the transformation prior to its execution. Benchmark calculations are provided for a set of small molecules.

Studying protein assembly with reversible Brownian dynamics of patchy particles
View Description Hide DescriptionAssembly of protein complexes like virus shells, the centriole, the nuclear pore complex, or the actin cytoskeleton is strongly determined by their spatial structure. Moreover, it is becoming increasingly clear that the reversible nature of protein assembly is also an essential element for their biological function. Here we introduce a computational approach for the Brownian dynamics of patchy particles with anisotropic assemblies and fully reversible reactions. Different particles stochastically associate and dissociate with microscopic reaction rates depending on their relative spatial positions. The translational and rotational diffusive properties of all protein complexes are evaluated onthefly. Because we focus on reversible assembly, we introduce a scheme which ensures detailed balance for patchy particles. We then show how the macroscopic rates follow from the microscopic ones. As an instructive example, we study the assembly of a pentameric ring structure, for which we find excellent agreement between simulation results and a macroscopic kinetic description without any adjustable parameters. This demonstrates that our approach correctly accounts for both the diffusive and reactive processes involved in protein assembly.

Complete spectrum of the infiniteU Hubbard ring using group theory
View Description Hide DescriptionWe present a full analytical solution of the multiconfigurational strongly correlated mixedvalence problem corresponding to the NHubbard ring filled with N−1 electrons, and infinite onsite repulsion. While the eigenvalues and the eigenstates of the model are known already, analytical determination of their degeneracy is presented here for the first time. The full solution, including degeneracy count, is achieved for each spin configuration by mapping the Hubbard model into a set of Hückelannulene problems for rings of variable size. The number and size of these effective Hückel annulenes, both crucial to obtain Hubbard states and their degeneracy, are determined by solving a wellknown combinatorial enumeration problem, the necklace problem for N−1 beads and two colors, within each subgroup of the C N−1 permutation group. Symmetryadapted solution of the necklace enumeration problem is finally achieved by means of the subduction of coset representation technique [S. Fujita, Theor. Chim. Acta76, 247 (1989)], which provides a general and elegant strategy to solve the onehole infiniteU Hubbard problem, including degeneracy count, for any ring size. The proposed group theoretical strategy to solve the infiniteU Hubbard problem for N−1 electrons is easily generalized to the case of arbitrary electron count L, by analyzing the permutation group C L and all its subgroups.

Using multiscale preconditioning to accelerate the convergence of iterative molecular calculations
View Description Hide DescriptionIterative procedures for optimizing properties of molecular models often converge slowly owing to the computational cost of accurately representing features of interest. Here, we introduce a preconditioning scheme that allows one to use a less expensive model to guide exploration of the energy landscape of a more expensive model and thus speed the discovery of locally stable states of the latter. We illustrate our approach in the contexts of energy minimization and the string method for finding transition pathways. The relation of the method to other multilevel simulation techniques and possible extensions are discussed.

Accelerating rare events while overcoming the lowbarrier problem using a temperature program
View Description Hide DescriptionWe present a hierarchical coarsegrained simulation technique called the temperature programmed molecular dynamics (TPMD) method for accelerating molecular dynamics (MD) simulations of rare events. The method is targeted towards materials where a system visits many times a collection of energy basins in the potential energy surface, called a superbasin, via lowbarrier moves before escaping to a new superbasin via a highbarrier move. The superbasin escape events are rare at the MD time scales. The lowbarrier moves become accessible to MD by employing a temperature program, i.e., the MD temperature changes during the simulation. Once a superbasin is detected, transitions within the superbasin are ignored, in effect causing coarsegraining of basins. The temperature program enables the system to escape from the superbasin with reduced computational cost thereby overcoming the “lowbarrier” problem. The main advantage of our approach is that the superbasintosuperbasin transitions are accurately obtained at the original temperature with a reasonable computational cost. We study surface diffusion in Ag/Ag(001) system and demonstrate the ability of the TPMD method to span a widerange of timescales.

An efficient extrapolation to the (T)/CBS limit
View Description Hide DescriptionWe extrapolate to the perturbative triples (T)/complete basis set (CBS) limit using double ζ basis sets without polarization functions (Wesleyan1Triples2ζ or “Wes1T2Z”) and triple ζ basis sets with a single level of polarization functions (Wesleyan1Triples3ζ or “Wes1T3Z”). These basis sets were optimized for 102 species representing the first two rows of the Periodic Table. The species include the entire set of neutral atoms, positive and negative atomic ions, as well as several homonuclear diatomic molecules, hydrides, rare gas dimers, polar molecules, such as oxides and fluorides, and a few transition states. The extrapolated Wes1T(2,3)Z triples energies agree with (T)/CBS benchmarks to within ±0.65 mE h, while the rms deviations of comparable model chemistries W1, CBSAPNO, and CBSQB3 for the same test set are ±0.23 mE h, ±2.37 mE h, and ±5.80 mE h, respectively. The Wes1T(2,3)Z triples calculation time for the largest hydrocarbon in the G2/97 test set, C6H5Me^{+}, is reduced by a factor of 25 when compared to W1. The costeffectiveness of the Wes1T(2,3)Z extrapolation validates the usefulness of the Wes1T2Z and Wes1T3Z basis sets which are now available for a more efficient extrapolation of the (T) component of any composite model chemistry.
 Advanced Experimental Techniques

Perturbation of nuclear spin polarizations in solid state NMR of nitroxidedoped samples by magicangle spinning without microwaves
View Description Hide DescriptionWe report solid state ^{13}C and ^{1}H nuclear magnetic resonance (NMR) experiments with magicangle spinning (MAS) on frozen solutions containing nitroxidebased paramagnetic dopants that indicate significant perturbations of nuclear spin polarizations without microwave irradiation. At temperatures near 25 K, ^{1}H and crosspolarized ^{13}C NMR signals from ^{15}N,^{13}Clabeled Lalanine in trinitroxidedoped glycerol/water are reduced by factors as large as six compared to signals from samples without nitroxide doping. Without MAS or at temperatures near 100 K, differences between signals with and without nitroxide doping are much smaller. We attribute most of the reduction of NMR signals under MAS near 25 K to nuclear spin depolarization through the crosseffect dynamic nuclear polarization mechanism, in which threespin flips drive nuclear polarizations toward equilibrium with spin polarization differences between electron pairs. When T1e is sufficiently long relative to the MAS rotation period, the distribution of electron spin polarization across the nitroxide electron paramagnetic resonance lineshape can be very different from the corresponding distribution in a static sample at thermal equilibrium, leading to the observed effects. We describe threespin and 3000spin calculations that qualitatively reproduce the experimental observations.
 Atoms, Molecules, and Clusters

The rate constant for radiative association of HF: Comparing quantum and classical dynamics^{a)}
View Description Hide DescriptionRadiative association for the formation of hydrogen fluoride through the A^{1}Π → X^{1}Σ^{+} and X^{1}Σ^{+} → X^{1}Σ^{+} transitions is studied using quantum and classical dynamics. The total thermal rate constant is obtained for temperatures from 10 K to 20 000 K. Agreement between semiclassical and quantum approaches is observed for the A^{1}Π → X^{1}Σ^{+}rate constant above 2000 K. The agreement is explained by the fact that the corresponding cross section is free of resonances for this system. At temperatures below 2000 K we improve the agreement by implementing a simplified semiclassical expression for the rate constant, which includes a quantum corrected pair distribution. The rate coefficient for the X^{1}Σ^{+} → X^{1}Σ^{+} transition is calculated using Breit–Wigner theory and a classical formula for the resonance and direct contributions, respectively. In comparison with quantum calculations the classical formula appears to overestimate the direct contribution to the rate constant by about 12% for this transition. Below about 450 K the resonance contribution is larger than the direct, and above that temperature the opposite holds. The biggest contribution from resonances is at the lowest temperature in the study, 10 K, where it is more than four times larger than the direct. Below 1800 K the radiative association rate constant due to X^{1}Σ^{+} → X^{1}Σ^{+} transitions dominates over A^{1}Π → X^{1}Σ^{+}, while above that temperature the situation is the opposite.

Highly correlated configuration interaction calculations on water with large orbital bases
View Description Hide DescriptionA priori selected configuration interaction (SCI) with truncation energy error [C. F. Bunge, J. Chem. Phys.125, 014107 (2006)] and CI by parts [C. F. Bunge and R. CarbóDorca, J. Chem. Phys.125, 014108 (2006)] are used to approximate the total nonrelativistic electronic ground state energy of water at fixed experimental geometry with CI up to sextuple excitations. Correlationconsistent polarized corevalence basis sets (ccpCVnZ) up to sextuple zeta and augmented correlationconsistent polarized corevalence basis sets (augccpCVnZ) up to quintuple zeta quality are employed. Truncation energy errors range between less than 1 μhartree, and 100 μhartree for the largest orbital set. Coupled cluster CCSD and CCSD(T) calculations are also obtained for comparison. Our best upper bound, −76.4343 hartree, obtained by SCI with up to sextuple excitations with a ccpCV6Z basis recovers more than 98.8% of the correlation energy of the system, and it is only about 3 kcal/mol above the “experimental” value. Despite that the present energy upper bounds are far below all previous ones, comparatively large dispersion errors in the determination of the extrapolated energies to the complete basis set do not allow to determine a reliable estimation of the full CI energy with an accuracy better than 0.6 mhartree (0.4 kcal/mol).

Elastic scattering of lowenergy electrons by 1,4dioxane
View Description Hide DescriptionWe report calculated cross sections for elastic collisions of lowenergyelectrons with 1,4dioxane. Our calculations employed the Schwinger multichannel method with pseudopotentials and were carried out in the staticexchange and staticexchange plus polarization approximations for energies up to 30 eV. Our results show the presence of three shape resonances belonging to the B u , A u , and B g symmetries and located at 7.0 eV, 8.4 eV, and 9.8 eV, respectively. We also report the presence of a RamsauerTownsend minimum located at around 0.05 eV. We compare our calculated cross sections with experimental data and Rmatrix and independent atom model along with the additivity rule corrected by using screening coefficients theoretical results for 1,4dioxane obtained by Palihawadana et al. [J. Chem. Phys.139, 014308 (2013)]. The agreement between the present and the Rmatrix theoretical calculations of Palihawadana et al. is relatively good at energies below 10 eV. Our calculated differential cross sections agree well with the experimental data, showing only some discrepancies at higher energies.