Volume 130, Issue 12, 28 March 2009
 ARTICLES

 Theoretical Methods and Algorithms

Structures and harmonic vibrational frequencies for excited states of diatomic molecules with CCSD(R12) and CCSD(F12) models
View Description Hide DescriptionThe equationofmotion coupledcluster method for excited states with the singlesanddoubles model (CCSD) has been implemented for ansatz 2 of the explicitly correlated CCSD(R12) and CCSD(F12) methods as part of the program package Dalton. In this model, an orthonormal complementary auxiliary basis set is used for the resolutionofidentity approximation in order to calculate the threeelectron integrals needed for CCSD(R12) and CCSD(F12). The additional CCSD(R12) or CCSD(F12) terms introduced within ansatz 2, which are not present in ansatz 1, are derived and discussed with regard to the extra costs needed for their computation. As a first application the basis set convergence of equilibrium bond lengths and harmonic vibrational frequencies has been investigated for some singlet excited states of the diatomic molecules , CO, BF, and BH. The calculated CCSD(F12) results show that the average absolute deviations of the bond lengths and frequencies from the basis set limits are below 0.1 pm and as well as 0.05 pm and for the triple and quadruplebasis sets, respectively. These deviations are shown to largely arise from the SCF basis set incompleteness errors.

Quadratic canonical transformation theory and higher order density matrices
View Description Hide DescriptionCanonical transformation (CT) theory provides a rigorously sizeextensive description of dynamic correlation in multireference systems, with an accuracy superior to and cost scaling lower than complete active space second order perturbation theory. Here we expand our previous theory by investigating (i) a commutator approximation that is applied at quadratic, as opposed to linear, order in the effective Hamiltonian, and (ii) incorporation of the threebody reduced density matrix in the operator and density matrix decompositions. The quadratic commutator approximation improves CT’s accuracy when used with a singledeterminant reference, repairing the previous formal disadvantage of the singlereference linear CT theory relative to singles and doubles coupled clustertheory. Calculations on the BH and HF binding curves confirm this improvement. In multireference systems, the threebody reduced density matrix increases the overall accuracy of the CT theory. Tests on the and binding curves yield results highly competitive with expensive stateoftheart multireference methods, such as the multireference Davidsoncorrected configuration interaction, averaged coupled pair functional, and averaged quadratic coupled clustertheories.

Relativistic allelectron molecular dynamics simulations
View Description Hide DescriptionThe scalarrelativistic Douglas–Kroll–Hess method is implemented in the Born–Oppenheimer molecular dynamics simulation package CP2K. Using relativistic densities in a nonrelativistic gradient routine is found to be a valid approximation of relativistic gradients. An excellent agreement between optimized structures and geometries obtained from numerical gradients is observed with an error smaller than 0.02 pm. Hydrogen halide dimers [, with , Cl, Br, I] serve as small test systems for firstprinciplesmolecular dynamics simulations.Relativistic effects are observed. That is, the amplitude of motion is larger, the frequency of motion is smaller, and the distances are larger in the relativistic picture. Several localization schemes are evaluated for different interatomic and intermolecular distances. The errors of these localization schemes are small for geometries which are similar to the equilibrium structure. They become larger for smaller distances, introducing a slight bias toward closed packed configurations.

Multivariate frequency domain analysis of protein dynamics
View Description Hide DescriptionMultivariate frequency domain analysis (MFDA) is proposed to characterize collective vibrational dynamics of protein obtained by a molecular dynamics (MD) simulation. MFDA performs principal component analysis (PCA) for a bandpass filtered multivariate time series using the multitaper method of spectral estimation. By applying MFDA to MD trajectories of bovine pancreatic trypsin inhibitor, we determined the collective vibrational modes in the frequency domain, which were identified by their vibrational frequencies and eigenvectors. At near zero temperature, the vibrational modes determined by MFDA agreed well with those calculated by normal modeanalysis. At 300 K, the vibrational modes exhibited characteristic features that were considerably different from the principal modes of the static distribution given by the standard PCA. The influences of aqueous environments were discussed based on two different sets of vibrational modes, one derived from a MD simulation in water and the other from a simulation in vacuum. Using the varimax rotation, an algorithm of the multivariate statistical analysis, the representative orthogonal set of eigenmodes was determined at each vibrational frequency.

An efficient algorithm for the densityfunctional theory treatment of dispersion interactions
View Description Hide DescriptionThe quasiselfconsistentfield dispersioncorrected densityfunctional theory formalism (QSCFDCDFT) is developed and presented as an efficient and reliable scheme for the DFT treatment of van der Waals dispersion complexes, including full geometry optimizations and frequency calculations with analytical energy derivatives in a routine way. For this purpose, the longrangecorrected Perdew–Burke–Ernzerhof exchange functional and the oneparameter progressive correlation functional of Hirao and coworkers are combined with the Andersson–Langreth–Lundqvist (ALL) longrange correlation functional. The timeconsuming selfconsistent incorporation of the ALL term in the DFT iterations needed for the calculation of forces and force constants is avoided by an a posteriori evaluation of the ALL term and its gradient based on an effective partitioning of the coordinate space into global and intramonomer coordinates. QSCFDCDFT is substantially faster than SCFDCDFT would be. QSCFDCDFT is used to explore the potential energy surface (PES) of the benzene dimer. The results for the binding energies and intermolecular distances agree well with coupledcluster calculations at the complete basisset limit. We identify 16 stationary points on the PES, which underlines the usefulness of analytical energy gradients for the investigation of the PES. Furthermore, the inclusion of analytically calculated zero point energies reveals that largeamplitude vibrations connect the eight most stable benzene dimer forms and make it difficult to identify a dominating complex form. The tilted T structure and the paralleldisplaced sandwich form have the same value of 2.40 kcal/mol, which agrees perfectly with the experimental value of .

Calculation of long time classical trajectories: Algorithmic treatment and applications for molecular systems
View Description Hide DescriptionWe study the problem of finding a path that joins a given initial state with a final one, where the evolution is governed by classical (Hamiltonian) dynamics. A new algorithm for the computation of long time transition trajectories connecting two configurations is presented. In particular, a strategy for finding transition paths between two stable basins is established. The starting point is the formulation of the equation of motion of classical mechanics in the framework of Jacobi’s principle; a shortening procedure inspired by Birkhoff’s method is then applied to find geodesic solutions. Numerical examples are given for Müller’s potential and the collinear reaction.

Nonadiabatic coupling vectors within linear response timedependent density functional theory
View Description Hide DescriptionA method is developed to compute the nonadiabatic coupling vectors (NACVs) between electronic ground and excited states as well as between any possible pair of excited states within the framework of linear response timedependent density functional theory (TDDFT) in the adiabatic approximation. The development is an extension to our previous work on surface hopping dynamics [E. Tapavicza et al., Phys. Rev. Lett.98, 023001 (2007)] for which we improve the description of the TDDFT approximation of the excited statewavefunctions by means of linear response orbitals. The method is first validated on the system that has a region of strong coupling near the conical intersection at the equilateral geometry. These results confirm the quality and the numerical efficiency of the approach, which has an accuracy comparable to the one achieved with wavefunctionbased methods. Finally, we apply the method to the calculation of the NACVs of protonated formaldimine along a surface hopping trajectory initiated in the second excited state.

A Chebychev propagator for inhomogeneous Schrödinger equations
View Description Hide DescriptionA propagation scheme for timedependent inhomogeneous Schrödinger equations is presented. Such equations occur in time dependent optimal control theory and in reactive scattering. A formal solution based on a polynomial expansion of the inhomogeneous term is derived. It is subjected to an approximation in terms of Chebychev polynomials. Different variants for the inhomogeneous propagator are demonstrated and applied to two examples from optimal control theory. Convergence behavior and numerical efficiency are analyzed.

Highly accurate tauleaping methods with random corrections
View Description Hide DescriptionWe aim to construct higher order tauleaping methods for numerically simulating stochastic chemical kineticsystems in this paper. By adding a random correction to the primitive tauleaping scheme in each time step, we greatly improve the accuracy of the tauleaping approximations. This gain in accuracy actually comes from the reduction in the local truncation error of the scheme in the order of , the marching time step size. While the local truncation error of the primitive tauleaping method is for all moments, our Poisson random correction tauleaping method, in which the correction term is a Poisson random variable, can reduce the local truncation error for the mean to , and both Gaussian random correction tauleaping methods, in which the correction term is a Gaussian random variable, can reduce the local truncation error for both the mean and covariance to . Numerical results demonstrate that these novel methods more accurately capture crucial properties such as the mean and variance than existing methods for simulating chemical reactionsystems. This work constitutes a first step to construct high order numerical methods for simulating jump processes. With further refinement and appropriately modified stepsize selection procedures, the random correction methods should provide a viable way of simulating chemical reactionsystems accurately and efficiently.

Electrostatic correlations in colloidal suspensions: Density profiles and effective charges beyond the Poisson–Boltzmann theory
View Description Hide DescriptionA theory is proposed which allows us to calculate the distribution of the multivalent counterions around a colloidal particle using the cell model. The results are compared with the Monte Carlo simulations and are found to be very accurate in the two asymptotic regimes, close to the colloidal particle and far from it. The theory allows to accurately calculate the osmotic pressure and the effective charge of colloidal particles with multivalent counterions.

Comparison between integrated and parallel tempering methods in enhanced sampling simulations
View Description Hide DescriptionRecently, we introduced an integrated tempering approach to enhance sampling in the energy and configuration space for large systems. In this paper, we show that this new method has a higher efficiency than bias potential and generalized ensemble methods, such as accelerated molecular dynamics and replicaexchange molecular dynamics (parallel tempering) methods, in yielding thermodynamic averages. Particularly, the sampling efficiencies in both energy and configuration spaces are compared in details between integrated and parallel tempering methods. Related issues regarding the efficiency involved in the usage of the parallel tempering method are also discussed.

Replica exchange statistical temperature Monte Carlo
View Description Hide DescriptionThe replica exchange statistical temperature Monte Carlo algorithm (RESTMC) is presented, extending the singlereplica STMC algorithm [J. Kim, J. E. Straub, and T. Keyes, Phys. Rev. Lett.97, 050601 (2006)] to alleviate the slow convergence of the conventional temperature replica exchange method (REM) with increasing system size. In contrast to the Gibbs–Boltzmann sampling at a specific temperature characteristic of the standard REM, RESTMC samples a range of temperatures in each replica and achieves a flat energy sampling employing the generalized sampling weight, which is automatically determined via the dynamic modification of the replicadependent statistical temperature. Faster weight determination, through the dynamic update of the statistical temperature, and the flat energy sampling, maximizing energy overlaps between neighboring replicas, lead to a considerable acceleration in the convergence of simulations even while employing significantly fewer replicas. The performance of RESTMC is demonstrated and quantitatively compared with that of the conventional REM under varying simulation conditions for LennardJones 19, 31, and 55 atomic clusters, exhibiting single and doublefunneled energy landscapes.

A new electronic structure method for doublet states: Configuration interaction in the space of ionized and determinants
View Description Hide DescriptionAn implementation of gradient and energy calculations for configuration interaction variant of equationofmotion coupled cluster with single and double substitutions for ionization potentials (EOMIPCCSD) is reported. The method (termed IPCISD) treats the ground and excited doublet electronic states of an electron system as ionizing excitations from a closedshell electron reference state. The method is naturally spin adapted, variational, and size intensive. The computational scaling is , in contrast with the scaling of EOMIPCCSD. The performance and capabilities of the new approach are demonstrated by application to the uracil cation and water and benzene dimer cations by benchmarking IPCISD against more accurate IPCCSD. The equilibrium geometries, especially relative differences between different ionized states, are well reproduced. The average absolute errors and the standard deviations averaged for all bond lengths in all electronic states (58 values in total) are 0.014 and 0.007 Å, respectively. IPCISD systematically underestimates intramolecular distances and overestimates intermolecular ones, because of the underlying uncorrelated Hartree–Fock reference wave function. The IPCISD excitation energies of the cations are of a semiquantitative value only, showing maximum errors of 0.35 eV relative to EOMIPCCSD. Trends in properties such as dipole moments, transition dipoles, and charge distributions are well reproduced by IPCISD.

Numerical evaluation of electron repulsion integrals for pseudoatomic orbitals and their derivatives
View Description Hide DescriptionA numerical method to calculate the fourcenter electronrepulsion integrals for strictly localized pseudoatomic orbital basis sets has been developed. Compared to the conventional Gaussian expansion method, this method has an advantage in the ease of combination with density functional calculations. Additional mathematical derivations are also presented including the analytic derivatives of the integrals with respect to atomic positions and spatial damping of the Coulomb interaction due to the screening effect. In the numerical test for a simple molecule, the convergence up to in energy is successfully obtained with a feasible cost of computation.

The staticexchange electronwater pseudopotential, in conjunction with a polarizable water model: A new Hamiltonian for hydratedelectron simulations
View Description Hide DescriptionPreviously, Turi and Borgis [J. Chem. Phys.117, 6186 (2002)] parametrized an electronwater interaction potential, intended for use in simulations of hydrated electrons, by considering in the “static exchange” (essentially, frozencore Hartree–Fock) approximation, then applying an approximate Phillips–Kleinman procedure to construct a oneelectron pseudopotential representing the electronwater interaction. To date, this pseudopotential has been used exclusively in conjunction with a simple point charge water model that is parametrized for bulk water and yields poor results for small, neutral water clusters. Here, we extend upon the work of Turi and Borgis by reparametrizing the electronwater pseudopotential for use with the AMOEBA water model, which performs well for neutral clusters. The result is a oneelectron model Hamiltonian for , in which the oneelectron wave function polarizes the water molecules, and vice versa, in a fully selfconsistent fashion. The new model is fully variational and analytic energy gradients are available. We have implemented the new model using a modified Davidson algorithm to compute eigenstates, with the unpaired electron represented on a realspace grid. Comparison to ab initio electronic structure calculations for cluster isomers ranging from to reveals that the new model is significantly more accurate than the Turi–Borgis model, for both relative isomer energies and for vertical electron detachment energies. Electronwater polarizationinteractions are found to be much more significant for cavity states of the unpaired electron than for surface states.

Bifurcation of noreturn transition states in manybody chemical reactions
View Description Hide DescriptionA new method is presented to study bifurcation of noreturn transition states (TSs) at potential saddles for systems of many degrees of freedom (dof). The method enables us to investigate analytically when and how the noreturn TS bifurcates. Our method reveals a new aspect of bifurcation for systems of many dof, i.e., the action variables of the bath dof play a role of control parameters as long as they remain approximately conserved. As an illustrative example, we demonstrate our new method by using a three atomic exchange reaction. The bifurcation of noreturn TSs gives rise to a shortlived intermediate state at the saddle, which results in the overestimation of the reaction rate. Hence, the understanding of the bifurcation of the noreturn TS is crucial to capture the complexity in kinetics and dynamics of the reactions. The definability of noreturn TSs in manybody chemical reactions is also addressed under the occurrence of bifurcation above the reaction threshold.

Phasespace surface hopping: Nonadiabatic dynamics in a superadiabatic basis
View Description Hide DescriptionIn this paper, we construct a phasespace surface hopping algorithm for use in systems that exhibit strong nonadiabatic coupling. The algorithm is derived from a representation of the electronic basis which is a function of the nuclear phasespace coordinates rather than the nuclear position coordinates. This phasespace adiabatic basis can be understood in the context of Berry’s superadiabatic basis formalism as the firstorder superadiabatic correction to the conventional positionspace adiabatic basis. This superadiabatic representation leads to nuclear dynamics described not by Newton’s equations of motion but by generalized Hamilton’s equations of motion. The phasespace surface hopping algorithm captures physical effects that cannot be described by traditional algorithms. For a simple model problem, we show that phasespace surface hopping is more accurate than positionspace surface hopping, especially when the nonadiabatic coupling is strong.

Coupled cluster and density functional studies on geometries and energies of excited states of ozone
View Description Hide DescriptionThe performance of singledeterminant methods for finding geometries and energies of excited states is tested on the ozone molecule. Geometries for lowlying singlet and triplet states of ozone were optimized by CCSD(T) and density functional theory(DFT) (with BPW91 functional) methods. DFT geometries were found to lie close to CCSD(T) values. Most CCSD(T) and DFT geometries and energies are in good agreement with available experimental and recent highlevel theoretical values, with deviations lying within 0.02 Å, 2°, and 0.3 eV. An exception is the state, having a larger deviation of bond distance and energy. A multiconfigurational treatment is required for this state. DFT geometry optimizations and calculations of vibrational frequencies were extended to higher states, covering over 30 excited states of ozone, with adiabatic excitation energies up to about 6 eV. Calculated harmonic frequencies showed several states, including , to be saddle points. Multireference configuration interaction (MRCI) bending potentials for first and second singlet and triplet states were used in verifying the CCSD(T) and DFT geometries and for locating additional minima. For first states, DFT bending potentials are compared with MRCI potentials. As a criterion for the quality of singledeterminant geometries and energies of excited states, comparison of their vertical excitation energies with MRCI or timedependent DFT values is recommended.

Rate constants calculation with a simple mixed quantum/classical implementation of the fluxflux correlation function method
View Description Hide DescriptionA simple mixed quantum/classical (mixedQ/C) implementation of the fluxflux correlation function method has been applied to evaluate rate constants for a twodimensional model system. The model consists of an Eckart barrier resembling the collinear reaction, linearly coupled to a harmonic oscillator. Results are presented for a broad range of parameters for temperatures between 140 and 300 K. It is found that the mixedQ/C method gives fairly accurate results as long as the reaction does not involve too many recrossings. This suggests that the methodology could be extended to treat direct polyatomic reactions in gas phase.

Charge asymmetry in pure vibrational states of the HD molecule
View Description Hide DescriptionVery accurate variational calculations of all rotationless states (also called pure vibrational states) of the HD molecule have been performed within the framework that does not assume the Born–Oppenheimer (BO) approximation. The nonBO wave functions of the states describing the internal motion of the proton, the deuteron, and the two electrons were expanded in terms of onecenter explicitly correlated Gaussian functions multiplied by even powers of the internuclear distance. Up to 6000 functions were used for each state. Both linear and nonlinear parameters of the wave functions of all 18 states were optimized with a procedure that employs the analytical gradient of the energy with respect to the nonlinear parameters of the Gaussians. These wave functions were used to calculate expectation values of the interparticle distances and some other related quantities. The results allow elucidation of the charge asymmetry in HD as a function of the vibrational excitation.