Quantum dynamics calculations using symmetrized, orthogonal Weyl-Heisenberg wavelets with a phase space truncation scheme. II. Construction and optimization

In this paper, we extend and elaborate upon a wavelet method first presented in a previous publication [B. Poirier, J. Theo. Comput. Chem. 2, 65 (2003)]. In particular, we focus on construction and optimization of the wavelet functions, from theoretical and numerical viewpoints, and also examine their localization properties. The wavelets used are modified Wilson-Daubechies wavelets, which in conjunction with a simple phase space truncation scheme, enable one to solve the multidimensional Schrodinger equation. This approach is ideally suited to rovibrational spectroscopy applications, but can be used in any context where differential equations are involved.     

Rovibrational spectroscopy calculations of neon dimer using a phase space truncated Weyl-Heisenberg wavelet basis

In a series of earlier articles [B. Poirier J. Theor. Comput. Chem. 2, 65 (2003); B. Poirier and A. Salam J. Chem. Phys. 121, 1690 (2004); B. Poirier and A. Salam J. Chem. Phys. 121, 1740 (2004)], a new method was introduced for performing exact quantum dynamics calculations in a manner that formally defeats exponential scaling with system dimensionality. The method combines an optimally localized, orthogonal Weyl-Heisenberg wavelet basis set with a simple phase space truncation scheme, and has already been applied to model systems up to 17 degrees of freedom (DOF's). In this paper, the approach is applied for the first time to a real molecular system (neon dimer), necessitating the development of an efficient numerical scheme for representing arbitrary potential energy functions in the wavelet representation. All bound rovibrational energy levels of neon dimer are computed, using both one DOF radial coordinate calculations and a three DOF Cartesian coordinate calculation. Even at such low dimensionalities, the approach is found to be competitive with another state-of-the-art method applied to the same system [J. Montgomery and B. Poirier J. Chem. Phys. 119, 6609 (2003)].     

Statistical mechanical theory for steady state systems. VI. Variational principles

Several variational principles that have been proposed for nonequilibrium systems are analyzed. These include the principle of minimum rate of entropy production due to Prigogine [Introduction to Thermodynamics of Irreversible Processes (Interscience, New York, 1967)], the principle of maximum rate of entropy production, which is common on the internet and in the natural sciences, two principles of minimum dissipation due to Onsager [Phys. Rev. 37, 405 (1931)] and to Onsager and Machlup [Phys. Rev. 91, 1505 (1953)], and the principle of maximum second entropy due to Attard [J. Chem.. Phys. 122, 154101 (2005); Phys. Chem. Chem. Phys. 8, 3585 (2006)]. The approaches of Onsager and Attard are argued to be the only viable theories. These two are related, although their physical interpretation and mathematical approximations differ. A numerical comparison with computer simulation results indicates that Attard's expression is the only accurate theory. The implications for the Langevin and other stochastic differential equations are discussed.     

Time-dependent wave packet study on trans-cis isomerization of HONO driven by an external field

The present paper is devoted to a full quantum mechanical study of the cis-->trans isomerization of HONO. In contrast to our previous study [Richter et al., J. Chem. Phys. 120, 6072 (2004)], the dynamics is now performed in the presence of an external time-dependent field in order to be closer to experimental conditions. A six-dimensional dipole surface is computed. Using a previously developed potential energy surface [Richter et al., J. Chem. Phys. 120, 1306 (2004)], all eigenstates up to 4000 cm(-1) are calculated. We simulate the dynamics during and after excitation by an electromagnetic pulse whose parameters are chosen to efficiently trigger the isomerization. Our investigations show that there is a selective isomerization pathway.     

Reconciling semiclassical and Bohmian mechanics. III. Scattering states for continuous potentials

In a previous paper [B. Poirier, J. Chem. Phys. 121, 4501 (2004)] a unique bipolar decomposition psi = psi1 + psi2 was presented for stationary bound states Psi of the one-dimensional Schrodinger equation, such that the components psi1 and psi2 approach their semiclassical WKB analogs in the large-action limit. The corresponding bipolar quantum trajectories, as defined in the usual Bohmian mechanical formulation, are classical-like and well behaved, even when Psi has many nodes or is wildly oscillatory. A modification for discontinuous potential stationary scattering states was presented in a second, companion paper [C. Trahan and B. Poirier, J. Chem. Phys.124, 034115 (2006), previous paper], whose generalization for continuous potentials is given here. The result is an exact quantum scattering methodology using classical trajectories. For additional convenience in handling the tunneling case, a constant-velocity-trajectory version is also developed.     

Definition of invariant molecular tensors in an origin change. III. Comparison of the methods suitable for the obtainment of invariant tensors for refraction and the Faraday effect

In a previous paper, denoted I [C. Duprez and J. Tillieu, Int. J. Quantum Chem. 17 , 155 (1980)], a global tensor, invariant in an origin change, has been introduced for all optical phenomena. The results obtained have been applied to two particular cases: refraction and the Faraday effect, in Paper II [C. Duprez, Int. J. Quantum Chem. 17 , 169 (1980)]. In this paper, we study the various methods suitable for the calculation of these tensors. Then, we compare the results obtained in II with those previously obtained by other authors who, when studying the two particular cases mentioned above, omitted some terms necessary to obtain the invariance.     

Non-Hamiltonian molecular dynamics implementation of the Gibbs ensemble method. I. Algorithm

In this paper (paper I) and the following paper (paper II) [C. Bratschi, H. Huber, and D. J. Searles, J. Chem. Phys. 126, 164105 (2007)], a new molecular dynamics algorithm implementing the Gibbs ensemble will be presented and then on applied to the liquid-vapor coexistence curve for two ab initio CO2 potentials. In paper I, the Gibbs ensemble molecular dynamics algorithm using non-Hamiltonian molecular dynamics techniques is introduced. It is shown that states of the correct probability density function are sampled and the correct exchange probability is generated. The extended system Nose-Hoover formalism is used to generate a constant temperature ensemble with equal pressures in the subsystems, combined with single particle transfers between the subsystems, over several time steps, to get equal chemical potentials.     

A practical treatment for the three-body interactions in the transcorrelated variational Monte Carlo method: application to atoms from lithium to neon

We suggest a practical solution to dealing with the three-body interactions in the transcorrelated variational Monte Carlo method (TC-VMC). In the TC-VMC method, which was suggested in our previous paper [N. Umezawa and S. Tsuneyuki, J. Chem. Phys. 119, 10015 (2003)], the Jastrow-Slater-type wave function is efficiently optimized through a self-consistent procedure by minimizing the variance of the local energy. The three-body terms in the transcorrelated self-consistent-field equation, which have been simply ignored in our previous works, are efficiently calculated by the Monte Carlo numerical integration. We found that our treatment for the three-body interactions is successful for atoms from Li to Ne.     

Linear scaling computation of the Fock matrix. VII. Periodic density functional theory at the Gamma point

Linear scaling quantum chemical methods for density functional theory are extended to the condensed phase at the Gamma point. For the two-electron Coulomb matrix, this is achieved with a tree-code algorithm for fast Coulomb summation [M. Challacombe and E. Schwegler, J. Chem. Phys. 106, 5526 (1997)], together with multipole representation of the crystal field [M. Challacombe, C. White, and M. Head-Gordon, J. Chem. Phys. 107, 10131 (1997)]. A periodic version of the hierarchical cubature algorithm [M. Challacombe, J. Chem. Phys. 113, 10037 (2000)], which builds a telescoping adaptive grid for numerical integration of the exchange-correlation matrix, is shown to be efficient when the problem is posed as integration over the unit cell. Commonalities between the Coulomb and exchange-correlation algorithms are discussed, with an emphasis on achieving linear scaling through the use of modern data structures. With these developments, convergence of the Gamma-point supercell approximation to the k-space integration limit is demonstrated for MgO and NaCl. Linear scaling construction of the Fockian and control of error is demonstrated for RBLYP6-21G* diamond up to 512 atoms.     

Complex autocorrelation function and energy spectrum by classical trajectory calculations

A quasiclassical method which enables evaluation of complex autocorrelation function from classical trajectory calculations is proposed. The method is applied for two highly excited nonlinearly coupled harmonic oscillators in regimes prevailed either by regular or chaotic classical motions. A good agreement of classical and quantum autocorrelation functions is found within short (Ehrnfest) time limit. Fourier transforms of the autocorrelation functions provide moderate resolved energy spectra, where classical and quantum results nearly coincide. The actual energy levels are obtained from approximate short-time autocorrelation functions with the help of filter diagonalization. This study is a follow up to our previous work [P. Zdanska and N. Moiseyev, J. Chem. Phys. 115, 10608 (2001)], where the complex autocorrelation has been obtained up to overall phase factors of recurrences.     

High precision variational calculations for the Born-Oppenheimer energies of the ground state of the hydrogen molecule

Born-Oppenheimer approximation Hylleraas variational calculations with up to 7034 expansion terms are reported for the 1sigma(g)+ ground state of neutral hydrogen at various internuclear distances. The nonrelativistic energy is calculated to be -1.174 475 714 220(1) hartree at R = 1.4 bohr, which is four orders of magnitude better than the best previous Hylleraas calculation, that of Wolniewicz [J. Chem. Phys. 103, 1792 (1995)]. This result agrees well with the best previous variational energy, -1.174 475 714 216 hartree, of Cencek (personal communication), obtained using explicitly correlated Gaussians (ECGs) [Cencek and Rychlewski, J. Chem. Phys. 98, 1252 (1993); Cencek et al., ibid. 95, 2572 (1995); Rychlewski, Adv. Quantum Chem. 31, 173 (1998)]. The uncertainty in our result is also discussed. The nonrelativistic energy is calculated to be -1.174 475 931 399(1) hartree at the equilibrium R = 1.4011 bohr distance. This result also agrees well with the best previous variational energy, -1.174 475 931 389 hartree, of Cencek and Rychlewski [Rychlewski, Handbook of Molecular Physics and Quantum Chemistry, edited by S. Wilson (Wiley, New York, 2003), Vol. 2, pp. 199-218; Rychlewski, Explicitly Correlated Wave Functions in Chemistry and Physics Theory and Applications, edited by J. Rychlewski (Kluwer Academic, Dordrecht, 2003), pp. 91-147.], obtained using ECGs.     

Local response dispersion method. II. Generalized multicenter interactions

Recently introduced local response dispersion method [T. Sato and H. Nakai, J. Chem. Phys. 131, 224104 (2009)], which is a first-principles alternative to empirical dispersion corrections in density functional theory, is implemented with generalized multicenter interactions involving both atomic and atomic pair polarizabilities. The generalization improves the asymptote of intermolecular interactions, reducing the mean absolute percentage error from about 30% to 6% in the molecular C(6) coefficients of more than 1000 dimers, compared to experimental values. The method is also applied to calculations of potential energy curves of molecules in the S22 database [P. Jure?ka et al., Phys. Chem. Chem. Phys. 8, 1985 (2006)]. The calculated potential energy curves are in a good agreement with reliable benchmarks recently published by Molnar et al. [J. Chem. Phys. 131, 065102 (2009)]. These improvements are achieved at the price of increasing complexity in the implementation, but without losing the computational efficiency of the previous two-center (atom-atom) formulation. A set of different truncations of two-center and three- or four-center interactions is shown to be optimal in the cost-performance balance.     

Problems in the comparison of theoretical and experimental hyperpolarizabilities revisited

The relationship between nonlinear susceptibilities and hyperpolarizabilities defined using different conventions is reexamined. In previous work [Willetts et al., J. Chem. Phys. 97, 7590 (1992)], relations between different conventions for microscopic hyperpolarizabilities have been derived, but the application of the corresponding conversion factors led to several inconsistencies. It is shown that different conventions for macroscopic susceptibilities have to be taken into account, too, in order to arrive at consistently comparable values. The complete set of conversion factors between several conventions are given for second harmonic generation, electric field induced second harmonic (EFISH) generation, and third harmonic generation. As an illustration, experimental EFISH and hyper-Rayleigh scattering results of p-nitroaniline are compared with each other and with recent results of ab initio computations including solvation effects. Several problems in the comparison of computational and experimental values are also discussed.     

Thermodynamic anomalies in a lattice model of water: solvation properties

We investigate a lattice-fluid model of water, defined on a three-dimensional body-centered-cubic lattice. Model molecules possess a tetrahedral symmetry, with four equivalent bonding arms. The model is similar to the one proposed by Roberts and Debenedetti [J. Chem. Phys. 105, 658 (1996)], simplified by removing distinction between "donors" and "acceptors." We focus on the solvation properties, mainly as far as an ideally inert (hydrophobic) solute is concerned. As in our previous analysis, devoted to neat water [J. Chem. Phys. 121, 11856 (2004)], we make use of a generalized first-order approximation on a tetrahedral cluster. We show that the model exhibits quite a coherent picture of water thermodynamics, reproducing qualitatively several anomalous properties observed both in pure water and in solutions of hydrophobic solutes. As far as supercooled liquid water is concerned, the model is consistent with the second critical-point scenario.     

Reconciling semiclassical and Bohmian mechanics. II. Scattering states for discontinuous potentials

In a previous paper [B. Poirier, J. Chem. Phys. 121, 4501 (2004)] a unique bipolar decomposition, psi = psi1 + psi2, was presented for stationary bound states Psi of the one-dimensional Schrodinger equation, such that the components psi1 and psi2 approach their semiclassical WKB analogs in the large action limit. Moreover, by applying the Madelung-Bohm ansatz to the components rather than to Psi itself, the resultant bipolar Bohmian mechanical formulation satisfies the correspondence principle. As a result, the bipolar quantum trajectories are classical-like and well behaved, even when psi has many nodes or is wildly oscillatory. In this paper, the previous decomposition scheme is modified in order to achieve the same desirable properties for stationary scattering states. Discontinuous potential systems are considered (hard wall, step potential, and square barrier/well), for which the bipolar quantum potential is found to be zero everywhere, except at the discontinuities. This approach leads to an exact numerical method for computing stationary scattering states of any desired boundary conditions, and reflection and transmission probabilities. The continuous potential case will be considered in a companion paper [C. Trahan and B. Poirier, J. Chem. Phys. 124, 034116 (2006), following paper].     

Molecular symmetry and ab initio calculations: IV. Symmetry-matrix and symmetry-supermatrix in calculations of two-electron repulsion integrals

The product of two Gaussians having different centers is itself a one-center Gaussian, thus multicenter integrals with a Cartesian Gaussian basis can be reduced to one-center integrals. Recurrence relations for overlap integrals and electron repulsion integrals (ERIs) are derived at these centers. The calculations of overlap integrals and ERIs are carried out step by step from the highest symmetry case (one center) to required cases (different centers) by using the translation of Cartesian Gaussians. Full exploitation of symmetry in calculation processes can result in optimal use of these recurrence relations. Compared with the recently published algorithms, based on the recurrence relations derived by Obara and Saika [J. Chem. Phys., 84 , 3963 (1986)], the floating point operations (FLOPs) for ERI calculations (having four different centers) can be reduced by a factor of ca. 2. A significant extra saving in calculations and storage can be obtained if atoms, linear, or planar molecules are discussed. © 1997 John Wiley & Sons, Inc.     

Vibrationally inelastic collisions in H+ +CO system: comparing quantum calculations with experiments

State-resolved cross beam experiments [H. Udseth et al., J. Chem. Phys. 60, 3051 (1974); J. Krutein and F. Linder, J. Chem. Phys. 71, 599 (1979); G. Niedner-Schatteburg and J. P. Toennies, Adv. Chem. Phys. LXXXII, 553 (1992)], coupled with proton energy loss spectroscopy for the inelastic scattering of H(+) from CO in the collision range of 10-30 eV show very low vibrational excitation of the target molecule. Stimulated by the experimentally observed low vibrational inelasticity in the system the ground and the first two low-lying excited electronic potential-energy surfaces have been computed using the ab initio multireference configuration interaction method. Quantum dynamics has been performed on the ground potential energy surface in the framework of vibrational close-coupling rotational infinite-order sudden approximation. The various computed dynamical attributes such as differential and integral cross sections, and average vibrational energy transfer are analyzed in detail, and compared successfully with the available experimental results.     

Multidimensional quantum trajectories: applications of the derivative propagation method

In a previous publication [J. Chem. Phys. 118, 9911 (2003)], the derivative propagation method (DPM) was introduced as a novel numerical scheme for solving the quantum hydrodynamic equations of motion (QHEM) and computing the time evolution of quantum mechanical wave packets. These equations are a set of coupled, nonlinear partial differential equations governing the time evolution of the real-valued functions C and S in the complex action, S=C(r,t) + iS(r,t)/Planck's over 2pi, where Psi(r,t)=exp(S). Past numerical solutions to the QHEM were obtained via ensemble trajectory propagation, where the required first- and second-order spatial derivatives were evaluated using fitting techniques such as moving least squares. In the DPM, however, equations of motion are developed for the derivatives themselves, and a truncated set of these are integrated along quantum trajectories concurrently with the original QHEM equations for C and S. Using the DPM quantum effects can be included at various orders of approximation; no spatial fitting is involved; there is no basis set expansion; and single, uncoupled quantum trajectories can be propagated (in parallel) rather than in correlated ensembles. In this study, the DPM is extended from previous one-dimensional (1D) results to calculate transmission probabilities for 2D and 3D wave packet evolution on coupled Eckart barrier/harmonic oscillator surfaces. In the 2D problem, the DPM results are compared to standard numerical integration of the time-dependent Schrodinger equation. Also in this study, the practicality of implementing the DPM for systems with many more degrees of freedom is discussed.     

Mixed Quantum Classical Reaction Rates based on the Phase Space Formulation of the Hierarchical Equations of Motion

In a previous work [J. Chem. Phys. 140 , 174105 (2014)], we have shown that a mixed quantum classical (MQC) rate theory can be derived to investigate the quantum tunneling effects in the proton transfer reactions. However, the method is based on the high temperature approximation of the hierarchical equation of motion (HEOM) with the Debye-Drude spectral density, and results in a multi-state Zusman type of equation. We now extend this theory to include quantum effects of the bath degrees of freedom. By writing the full HEOM into a multidimensional partial differential equation in phase space, we can define a new reaction coordinate, and the previous method can be generalized to the full quantum regime. The validity of the new method is demonstrated by using numerical examples, including the spin-Boson model, and the double well model for proton transfer reaction. The new method is found to resolve some key problems of the previous theory based on high temperature approximation, including possible numerical instability in long time simulation and wrong rate constant at low temperatures.     

A justification for a nonadiabatic surface hopping Herman-Kluk semiclassical initial value representation of the time evolution operator

A justification is given for the validity of a nonadiabatic surface hopping Herman-Kluk (HK) semiclassical initial value representation (SC-IVR) method. The method is based on a propagator that combines the single surface HK SC-IVR method [J. Chem. Phys. 84, 326 (1986)] and Herman's nonadiabatic semiclassical surface hopping theory [J. Chem. Phys. 103, 8081 (1995)], which was originally developed using the primitive semiclassical Van Vleck propagator. We show that the nonadiabatic HK SC-IVR propagator satisfies the time-dependent Schrodinger equation to the first order of variant Planck's over 2pi and the error is O(variant Planck's over 2pi(2)). As a required lemma, we show that the stationary phase approximation, under current assumptions, has an error term variant Planck's over 2pi(1) order higher than the leading term. Our derivation suggests some changes to the previous development, and it is shown that the numerical accuracy in applications to Tully's three model systems in low energies is improved.