Volume 135, Issue 8, 28 August 2011

A series of quantum molecular dynamics simulations have been performed to investigate the energetic, structural, dynamic, and spectroscopic properties of methanol cluster anions, [(CH_{3}OH)_{ n }]^{−}, (n = 50–500). Consistent with the inference from photoelectron imaging experiments, we find two main localization modes of the excess electron in equilibrated methanol clusters at ∼200 K. The two different localization patterns have strikingly different physical properties, consistent with experimental observations, and are manifest in comparable cluster sizes to those observed. Smaller clusters (n ≤ 128) tend to localize the electron in very weakly bound, diffuse electronic states on the surface of the cluster, while in larger ones the electron is stabilized in solvent cavities, in compact interiorbound states. The interior states exhibit properties that largely resemble and smoothly extrapolate to those simulated for a solvated electron in bulk methanol. The surface electronic states of methanol cluster anions are significantly more weakly bound than the surface states of the anionic water clusters. The key source of the difference is the lack of stabilizing free hydroxyl groups on a relaxed methanol clustersurface. We also provide a mechanistic picture that illustrates the essential role of the interactions of the excess electron with the hydroxyl groups in the dynamic process of the transition of the electron from surfacebound states to interiorbound states.
 ARTICLES

 Theoretical Methods and Algorithms

Reciprocity in the degeneracies of some tetraatomic molecular ions
View Description Hide DescriptionVarious ab initio computations, as, e.g., in G. J. Halász and Á. Vibók, Int. J. Quantum Chem.111, 342 (2011), have shown that in molecules of the type (HCCH)^{+}, when the extremal H atoms are distorted from a linear form but maintain a planar geometry,a pair of conical intersections (ci) occur at such positions that the ratios of the distortional coordinates of the two atoms are in the two ci's reciprocals of each other. These computations have here been extended to locate the ci's also for HCNH. The two groups of results are explained by simple analytic perturbational expressions for the energy differences of the lowest adjacent electronic states, with inclusion of excited state effects.

An orbitalinvariant and strictly size extensive postHartreeFock correlation functional
View Description Hide DescriptionA strictly size extensive postHartreeFock correlation functional being invariant with respect to orbital transformations within the occupied and virtual subspaces is presented. While avoiding the necessity to solve additional Z vector equations for the calculation of properties and energy gradients, this functional reproduces almost exactly the results of coupledcluster singles doubles (CCSD) calculations. In particular, it is demonstrated that the method is rigorous in the sense that it can be systematically improved by the perturbative inclusion of triple excitations in the same way as CCSD. As to the computational cost, the presented approach is somewhat more expensive than the CCSD if the energy is variationally optimized with respect to both the orbitals and the excitation amplitudes. Replacement of orbital optimization by the Brueckner condition reduces the computational cost by a factor of two, thus making the method less expensive than CCSD.

How accurate are the nonlinear chemical FokkerPlanck and chemical Langevin equations?
View Description Hide DescriptionThe chemical FokkerPlanck equation and the corresponding chemical Langevin equation are commonly used approximations of the chemical master equation. These equations are derived from an uncontrolled, secondorder truncation of the KramersMoyal expansion of the chemical master equation and hence their accuracy remains to be clarified. We use the systemsize expansion to show that chemical FokkerPlanck estimates of the mean concentrations and of the variance of the concentration fluctuations about the mean are accurate to order Ω^{−3/2} for reactionsystems which do not obey detailed balance and at least accurate to order Ω^{−2} for systems obeying detailed balance, where Ω is the characteristic size of the system. Hence, the chemical FokkerPlanck equation turns out to be more accurate than the linearnoise approximation of the chemical master equation (the linear FokkerPlanck equation) which leads to mean concentration estimates accurate to order Ω^{−1/2} and variance estimates accurate to order Ω^{−3/2}. This higher accuracy is particularly conspicuous for chemical systems realized in small volumes such as biochemical reactions inside cells. A formula is also obtained for the approximate size of the relative errors in the concentration and variance predictions of the chemical FokkerPlanck equation, where the relative error is defined as the difference between the predictions of the chemical FokkerPlanck equation and the master equation divided by the prediction of the master equation. For dimerization and enzymecatalyzed reactions, the errors are typically less than few percent even when the steadystate is characterized by merely few tens of molecules.

Breaking the carbon dimer: The challenges of multiple bond dissociation with full configuration interaction quantum Monte Carlo methods
View Description Hide DescriptionThe full configuration interaction quantum Monte Carlo (FCIQMC) method, as well as its “initiator” extension (iFCIQMC), is used to tackle the complex electronic structure of the carbon dimer across the entire dissociation reaction coordinate, as a prototypical example of a strongly correlated molecular system. Various basis sets of increasing size up to the large ccpVQZ are used, spanning a fully accessible Nelectron basis of over 10^{12} Slater determinants, and the accuracy of the method is demonstrated in each basis set. Convergence to the FCI limit is achieved in the largest basis with only walkers within random errorbars of a few tenths of a millihartree across the binding curve, and extensive comparisons to FCI, CCSD(T), MRCI, and CEEIS results are made where possible. A detailed exposition of the convergence properties of the FCIQMC methods is provided, considering convergence with elapsed imaginary time, number of walkers and size of the basis. Various symmetries which can be incorporated into the stochastic dynamic, beyond the standard abelian point group symmetry and spin polarisation are also described. These can have significant benefit to the computational effort of the calculations, as well as the ability to converge to various excited states. The results presented demonstrate a new benchmark accuracy in basisset energies for systems of this size, significantly improving on previous state of the art estimates.

An approximate densityfunctional method using the HarrisFoulkes functional
View Description Hide DescriptionWe present a method which uses the results of a molecular KohnSham calculation at a reference geometry to approximate the energy at many different geometries. The KohnSham electron density of the reference geometry is decomposed into atomic fragments, which move with the nuclei to approximate the density at a new geometry and the energy is evaluated with the HarrisFoulkes functional. Preliminary results for a biological quantummechanics/molecularmechanics trajectory are promising: the errors of referencegeometry HarrisFoulkes (compared to full selfconsistent KohnSham) for the PBE exchangecorrelation functional have the same magnitude as the difference between the energies of PBE and BLYP.

On the accuracy of the state space restriction approximation for spin dynamics simulations
View Description Hide DescriptionWe present an algebraic foundation for the state space restriction approximation in spin dynamics simulations and derive applicability criteria as well as minimal basis set requirements for practically encountered simulation tasks. The results are illustrated with nuclear magnetic resonance(NMR),electron spin resonance(ESR),dynamic nuclear polarization (DNP), and spin chemistry simulations. It is demonstrated that state space restriction yields accurate results in systems where the time scale of spin relaxation processes approximately matches the time scale of the experiment. Rigorous error bounds and basis set requirements are derived.

Incorporation of charge transfer into the explicit polarization fragment method by grand canonical density functional theory
View Description Hide DescriptionMolecular fragmentation algorithms provide a powerful approach to extending electronic structure methods to very large systems. Here we present a method for including charge transfer between molecular fragments in the explicit polarization (XPol) fragment method for calculating potential energy surfaces. In the conventional XPol method, the total charge of each fragment is preserved, and charge transfer between fragments is not allowed. The description of charge transfer is made possible by treating each fragment as an open system with respect to the number of electrons. To achieve this, we applied Mermin's finite temperature method to the XPol wave function. In the application of this method to XPol, the fragments are open systems that partially equilibrate their number of electrons through a quasithermodynamics electron reservoir. The number of electrons in a given fragment can take a fractional value, and the electrons of each fragment obey the Fermi–Dirac distribution. The equilibrium state for the electrons is determined by electronegativity equalization with conservation of the total number of electrons. The amount of charge transfer is controlled by reinterpreting the temperature parameter in the Fermi–Dirac distribution function as a coupling strength parameter. We determined this coupling parameter so as to reproduce the charge transferenergy obtained by block localized energy decomposition analysis. We apply the new method to ten systems, and we show that it can yield reasonable approximations to potential energy profiles, to charge transfer stabilization energies, and to the direction and amount of charge transferred.

The theoretical currentvoltage dependence of a nondegenerate disordered organic material obtained with conductive atomic force microscopy
View Description Hide DescriptionWe develop a simple continuum model for the current voltage characteristics of a material as measured by the conducting atomic force microscopy, including space charge effects. We address the effect of the point contact on the magnitude of the current and on the transition voltages between the different current regimes by comparing these with the corresponding expressions obtained with planar electrodes.

Using the chargestabilization technique in the double ionization potential equationofmotion calculations with dianion references
View Description Hide DescriptionThe chargestabilization method is applied to double ionization potential equationofmotion (EOMDIP) calculations to stabilize unstable dianion reference functions. The autoionizing character of the dianionic reference states spoils the numeric performance of EOMDIP limiting applications of this method. We demonstrate that reliable excitation energies can be computed by EOMDIP using a stabilized resonance wave function instead of the lowest energy solution corresponding to the neutral + free electron(s) state of the system. The details of chargestabilization procedure are discussed and illustrated by examples. The choice of optimal stabilizing Coulomb potential, which is strong enough to stabilize the dianion reference, yet, minimally perturbs the target states of the neutral, is the crux of the approach. Two algorithms of choosing optimal parameters of the stabilization potential are presented. One is based on the orbital energies, and another – on the basis set dependence of the total HartreeFock energy of the reference. Our benchmark calculations of the singlettriplet energy gaps in several diradicals show a remarkable improvement of the EOMDIP accuracy in problematic cases. Overall, the excitation energies in diradicals computed using the stabilized EOMDIP are within 0.2 eV from the reference EOM spinflip values.

Least constraint approach to the extraction of internal motions from molecular dynamics trajectories of flexible macromolecules
View Description Hide DescriptionWe propose a rigorous method for removing rigidbody motions from a given molecular dynamics trajectory of a flexible macromolecule. The method becomes exact in the limit of an infinitesimally small sampling step for the input trajectory. In a recent paper [G. Kneller, J. Chem. Phys.128, 194101 (2008)]10.1063/1.2902290, one of us showed that virtual internal atomic displacements for small time increments can be derived from Gauss’ principle of least constraint, which leads to a rotational superposition problem for the atomic coordinates in two consecutive time frames of the input trajectory. Here, we demonstrate that the accumulation of these displacements in a molecularfixed frame, which evolves in time according to the virtual rigidbody motions, leads to the desired trajectory for internal motions. The atomic coordinates in the input and output trajectory are related by a rototranslation, which guarantees that the internal energy of the molecule is left invariant. We present a convenient implementation of our method, in which the accumulation of the internal displacements is performed implicitly. Two numerical examples illustrate the difference to the classical approach for removing macromolecular rigidbody motions, which consists of aligning its configurations in the input trajectory with a fixed reference structure.

Extension of the invariant environment refinement technique + reverse Monte Carlo method of structural modelling for interpreting experimental structure factors: The cases of amorphous silicon, phosphorus, and liquid argon
View Description Hide DescriptionThe invariant environment refinement technique, as applied to reverse Monte Carlomodelling [invariant environment refinement technique + reverse Monte Carlo (INVERT + RMC); M. J. Cliffe, M. T. Dove, D. A. Drabold, and A. L. Goodwin, Phys. Rev. Lett.104, 125501 (2010)10.1103/PhysRevLett.104.125501], is extended so that it is now applicable for interpreting the structure factor (instead of the pair distribution function). The new algorithm, called the local invariance calculation, is presented by the examples of amorphous silicon, phosphorus, and liquid argon. As a measure of the effectiveness of the new algorithm, the ratio of exactly fourfold coordinated Si atoms was larger than obtained previously by the INVERTRMC scheme.

Calculation of multiple initial state selected reaction probabilities from Chebyshev fluxflux correlation functions: Influence of reactant internal excitations on H + H_{2}O → OH + H_{2}
View Description Hide DescriptionA Chebyshevbased fluxflux correlation function approach is introduced for calculating multiple initial state selected reaction probabilities for bimolecular reactions. Based on the quantum transitionstate theory, this approach propagates, with the exact Chebyshev propagator, transitionstate wave packets towards the reactant asymptote. It is accurate and efficient if many initial state selected reaction probabilities are needed. This approach is applied to the title reaction to elucidate the influence of the H_{2}O rovibrational states on its reactivity. Results from several potential energy surfaces are compared.

Pair correlation function integrals: Computation and use
View Description Hide DescriptionWe describe a method for extending radial distribution functions obtained from molecular simulations of pure and mixed molecular fluids to arbitrary distances. The method allows total correlation function integrals to be reliably calculated from simulations of relatively small systems. The longdistance behavior of radial distribution functions is determined by requiring that the corresponding direct correlation functions follow certain approximations at long distances. We have briefly described the method and tested its performance in previous communications [R. Wedberg, J. P. O’Connell, G. H. Peters, and J. Abildskov, Mol. Simul.36, 1243 (2010);10.1080/08927020903536366Fluid Phase Equilib.302, 32 (2011)]10.1016/j.fluid.2010.10.004, but describe here its theoretical basis more thoroughly and derive longdistance approximations for the direct correlation functions. We describe the numerical implementation of the method in detail, and report numerical tests complementing previous results. Pure molecular fluids are here studied in the isothermalisobaric ensemble with isothermal compressibilities evaluated from the total correlation function integrals and compared with values derived from volume fluctuations. For systems where the radial distribution function has structure beyond the sampling limit imposed by the system size, the integration is more reliable, and usually more accurate, than simple integral truncation.

Analytic energy gradients for the spinfree exact twocomponent theory using an exact block diagonalization for the oneelectron Dirac Hamiltonian
View Description Hide DescriptionWe report the implementation of analytic energy gradients for the evaluation of firstorder electrical properties and nuclear forces within the framework of the spinfree (SF) exact twocomponent (X2c) theory. In the scheme presented here, referred to in the following as SFX2c1e, the decoupling of electronic and positronic solutions is performed for the oneelectron Dirac Hamiltonian in its matrix representation using a single unitary transformation. The resulting twocomponent oneelectron matrix Hamiltonian is combined with untransformed twoelectron interactions for subsequent selfconsistentfield and electroncorrelated calculations. The “picturechange” effect in the calculation of properties is taken into account by considering the full derivative of the twocomponent Hamiltonian matrix with respect to the external perturbation. The applicability of the analyticgradient scheme presented here is demonstrated in benchmark calculations. SFX2c1e results for the dipole moments and electricfield gradients of the hydrogen halides are compared with those obtained from nonrelativistic, SF highorder DouglasKrollHess, and SF DiracCoulomb calculations. It is shown that the use of untransformed twoelectron interactions introduces rather small errors for these properties. As a first application of the analytic geometrical gradient, we report the equilibrium geometry of methylcopper (CuCH_{3}) determined at various levels of theory.

Velocityscaling optimized replica exchange molecular dynamics of proteins in a hybrid explicit/implicit solvent
View Description Hide DescriptionWe propose a scheme for replica exchange molecular dynamics of proteins in explicit solvent that minimizes the number of required replicas using velocity rescaling. Our approach relies on a hybrid method where the protein evolves at each temperature in an explicit solvent, but replica exchange moves utilize an implicit solvent term. The two terms are coupled through the velocity rescaling. We test the efficiency of this approach for a common test case, the trpcage protein.

Comparison of Brownian dynamics algorithms with hydrodynamic interaction
View Description Hide DescriptionThe hydrodynamic interaction is an essential effect to consider in Brownian dynamics simulations of polymer and nanoparticle dilute solutions. Several mathematical approaches can be used to build Brownian dynamics algorithms with hydrodynamic interaction, the most common of them being the exact but time demanding Cholesky decomposition and the Chebyshev polynomial expansion. Recently, Geyer and Winter [J. Chem. Phys.130, 1149051 (2009)]10.1063/1.3089668 have proposed a new approximation to treat the hydrodynamic interaction that seems quite efficient and is increasingly used. So far, a systematic comparison among those approaches has not been clearly made. In this paper, several features and the efficiency of typical implementations of those approaches are evaluated by using beadandspring chain models. The different sensitivity to the bead overlap detected for the different implementations may be of interest to select the suitable algorithm for a given simulation.

A systematic formulation of the virial expansion for nonadditive interaction potentials
View Description Hide DescriptionA new formulation of the virial expansion for a classical gas is derived without the restriction to pairwiseadditive interaction potentials. Explicit expressions up to the eighth virial coefficient, suitable for numerical evaluation, are given in the form of integrals over sums of cluster diagrams. Compared with previous approaches, the number of cluster diagrams increases more moderately with increasing order of the virial coefficient. Thus, the new formulation should be particularly useful for the computation of highorder virial coefficients.

Application of an efficient multireference approach to freebase porphin and metalloporphyrins: Ground, excited, and positive ion states
View Description Hide DescriptionThe improved virtual orbitalcomplete active space configuration interaction (IVOCASCI) method is applied to determine the geometries of the ground state of freebase porphin and its metal derivatives, magnesium and zincporphyrins. The vertical excitation energies and ionization potentials are computed at these optimized geometries using an IVObased version of multireference MöllerPlesset (IVOMRMP) perturbation theory. The geometries and excitation energies obtained from the IVOCASCI and IVOMRMP methods agree well with experiment and with other correlated manybody methods. We also provide the ground state vibrational frequencies for freebase porphin and Mgporphyrin. All frequencies are real in contrast to selfconsistent field treatments which yield an imaginary frequency. Ground state normal mode frequencies (scaled) of freebase porphin and magnesiumporphyrin from IVOCASCI and complete active space selfconsistent field methods are quite similar and are consistent with BeckeSlaterHartreeFock exchange and LeeYangParr correlationdensity functional theory calculations and with experiment. In addition, geometries are determined for lowlying excited state triplets and for positive ion states of the molecules. To our knowledge, no prior experimental and theoretical data are available for these excited state geometries of magnesium and zincporphyrins. Given that the IVOCASCI and IVOMRMP computed geometries and excitation energies agree favorably with experiment and with available theoretical data, our predicted excited state geometries should be equally accurate.

Closedshell ring coupled cluster doubles theory with range separation applied on weak intermolecular interactions
View Description Hide DescriptionWe explore different variants of the random phase approximation to the correlation energy derived from closedshell ringdiagram approximations to coupled cluster doubles theory. We implement these variants in rangeseparated densityfunctional theory, i.e., by combining the longrange random phase approximations with shortrange densityfunctional approximations. We perform tests on the raregas dimers He_{2}, Ne_{2}, and Ar_{2}, and on the weakly interacting molecular complexes of the S22 set of Jurečka et al.[P. Jurečka, J. Šponer, J. Černý, and P. Hobza, Phys. Chem. Chem. Phys.8, 1985 (2006)10.1039/b600027d]. The two best variants correspond to the ones originally proposed by Szabo and Ostlund [A. Szabo and N. S. Ostlund, J. Chem. Phys.67, 4351 (1977)10.1063/1.434580]. With range separation, they reach mean absolute errors on the equilibrium interaction energies of the S22 set of about 0.4 kcal/mol, corresponding to mean absolute percentage errors of about 4%, with the augccpVDZ basis set.

Selfconsistency in frozendensity embedding theory based calculations
View Description Hide DescriptionThe bifunctional for the nonelectrostatic part of the exact embedding potential of frozendensity embedding theory (FDET) depends on whether the embedded part is described by means of a real interacting manyelectron system or the reference system of noninteracting electrons (see [Wesolowski, Phys. Rev. A.77, 11444 (2008)]). The difference , where ΔF ^{ MD }[ρ_{ A }] is the functional bound from below by the correlation functional E _{ c }[ρ_{ A }] and from above by zero. Taking into account ΔF ^{ MD }[ρ_{ A }] in both the embedding potential and in energy is indispensable for assuring that all calculated quantities are selfconsistent and that FDET leads to the exact energy and density in the limit of exact functionals. Since not much is known about good approximations for ΔF ^{ MD }[ρ_{ A }], we examine numerically the adequacy of neglecting ΔF ^{ MD }[ρ_{ A }] entirely. To this end, we analyze the significance of in the case where the magnitude of ΔF ^{ MD }[ρ_{ A }] is the largest, i.e., for HartreeFock wavefunction. In hydrogen bonded model systems, neglecting in the embedding potential marginally affects the total energy (less than 5% change in the interaction energy) but results in qualitative changes in the calculated hydrogenbonding induced shifts of the orbital energies. Based on this estimation, we conclude that neglecting may represent a good approximation for multireference variational methods using adequate choice for the active space. Doing the same for singlereference perturbative methods is not recommended. Not only it leads to violation of selfconsistency but might result in large effect on orbital energies. It is shown also that the errors in total energy due to neglecting do not cancel but rather add up to the errors due to approximation for the bifunctional of the nonadditive kinetic potential.