Ab initio molecular dynamics with discrete variable representation basis sets: techniques and application to liquid water

Finite temperature ab initio molecular dynamics (AIMD), in which forces are obtained from "on-the-fly" electronic structure calculations, is a widely used technique for studying structural and dynamical properties of chemically active systems. Recently, we introduced an AIMD scheme based on discrete variable representation (DVR) basis sets, which was shown to have improved convergence properties over the conventional plane wave (PW) basis set [Liu,Y.; et al. Phys. Rev. B 2003, 68, 125110]. In the present work, the numerical algorithms for the DVR based AIMD scheme (DVR/AIMD) are provided in detail, and the latest developments of the approach are presented. The accuracy and stability of the current implementation of the DVR/AIMD scheme are tested by performing a simulation of liquid water at ambient conditions. The structural information obtained from the present work is in good agreement with the result of recent AIMD simulations with a PW basis set (PW/AIMD). Advantages of using the DVR/AIMD scheme over the PW/AIMD method are discussed. In particular, it is shown that a DVR/AIMD simulation of liquid water in the complete basis set limit is possible with a relatively small number of grid points.     

Quantum-classical dynamics of scattering processes in adiabatic and diabatic representations

We demonstrate the workability of a TDDVR based [J. Chem. Phys. 118, 5302 (2003)], novel quantum-classical approach, for simulating scattering processes on a quasi-Jahn-Teller model [J. Chem. Phys. 105, 9141 (1996)] surface. The formulation introduces a set of DVR grid points defined by the Hermite part of the basis set in each dimension and allows the movement of grid points around the central trajectory. With enough trajectories (grid points), the method converges to the exact quantum formulation whereas with only one grid point, we recover the conventional molecular dynamics approach. The time-dependent Schrodinger equation and classical equations of motion are solved self-consistently and electronic transitions are allowed anywhere in the configuration space among any number of coupled states. Quantum-classical calculations are performed on diabatic surfaces (two and three) to reveal the effects of symmetry on inelastic and reactive state-to-state transition probabilities, along with calculations on an adiabatic surface with ordinary Born-Oppenheimer approximation. Excellent agreement between TDDVR and DVR results is obtained in both the representations.     

Optimal grids for generalized finite basis and discrete variable representations: definition and method of calculation

The method of optimal generalized finite basis and discrete variable representations (FBR and DVR) generalizes the standard, Gaussian quadrature grid-classical orthonormal polynomial basis-based FBR/DVR method to general sets of grid points and to general, nondirect product, and/or nonpolynomial bases. Here, it is shown how an optimal set of grid points can be obtained for an optimal generalized FBR/DVR calculation with a given truncated basis. Basis set optimized and potential optimized grids are defined. The optimized grids are shown to minimize a function of grid points derived by relating the optimal generalized FBR of a Hamiltonian operator to a non-Hermitian effective Hamiltonian matrix. Locating the global minimum of this function can be reduced to finding the zeros of a function in the case of one dimensional problems and to solving a system of D nonlinear equations repeatedly in the case of D>1 dimensional problems when there is an equal number of grid points and basis functions. Gaussian quadrature grids are shown to be basis optimized grids. It is demonstrated by a numerical example that an optimal generalized FBR/DVR calculation of the eigenvalues of a Hamiltonian operator with potential optimized grids can have orders of magnitude higher accuracy than a variational calculation employing the same truncated basis. Nevertheless, for numerical integration with the optimal generalized FBR quadrature rule basis optimized grids are the best among grids of the same number of points. The notions of Gaussian quadrature and Gaussian quadrature accuracy are extended to general, multivariable basis functions.     

How to choose one-dimensional basis functions so that a very efficient multidimensional basis may be extracted from a direct product of the one-dimensional functions: energy levels of coupled systems with as many as 16 coordinates

In this paper we propose a scheme for choosing basis functions for quantum dynamics calculations. Direct product bases are frequently used. The number of direct product functions required to converge a spectrum, compute a rate constant, etc., is so large that direct product calculations are impossible for molecules or reacting systems with more than four atoms. It is common to extract a smaller working basis from a huge direct product basis by removing some of the product functions. We advocate a build and prune strategy of this type. The one-dimensional (1D) functions from which we build the direct product basis are chosen to satisfy two conditions: (1) they nearly diagonalize the full Hamiltonian matrix; (2) they minimize off-diagonal matrix elements that couple basis functions with diagonal elements close to those of the energy levels we wish to compute. By imposing these conditions we increase the number of product functions that can be removed from the multidimensional basis without degrading the accuracy of computed energy levels. Two basic types of 1D basis functions are in common use: eigenfunctions of 1D Hamiltonians and discrete variable representation (DVR) functions. Both have advantages and disadvantages. The 1D functions we propose are intermediate between the 1D eigenfunction functions and the DVR functions. If the coupling is very weak, they are very nearly 1D eigenfunction functions. As the strength of the coupling is increased they resemble more closely DVR functions. We assess the usefulness of our basis by applying it to model 6D, 8D, and 16D Hamiltonians with various coupling strengths. We find approximately linear scaling.     

Basis set construction for molecular electronic structure theory: natural orbital and Gauss-Slater basis for smooth pseudopotentials

A simple yet general method for constructing basis sets for molecular electronic structure calculations is presented. These basis sets consist of atomic natural orbitals from a multiconfigurational self-consistent field calculation supplemented with primitive functions, chosen such that the asymptotics are appropriate for the potential of the system. Primitives are optimized for the homonuclear diatomic molecule to produce a balanced basis set. Two general features that facilitate this basis construction are demonstrated. First, weak coupling exists between the optimal exponents of primitives with different angular momenta. Second, the optimal primitive exponents for a chosen system depend weakly on the particular level of theory employed for optimization. The explicit case considered here is a basis set appropriate for the Burkatzki-Filippi-Dolg pseudopotentials. Since these pseudopotentials are finite at nuclei and have a Coulomb tail, the recently proposed Gauss-Slater functions are the appropriate primitives. Double- and triple-zeta bases are developed for elements hydrogen through argon. These new bases offer significant gains over the corresponding Burkatzki-Filippi-Dolg bases at various levels of theory. Using a Gaussian expansion of the basis functions, these bases can be employed in any electronic structure method. Quantum Monte Carlo provides an added benefit: expansions are unnecessary since the integrals are evaluated numerically.     

A multidimensional discrete variable representation basis obtained by simultaneous diagonalization

Direct product basis functions are frequently used in quantum dynamics calculations, but they are poor in the sense that many such functions are required to converge a spectrum, compute a rate constant, etc. Much better, contracted, basis functions, that account for coupling between coordinates, can be obtained by diagonalizing reduced dimension Hamiltonians. If a direct product basis is used, it is advantageous to use discrete variable representation (DVR) basis functions because matrix representations of functions of coordinates are diagonal in the DVR. By diagonalizing matrices representing coordinates it is straightforward to obtain the DVR that corresponds to any direct product basis. Because contracted basis functions are eigenfunctions of reduced dimension Hamiltonians that include coupling terms they are not direct product functions. The advantages of contracted basis functions and the advantages of the DVR therefore appear to be mutually exclusive. A DVR that corresponds to contracted functions is unknown. In this paper we propose such a DVR. It spans the same space as a contracted basis, but in it matrix representations of coordinates are diagonal. The DVR basis functions are chosen to achieve maximal diagonality of coordinate matrices. We assess the accuracy of this DVR by applying it to model four-dimensional problems.     

All‐electron LCAO calculations of the LiF crystal phonon spectrum: Influence of the basis set,the exchange‐correlation functional,and the supercell size

For the first time the convergence of the phonon frequencies and dispersion curves in terms of the supercell size is studied in ab initio frozen phonon calculations on LiF crystal. Helmann–Feynman forces over atomic displacements are found in all‐electron calculations with the localized atomic functions (LCAO) basis using CRYSTAL06 program. The Parlinski–Li–Kawazoe method and FROPHO program are used to calculate the dynamical matrix and phonon frequencies of the supercells. For fcc lattice, it is demonstrated that use of the full supercell space group (including the supercell inner translations) enables to reduce essentially the number of the displacements under consideration. For Hartree–Fock (HF), PBE and hybrid PBE0, B3LYP, and B3PW exchange‐correlation functionals the atomic basis set optimization is performed. The supercells up to 216 atoms (3 × 3 × 3 conventional unit cells) are considered. The phonon frequencies using the supercells of different size and shape are compared. For the commensurate with supercell k ‐points the best agreement of the theoretical results with the experimental data is found for B3PW exchange‐correlation functional calculations with the optimized basis set. The phonon frequencies at the most non‐commensurate k ‐points converged for the supercell consisting of 4 × 4 × 4 primitive cells and ensures the accuracy 1–2% in the thermodynamic properties calculated (the Helmholtz free energy, entropy, and heat capacity at the room temperature). © 2009 Wiley Periodicals, Inc. J Comput Chem 2009     

Structure of liquid water at ambient temperature from ab initio molecular dynamics performed in the complete basis set limit

Structural properties of liquid water at ambient temperature were studied using Car-Parrinello [Phys. Rev. Lett. 55, 2471 (1985)] ab initio molecular dynamics (CPAIMD) simulations combined with the Kohn-Sham (KS) density functional theory and the BLYP exchange-correlation functional for the electronic structure. Unlike other recent work on the same subject, where plane-wave (PW) or hybrid Gaussian/plane-wave basis sets were employed, in the present paper, a discrete variable representation (DVR) basis set is used to expand the KS orbitals, so that with the real-space grid adapted in the present work, the properties of liquid water could be obtained very near the complete basis set limit. Structural properties of liquid water were extracted from a 30 ps CPAIMD-BLYP/DVR trajectory at 300 K. The radial distribution functions (RDFs), spatial distribution functions, and hydrogen bond geometry obtained from the CPAIMD-BLYP/DVR simulation are generally in good agreement with the most up to date experimental measurements. Compared to recent ab initio MD simulations based on PW basis sets, less significant overstructuring was found in the RDFs and the distributions of hydrogen bond angles, suggesting that previous plane-wave and Gaussian basis set calculations have exaggerated the tendency toward overstructuring.     

Cubic-grid Gaussian basis sets for electron scattering calculations. I. Definition and construction

We propose a new type of Gaussian basis sets for use in calculations of electron scattering by molecules. Instead of locating the basis-set functions on the atomic centers of the target molecule, we place primitive s-type Gaussians at the positions of a cubic lattice with a regular grid. The grid and the Gaussian exponent are fixed so as to give the best representation of the plane-wave function. Plane-wave functions and Green functions obtained by means of the cubic-grid basis set are tested graphically against exact functions and functions expressed by means of a conventional Gaussian basis set. © 1995 John Wiley & Sons, Inc.     

Degeneracy in discrete variable representations: general considerations and application to the multiconfigurational time-dependent Hartree approach

Problems appear in discrete variable representations (DVRs) based on general basis sets when the coordinate matrix has degenerate eigenvalues. Then the DVR is not uniquely defined. This paper shows that this problem can be caused by symmetry. Taking the symmetry into account when constructing the DVR solves the problem. The symmetry effect can be particularly important for the time-dependent DVR used in multiconfigurational time-dependent Hartree calculations employing the correlation DVR (CDVR) approach. Problems reported previously for the initial-state selected treatment of the H+H(2) reaction can be attributed to this symmetry effect. They can be solved by using a symmetry-adapted approach to construct the time-dependent DVR. Thus, the present paper shows that the CDVR scheme can be employed also in initial-state selected scattering calculations if the symmetry of the system is properly taken into account in the construction of the time-dependent DVR.     

Variational properties of the discrete variable representation: Discrete variable representation via effective operators

A variational finite basis representation/discrete variable representation (FBR/DVR) Hamiltonian operator has been introduced. By calculating its matrix elements exactly one obtains, depending on the choice of the basis set, either a variational FBR or a variational DVR. The domain of grid points on which the FBR/DVR is variational has been shown to consist of the subsets of the set of grid points one obtains by diagonalizing commuting variational basis representations of the coordinate operators. The variational property implies that the optimal of the subsets of a fixed number of points, i.e., the subset which gives the possible highest accuracy eigenpairs, gives the DVR of the smallest trace. The symmetry properties of the variational FBR/DVR Hamiltonian operator are analyzed and methods to incorporate symmetry into FBR/DVR calculations are discussed. It is shown how the Fourier-basis FBR/DVR suitable to solving periodic systems arise within the theory presented. Numerical examples are given to illustrate the theoretical results. The use of variational effective Hamiltonian and coordinate operators has been instrumental in this study. They have been introduced in a novel way by exploiting quasi-Hermiticity.     

Treating singularities present in the Sutcliffe-Tennyson vibrational Hamiltonian in orthogonal internal coordinates

Two methods are developed, when solving the related time-independent Schrodinger equation (TISE), to cope with the singular terms of the vibrational kinetic energy operator of a triatomic molecule given in orthogonal internal coordinates. The first method provides a mathematically correct treatment of all singular terms. The vibrational eigenfunctions are approximated by linear combinations of functions of a three-dimensional nondirect-product basis, where basis functions are formed by coupling Bessel-DVR functions, where DVR stands for discrete variable representation, depending on distance-type coordinates and Legendre polynomials depending on angle bending. In the second method one of the singular terms related to a distance-type coordinate, deemed to be unimportant for spectroscopic applications, is given no special treatment. Here the basis set is obtained by taking the direct product of a one-dimensional DVR basis with a two-dimensional nondirect-product basis, the latter formed by coupling Bessel-DVR functions and Legendre polynomials. With the basis functions defined, matrix representations of the TISE are set up and solved numerically to obtain the vibrational energy levels of H3+. The numerical calculations show that the first method treating all singularities is computationally inefficient, while the second method treating properly only the singularities having physical importance is quite efficient.     

Constant pressure ab initio molecular dynamics with discrete variable representation basis sets

The use of discrete variable representation (DVR) basis sets within ab initio molecular dynamics calculations allows the latter to be performed with converged energies and, more importantly, converged forces. In this paper, we show how to carry out ab initio molecular dynamics calculations in the isothermal-isobaric ensemble with fully flexible simulation boxes within the DVR basis set framework. In particular, we derive the appropriate DVR based expression for the pressure tensor when the electronic structure is represented using Kohn-Sham density functional theory, and we examine the convergence of this expression as a function of the basis set size. An illustrative example using 64 silicon atoms in a fully flexible box using a combination of the Martyna-Tobias-Klein [Martyna et al., J. Chem. Phys. 101, 4177 (1994)] and Car-Parrinello [Car and Parinello, Phys. Rev. Lett. 55, 2471 (1985)] algorithms is presented to demonstrate the efficacy of the approach.     

Finite basis representations with nondirect product basis functions having structure similar to that of spherical harmonics

The currently most efficient finite basis representation (FBR) method [Corey et al., in Numerical Grid Methods and Their Applications to Schrodinger Equation, NATO ASI Series C, edited by C. Cerjan (Kluwer Academic, New York, 1993), Vol. 412, p. 1; Bramley et al., J. Chem. Phys. 100, 6175 (1994)] designed specifically to deal with nondirect product bases of structures phi(n) (l)(s)f(l)(u), chi(m) (l)(t)phi(n) (l)(s)f(l)(u), etc., employs very special l-independent grids and results in a symmetric FBR. While highly efficient, this method is not general enough. For instance, it cannot deal with nondirect product bases of the above structure efficiently if the functions phi(n) (l)(s) [and/or chi(m) (l)(t)] are discrete variable representation (DVR) functions of the infinite type. The optimal-generalized FBR(DVR) method [V. Szalay, J. Chem. Phys. 105, 6940 (1996)] is designed to deal with general, i.e., direct and/or nondirect product, bases and grids. This robust method, however, is too general, and its direct application can result in inefficient computer codes [Czako et al., J. Chem. Phys. 122, 024101 (2005)]. It is shown here how the optimal-generalized FBR method can be simplified in the case of nondirect product bases of structures phi(n) (l)(s)f(l)(u), chi(m) (l)(t)phi(n) (l)(s)f(l)(u), etc. As a result the commonly used symmetric FBR is recovered and simplified nonsymmetric FBRs utilizing very special l-dependent grids are obtained. The nonsymmetric FBRs are more general than the symmetric FBR in that they can be employed efficiently even when the functions phi(n) (l)(s) [and/or chi(m) (l)(t)] are DVR functions of the infinite type. Arithmetic operation counts and a simple numerical example presented show unambiguously that setting up the Hamiltonian matrix requires significantly less computer time when using one of the proposed nonsymmetric FBRs than that in the symmetric FBR. Therefore, application of this nonsymmetric FBR is more efficient than that of the symmetric FBR when one wants to diagonalize the Hamiltonian matrix either by a direct or via a basis-set contraction method. Enormous decrease of computer time can be achieved, with respect to a direct application of the optimal-generalized FBR, by employing one of the simplified nonsymmetric FBRs as is demonstrated in noniterative calculations of the low-lying vibrational energy levels of the H3+ molecular ion. The arithmetic operation counts of the Hamiltonian matrix vector products and the properties of a recently developed diagonalization method [Andreozzi et al., J. Phys. A Math. Gen. 35, L61 (2002)] suggest that the nonsymmetric FBR applied along with this particular diagonalization method is suitable to large scale iterative calculations. Whether or not the nonsymmetric FBR is competitive with the symmetric FBR in large-scale iterative calculations still has to be investigated numerically.     

A discrete look at localization

This article introduces a set of localized orthonormal functions to serve as basis functions for quantum calculations. They are defined to be eigenfunctions of the position operator in a function space. Their properties, including their variances, for a one-dimensional system are developed. The application to simple harmonic motion is considered as an example and, in particular, the time evolution of an initially localized function is calculated and shown to be periodic. The theory can be interpreted as producing a discrete quantization of space with Hamiltonian interactions that are predominantly between nearest neighbors. These functions can also be used in approximate calculations. To illustrate their accuracy, the example of a Morse oscillator treated as a perturbation of a harmonic oscillator is reconsidered. It is shown that the localized functions in a variational calculation lead to a result that is a good approximation for the lowest states. Furthermore, the use of a wave function that is defined only at discrete points can be justified as the first approximation to this, so that its accuracy can also be discussed. © 1995 John Wiley & Sons, Inc.     

Propagation with distributed Gaussians as a sparse, adaptive basis for higher-dimensional quantum dynamics

A simple quantum wavepacket propagation algorithm is presented, designed to produce a very compact, non-direct product representation in higher-dimensional cases. Instead of moving basis functions around, localized basis functions at pre-defined centers are added to and deleted from the representation, generating an active basis function set strictly localized to the region where the moving wavepacket has significantly non-zero values. Simple one-dimensional examples prove this property, as well as the ability of the algorithm to accommodate splitting and rejoining of an arbitrary number of wavefunction pieces, and tunnelling through potential energy barriers. It is argued that future applications to higher-dimensional examples will be less expensive than with traditional direct-product bases, since making the basis adaptive has a lower scaling than the elementary steps necessary for any propagation algorithm itself.     

Superior performance of Mukherjee’s state-specific multi-reference coupled-cluster theory at the singles and doubles truncation scheme with localized active orbitals

The state-specific multi-reference coupled-cluster (SS-MRCC) theory of Mukherjee et al., in its singles and doubles truncation scheme (SS-MRCCSD), misses important couplings between the virtual functions reached by single and double excitations from different model functions. Since the SS-MRCC theory is not invariant with respect to the transformations among the active orbitals, the results are dependent on the active orbitals chosen. We demonstrate in this paper with results for potential energy curves for several example molecules involving single and multiple bond dissociation that the performance of SS-MRCCSD is significantly improved if localized active orbitals are used. The improvement is remarkable both in terms of the non-parallelity error and the magnitude of correlation energy recovered vis-a-vis the full configuration interaction results with the same basis set. The results bolster our claim that SS-MRCCSD with localized orbitals is an accurate general theory for potential energy surfaces.     

Numerically Integrable Confined Basis Functions for Band Structure Calculations

Numerical atomic and Slater-type radial functions are matched with polynomials which vanish at a predetermined radius. Outside this radius the functions are identically equal to zero. Such functions with a different number of derivatives continuous at the cutoff and matching radii were used for band structure calculations of copper. The utilization of Gauss quadratures, constructed without accounting for discontinuity of the basis functions, was shown to ensure the required accuracy of numerical integration for matrix element calculations with the basis functions that have three lowest derivatives discontinuous at the cutoff and matching radii. In the case of matching radius variations over a wide range, the quality of the basis set is virtually independent of this value. Increasing the cutoff radius improves the basis quality. Moreover, a basis limit is reached even at the cutoff radius smaller than a two-fold distance between the nearest neighbors. The resulting basis set is not inferior in quality to conventional bases of atomic or Slater-type functions and, in contrast to conventional functions, requires no special techniques of numerical integration. The proposed basis set enables a reduction of computational time by a factor of 2 approximately.     

Quantum theory of chemical reactions in the presence of electromagnetic fields

We present a theory for rigorous quantum scattering calculations of probabilities for chemical reactions of atoms with diatomic molecules in the presence of an external electric field. The approach is based on the fully uncoupled basis set representation of the total wave function in the space-fixed coordinate frame, the Fock-Delves hyperspherical coordinates, and the adiabatic partitioning of the total Hamiltonian of the reactive system. The adiabatic channel wave functions are expanded in basis sets of hyperangular functions corresponding to different reaction arrangements, and the interactions with external fields are included in each chemical arrangement separately. We apply the theory to examine the effects of electric fields on the chemical reactions of LiF molecules with H atoms and HF molecules with Li atoms at low temperatures and show that electric fields may enhance the probability of chemical reactions and modify reactive scattering resonances by coupling the rotational states of the reactants. Our preliminary results suggest that chemical reactions of polar molecules at temperatures below 1 K can be selectively manipulated with dc electric fields and microwave laser radiation.     

An ab initio method for approximation of the frozen molecular fragment

An ab initio method for calculation on many-electron molecular systems with the approximation of the inactive part of a molecule by frozen molecular fragment is presented. In the following method the SCF calculations are performed in two series. First the molecular orbitals resulting from the first SCF calculation (modest basis set) are localized. In the second SCF run, the basis set is extended for the active part of the molecule, while molecular orbitals of the inactive part, selected from the localized set, are kept frozen. The results are in good agreement with the extended basis set calculation.