Volume 138, Issue 10, 14 March 2013
 ARTICLES

 Theoretical Methods and Algorithms

Appraisal of molecular tailoring approach for large clusters
View Description Hide DescriptionHigh level ab initio investigations on molecular clusters are generally restricted to those of small size essentially due to the nonlinear scaling of corresponding computational cost. Molecular tailoring approach (MTA) is a fragmentationbased method, which offers an economical and efficient route for studying larger clusters. However, due to its approximate nature, the MTAenergies carry some errors visàvis their full calculation counterparts. These errors in the MTAenergies are reduced by grafting the correction at a lower basis set (e.g., 631+G(d)) onto a higher basis set (e.g., augccpvdz or augccpvtz) calculation at MP2 level of theory. Further, better estimates of energies are obtained by making use of manybody interaction analysis. For this purpose, Rgoodness (Rg) parameters for the three and fourbody interactions in a fragmentation scheme are proposed. The procedure employing grafting and manybody analysis has been tested out on molecular clusters of water, benzene, acetylene and carbon dioxide. It is found that for the fragmentation scheme having higher three and fourbody Rgvalues, the errors in MTAgrafted energies are reduced typically to ∼0.2 mH at MP2 level calculation. Coupled with the advantage in terms of computational resources and CPU time, the present method opens a possibility of accurate treatment of large molecular clusters.

Derivation and assessment of relativistic hyperfinecoupling tensors on the basis of orbitaloptimized secondorder Møller–Plesset perturbation theory and the secondorder Douglas–Kroll–Hess transformation
View Description Hide DescriptionThe accurate calculation of hyperfinecoupling tensors requires a good description of the electronic spin density, especially close to and at the nucleus. Thus, dynamic correlation as well as relativistic effects have to be included in the quantumchemical calculation of this quantity. In this paper, orbitaloptimized secondorder Møller–Plesset perturbation theory (MP2) is combined with the secondorder Douglas–Kroll–Hess (DKH) transformation to yield an efficient and accurate ab initio method for the calculation of hyperfine couplings for larger molecules including heavy elements. Particular attention is paid to the derivation of the hyperfinecoupling tensor in the DKH framework. In the presence of a magnetic field, the DKHtransformation is not unique. Two different versions can be found in the literature. In this paper, a detailed derivation of oneelectron contributions to the hyperfinecoupling tensor as they arise in linearresponse theory is given for both DKHtransformations. It turns out that one of the two variants produces divergent hyperfinecoupling constants. The possibility to remove this divergence through a physically motivated finitenucleus model taking into account the different extent of charge and magnetization distribution is discussed. Hyperfinecoupling values obtained at the orbitaloptimized MP2 level with secondorder DKH corrections for the nondivergent variant are presented. The influence of a Gaussian nucleus model is studied. The method is compared to fourcomponent, highaccuracy calculations for a number of cations and atoms. Comparison to B3LYP and B2PLYP is made for a set of transitionmetal complexes of moderate size.

Nonadiabaticity in a JahnTeller system probed by absorption and resonance Raman scattering
View Description Hide DescriptionA theory of absorption and resonance Raman scattering of impurity centers in crystals with E⊗etype JahnTeller effect in the excited state is presented. The vibronic interaction with nontotally symmetric local or pseudolocal modes and with a continuum of bath modes (phonons) is considered. A number of specific quantum effects, such as the nonadiabaticityinduced enhancement of the Raman scattering at highenergy excitation, the size effect of the final state, the interference of different channels of scattering, the Fermi resonances in the conical intersection, and others, were shown to become apparent in the calculated spectra. The vibronic interaction with phonons essentially determines the structure of the spectra.

Analytical energy gradients for secondorder multireference perturbation theory using density fitting
View Description Hide DescriptionWe present algorithms for computing analytical energy gradients for multiconfiguration selfconsistent field methods and partially internally contracted complete active space secondorder perturbation theory (CASPT2) using density fitting (DF). Our implementation is applicable to both singlestate and multistate CASPT2 analytical gradients. The accuracy of the new methods is demonstrated for structures and excitation energies of valence and Rydberg states of pyrrole, as well as for structures and adiabatic singlettriplet energy splittings for the hydro, the O,O^{′}formato, and the N,N^{′}diiminatocopperdioxygen complexes. It is shown that the effects of density fitting on optimized structures and relative energies are negligible. For cases in which the total cost is dominated by the integral evaluations and transformations, the DFCASPT2 gradient calculations are found to be faster than the corresponding conventional calculations by typically a factor of three to five using tripleζ basis sets, and by about a factor of ten using quadrupleζ basis sets.

Thermal activation at moderatetohigh and high damping: Finite barrier effects and force spectroscopy
View Description Hide DescriptionWe study the thermal escape problem in the moderatetohigh and high damping regime of a system with a parabolic barrier. We present a formula that matches our numerical results accounting for finite barrier effects, and compare it with previous works. We also show results for the full damping range. We quantitatively study some aspects on the relation between mean first passage time and the definition of an escape rate. To finish, we apply our results and considerations in the framework of force spectroscopy problems. We study the differences on the predictions using the different theories and discuss the role of as the relevant parameter at high damping.

Computing rovibrational levels of methane with curvilinear internal vibrational coordinates and an Eckart frame
View Description Hide DescriptionWe present a new procedure for computing a rovibrational spectrum of a polyatomic molecule and apply it to methane. The Schrödinger equation is solved, numerically exactly, by using a nested contracted basis. Rovibrational wavefunctions are computed in a v⟩JKM⟩ basis, where v⟩ is a vibrational wavefunction and JKM⟩ is a symmetric top wavefunction. In turn, the v⟩ are obtained by solving a vibrational Schrödinger equation with basis functions that are products of contracted bend and stretch functions. At all stages of the calculation we exploit parity symmetry. The calculations are done in internal coordinates that facilitate the treatment of large amplitude motion. An Eckart moleculefixed frame is used by numerically computing coefficients of the kinetic energy operator. The efficacy of the method is demonstrated by calculating a large number of converged J = 10 methane rovibrational levels in the Tetradecad polyad. No previous calculation of rovibrational levels of methane includes as many levels as we report in this paper.

An algorithm for quantum mechanical finitenuclearmass variational calculations of atoms with L = 3 using allelectron explicitly correlated Gaussian basis functions
View Description Hide DescriptionA new algorithm for quantummechanical nonrelativistic calculation of the Hamiltonian matrix elements with allelectron explicitly correlated Gaussian functions for atoms with an arbitrary number of s electrons and with three p electrons, or one p electron and one d electron, or one f electron is developed and implemented. In particular the implementation concerns atomic states with L = 3 and M = 0. The Hamiltonian used in the approach is obtained by rigorously separating the centerofmass motion from the laboratoryframe all particle Hamiltonian, and thus it explicitly depends on the finite mass of the nucleus. The approach is employed to perform test calculations on the lowest ^{2}F state of the two main isotopes of the lithium atom, ^{7}Li and ^{6}Li.

Adaptive molecular decomposition: Largescale quantum chemistry for liquids
View Description Hide DescriptionWe present a linearscaling method based on selfconsistent charge nonorthogonal tightbinding. Linear scaling is achieved using a manybody expansion, which is adjusted dynamically to the instantaneous molecular configuration of a liquid. The method is capable of simulating liquids over large length and time scales, and also handles reactions correctly. Benchmarking on typical carbonate electrolytes used in Liion batteries displays excellent agreement with results from full tightbinding calculations. The decomposition slightly breaks the HellmannFeynman theorem, which is demonstrated by application to water. However, an additional correction also enables dynamical simulation in this case.

Normconserving pseudopotentials with chemical accuracy compared to allelectron calculations
View Description Hide DescriptionBy adding a nonlinear core correction to the well established dual space Gaussian type pseudopotentials for the chemical elements up to the third period, we construct improved pseudopotentials for the PerdewBurkeErnzerhof [J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett.77, 3865 (Year: 1996)10.1103/PhysRevLett.77.3865] functional and demonstrate that they exhibit excellent accuracy. Our benchmarks for the G21 test set show average atomization energy errors of only half a kcal/mol. The pseudopotentials also remain highly reliable for high pressure phases of crystalline solids. When supplemented by empirical dispersion corrections [S. Grimme, J. Comput. Chem.27, 1787 (Year: 2006)10.1002/jcc.20495; S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys.132, 154104 (Year: 2010)10.1063/1.3382344] the average error in the interaction energy between molecules is also about half a kcal/mol. The accuracy that can be obtained by these pseudopotentials in combination with a systematic basis set is well superior to the accuracy that can be obtained by commonly used medium size Gaussian basis sets in allelectron calculations.

A Cartesian quasiclassical model to nonequilibrium quantum transport: The Anderson impurity model
View Description Hide DescriptionWe apply the recently proposed quasiclassical approach for a second quantized manyelectron Hamiltonian in Cartesian coordinates [B. Li and W. H. Miller, J. Chem. Phys.137, 154107 (Year: 2012)10.1063/1.4757935] to correlated nonequilibrium quantum transport. The approach provides accurate results for the resonant level model for a wide range of temperatures, bias, and gate voltages, correcting the flaws of our recently proposed mapping using actionangle variables. When electronelectron interactions are included, a Gaussian function scheme is required to map the twoelectron integrals, leading to quantitative results for the Anderson impurity model. In particular, we show that the current mapping is capable of capturing quantitatively the Coulomb blockade effect and the temperature dependence of the current below and above the blockade.

Mathematics of small stochastic reaction networks: A boundary layer theory for eigenstate analysis
View Description Hide DescriptionWe study and analyze the stochastic dynamics of a reversible bimolecular reaction A + B ↔ C called the “trivalent reaction.” This reaction is of a fundamental nature and is part of many biochemical reaction networks. The stochastic dynamics is given by the stochastic master equation, which is difficult to solve except when the equilibrium state solution is desired. We present a novel way of finding the eigenstates of this system of differencedifferential equations, using perturbation analysis of ordinary differential equations arising from approximation of the difference equations. The time evolution of the state probabilities can then be expressed in terms of the eigenvalues and the eigenvectors.

Theory of reversible diffusioninfluenced reactions with nonMarkovian dissociation in two space dimensions
View Description Hide DescriptionWe investigate the reversible diffusioninfluenced reaction of an isolated pair in the presence of a nonMarkovian generalization of the backreaction boundary condition in two space dimensions. Following earlier work by Agmon and Weiss, we consider residence time probability densities that decay slower than an exponential and that are characterized by a single parameter 0 < σ ⩽ 1. We calculate an exact expression for a Green's function of the twodimensional diffusion equation subject to a nonMarkovian backreaction boundary condition that is valid for arbitrary σ and for all times. We use the obtained expression to derive the survival probability for the initially unbound pair and we calculate an exact expression for the probability S(t*) that the initially bound particle is unbound. Finally, we obtain an approximate solution for long times. In particular, we show that the ultimate fate of the bound state is complete dissociation, as in the Markovian case. However, the limiting value is approached quite differently: Instead of a ∼t ^{−1} decay, we obtain 1 − S(t*) ∼ t ^{−σ}ln t. The derived expressions should be relevant for a better understanding of reversible membranebound reactions in cell biology.

A fully variational spinorbit coupled complete active space selfconsistent field approach: Application to electron paramagnetic resonance gtensors
View Description Hide DescriptionIn this work, a relativistic version of the stateaveraged complete active space selfconsistent field method is developed (spinorbit coupled stateaveraged complete active space selfconsistent field; CASSOC). The program follows a “onestep strategy” and treats the spinorbit interaction (SOC) on the same footing as the electronelectron interaction. As opposed to other existing approaches, the program employs an intermediate coupling scheme in which spin and space symmetry adapted configuration space functions are allowed to interact via SOC. This adds to the transparency and computational efficiency of the procedure. The approach requires the utilization of complexvalued configuration interaction coefficients, but the molecular orbital coefficients can be kept realvalued without loss of generality. Hence, expensive arithmetic associated with evaluation of complexvalued transformed molecular integrals is completely avoided. In order to investigate the quality of the calculated wave function, we extended the method to the calculation of electronic gtensors. As the SOC is already treated to all orders in the SACASSCF process, first order perturbation theory with the Zeeman operator is sufficient to accomplish this task. As a testset, we calculated gtensors of a set of diatomics, a set of d^{1} transition metal complexes MOX_{4} ^{n−}, and a set of 5f^{1} actinide complexes AnX_{6} ^{n−}. These calculations reveal that the effect of the wavefunction relaxation due to variation inclusion of SOC is of the same order of magnitude as the effect of inclusion of dynamic correlation and hence cannot be neglected for the accurate prediction of electronic gtensors.

Reduction of chemical reaction networks through delay distributions
View Description Hide DescriptionAccurate modelling and simulation of dynamic cellular events require two main ingredients: an adequate description of key chemical reactions and simulation of such chemical events in reasonable time spans. Quite logically, posing the right model is a crucial step for any endeavour in Computational Biology. However, more often than not, it is the associated computational costs which actually limit our capabilities of representing complex cellular behaviour. In this paper, we propose a methodology aimed at representing chains of chemical reactions by much simpler, reduced models. The abridgement is achieved by generation of modelspecific delay distribution functions, consecutively fed to a delay stochastic simulation algorithm. We show how such delay distributions can be analytically described whenever the system is solely composed of consecutive firstorder reactions, with or without additional “backward” bypass reactions, yielding an exact reduction. For models including other types of monomolecular reactions (constitutive synthesis, degradation, or “forward” bypass reactions), we discuss why one must adopt a numerical approach for its accurate stochastic representation, and propose two alternatives for this. In these cases, the accuracy depends on the respective numerical sample size. Our model reduction methodology yields significantly lower computational costs while retaining accuracy. Quite naturally, computational costs increase alongside network size and separation of time scales. Thus, we expect our model reduction methodologies to significantly decrease computational costs in these instances. We anticipate the use of delays in model reduction will greatly alleviate some of the current restrictions in simulating large sets of chemical reactions, largely applicable in pharmaceutical and biological research.

Perturbative wavepacket spawning procedure for nonadiabatic dynamics in diabatic representation
View Description Hide DescriptionI present a new formulation of wavepacket spawning procedure based on a second order perturbation theory expression for population transfer between different diabatic electronic states. The employed perturbation theory (PT) expansion is based on an assumption that diabatic states can be represented locally with their Taylor series up to quadratic terms in nuclear coordinates (local harmonic approximation). The corresponding local harmonic basis of vibrational states makes infinite summation over excited states in PT expressions possible, and thus, it provides a complete basis set expression for the population transfer. This allows me to detect when a finite basis set expansion employed in variational wave packet propagation does not adequately describe the interstate population transfer. Also, it suggests a rigorous criterion for basis set expansion (spawning). The proposed procedure is illustrated for the variational multiconfigurational Gaussian wave packet method applied to 1D and 2D model examples, and it also can be extended to direct onthefly dynamics with any Gaussian wave packet propagation method.

Spectroscopic properties of Ar_{ x }–Zn and Ar_{ x }–Ag^{+} (x = 1,2) van der Waals complexes
View Description Hide DescriptionPotential energy curves have been constructed using coupled cluster with singles, doubles, and perturbative triple excitations (CCSD(T)) in combination with allelectron and pseudopotentialbased multiply augmented correlation consistent basis sets [maugccpV(n + d)Z; m = singly, doubly, triply, n = D,T,Q,5]. The effect of basis set superposition error on the spectroscopic properties of Ar–Zn, Ar_{2}–Zn, Ar–Ag^{+}, and Ar_{2}–Ag^{+} van der Waals complexes was examined. The diffuse functions of the doubly and triply augmented basis sets have been constructed using the eventempered expansion. The a posteriori counterpoise scheme of Boys and Bernardi and its generalized variant by Valiron and Mayer has been utilized to correct for basis set superposition error (BSSE) in the calculated spectroscopic properties for diatomic and triatomic species. It is found that even at the extrapolated complete basis set limit for the energetic properties, the pseudopotentialbased calculations still suffer from significant BSSE effects unlike the allelectron basis sets. This indicates that the quality of the approximations used in the design of pseudopotentials could have major impact on a seemingly valenceexclusive effect like BSSE. We confirm the experimentally determined equilibrium internuclear distance (r _{e}), binding energy (D _{e}), harmonic vibrational frequency (ω_{e}), and C^{1}Π ← X ^{1}Σ transition energy for ArZn and also predict the spectroscopic properties for the lowlying excited states of linear Ar_{2}–Zn (X ^{1}Σ_{g}, ^{3}Π_{g}, ^{1}Π_{g}), Ar–Ag^{+} (X ^{1}Σ, ^{3}Σ, ^{3}Π, ^{3}Δ, ^{1}Σ, ^{1}Π, ^{1}Δ), and Ar_{2}–Ag^{+} (X ^{1}Σ_{g}, ^{3}Σ_{g}, ^{3}Π_{g}, ^{3}Δ_{g}, ^{1}Σ_{g}, ^{1}Π_{g}, ^{1}Δ_{g}) complexes, using the CCSD(T) and MRCISD + Q methods, to aid in their experimental characterizations.

Fractional diffusionreaction stochastic simulations
View Description Hide DescriptionA novel method is presented for the simulation of a discrete state space, continuous time Markov process subject to fractional diffusion. The method is based on LieTrotter operator splitting of the diffusion and reaction terms in the master equation. The diffusion term follows a multinomial distribution governed by a kernel that is the discretized solution of the fractional diffusion equation. The algorithm is validated and simulations are provided for the FisherKPP wavefront. It is shown that the wave speed is dictated by the order of the fractional derivative, where lower values result in a faster wave than in the case of classical diffusion. Since many physical processes deviate from classical diffusion, fractional diffusion methods are necessary for accurate simulations.

Marcus canonical integral for nonGaussian processes and its computation: Pathwise simulation and tauleaping algorithm
View Description Hide DescriptionThe stochastic integral ensuring the NewtonLeibnitz chain rule is essential in stochastic energetics. Marcus canonical integral has this property and can be understood as the WongZakai type smoothing limit when the driving process is nonGaussian. However, this important concept seems not wellknown for physicists. In this paper, we discuss Marcus integral for nonGaussian processes and its computation in the context of stochastic energetics. We give a comprehensive introduction to Marcus integral and compare three equivalent definitions in the literature. We introduce the exact pathwise simulation algorithm and give the error analysis. We show how to compute the thermodynamic quantities based on the pathwise simulation algorithm. We highlight the information hidden in the Marcus mapping, which plays the key role in determining thermodynamic quantities. We further propose the tauleaping algorithm, which advance the process with deterministic time steps when tauleaping condition is satisfied. The numerical experiments and its efficiency analysis show that it is very promising.

Order parameter free enhanced sampling of the vaporliquid transition using the generalized replica exchange method
View Description Hide DescriptionThe generalized Replica Exchange Method (gREM) is extended into the isobaricisothermal ensemble, and applied to simulate a vaporliquid phase transition in LennardJones fluids. Merging an optimally designed generalized ensemble sampling with replica exchange, gREM is particularly well suited for the effective simulation of firstorder phase transitions characterized by “backbending” in the statistical temperature. While the metastable and unstable states in the vicinity of the firstorder phase transition are masked by the enthalpy gap in temperature replica exchange method simulations, they are transformed into stable states through the parameterized effective sampling weights in gREM simulations, and join vapor and liquid phases with a succession of unimodal enthalpy distributions. The enhanced sampling across metastable and unstable states is achieved without the need to identify a “good” order parameter for biased sampling. We performed gREM simulations at various pressures below and near the critical pressure to examine the change in behavior of the vaporliquid phase transition at different pressures. We observed a crossover from the firstorder phase transition at low pressure, characterized by the backbending in the statistical temperature and the “kink” in the Gibbs free energy, to a continuous secondorder phase transition near the critical pressure. The controlling mechanisms of nucleation and continuous phase transition are evident and the coexistence properties and phase diagram are found in agreement with literature results.

Charge carrier transport in molecularly doped polycarbonate as a test case for the dipolar glass model
View Description Hide DescriptionWe present the results of Monte Carlo simulations of the charge carrier transport in a disordered molecular system containing spatial and energetic disorders using the dipolar glass model. Model parameters of the material were chosen to fit a typical polar organic photoconductor polycarbonate doped with 30% of aromatic hydrazone, whose transport properties are well documented in literature. Simulated carrier mobility demonstrates a usual PooleFrenkel field dependence and its slope is very close to the experimental value without using any adjustable parameter. At room temperature transients are universal with respect to the electric field and transport layer thickness. At the same time, carrier mobility does not depend on the layer thickness and transients develop a welldefined plateau where the current does not depend on time, thus demonstrating a nondispersive transport regime. Tails of the transients decay as power law with the exponent close to −2. This particular feature indicates that transients are close to the boundary between dispersive and nondispersive transport regimes. Shapes of the simulated transients are in very good agreement with the experimental ones. In summary, we provide a first verification of a selfconsistency of the dipolar glass transport model, where major transport parameters, extracted from the experimental transport data, are then used in the transport simulation, and the resulting mobility field dependence and transients are in very good agreement with the initial experimental data.