Ab initio molecular dynamics using hybrid density functionals

Ab initio molecular dynamics simulations with hybrid density functionals have so far found little application due to their computational cost. In this work, an implementation of the Hartree-Fock exchange is presented that is specifically targeted at ab initio molecular dynamics simulations of medium sized systems. We demonstrate that our implementation, which is available as part of the CP2K/Quickstep program, is robust and efficient. Several prescreening techniques lead to a linear scaling cost for integral evaluation and storage. Integral compression techniques allow for in-core calculations on systems containing several thousand basis functions. The massively parallel implementation respects integral symmetry and scales up to hundreds of CPUs using a dynamic load balancing scheme. A time-reversible multiple time step scheme, exploiting the difference in computational efficiency between hybrid and local functionals, brings further time savings. With extensive simulations of liquid water, we demonstrate the ability to perform, for several tens of picoseconds, ab initio molecular dynamics based on hybrid functionals of systems in the condensed phase containing a few thousand Gaussian basis functions.     

Density functional investigations of the properties and thermochemistry of UF6 and UF5 using valence-electron and all-electron approaches

The structural properties and thermochemistry of UF6 and UF5 have been investigated using both Hartree-Fock and density functional theory (DFT) approximations. Within the latter approach, the local spin-density approximation, the generalized gradient approximation, and hybrid density functionals were considered. To describe the uranium atom we employed small-core (60 electrons) and large-core (78 electrons) relativistic effective core potentials (RECPs), as well as the all-electron approximation based on the two-component third-order Douglas-Kroll-Hess Hamiltonian. For structural properties, we obtained very good agreement with experiment with DFT and both large and small-core RECPs. The best match with experiment is given by the hybrid functionals with the small-core RECP. The bond dissociation energy (BDE) was obtained from the relative energies of the fragments [UF6 --> UF5 + F], corrected for zero-point energy and spin-orbit interaction. Very good agreement was found between the BDE obtained from all-electron calculations and those calculated with the small-core RECP, while those from the large-core RECP are off by more than 50%. In order to obtain good agreement with experiment in the BDE it is imperative to work with hybrid density functionals and a small-core RECP.     

Assessment of new meta and hybrid meta density functionals for predicting the geometry and binding energy of a challenging system: the dimer of H2S and benzene

Noncovalent interactions of a hydrogen bond donor with an aromatic pi system present a challenge for density functional theory, and most density functionals do not perform well for this kind of interaction. Here we test seven recent density functionals from our research group, along with the popular B3LYP functional, for the dimer of H 2S with benzene. The functionals considered include the four new meta and hybrid meta density functionals of the M06 suite, three slightly older hybrid meta functionals, and the B3LYP hybrid functional, and they were tested for their abilities to predict the dissociation energies of three conformations of the H 2S-benzene dimer and to reproduce the key geometric parameters of the equilibrium conformation of this dimer. All of the functionals tested except B3LYP correctly predict which of the three conformations of the dimer is the most stable. The functionals that are best able to reproduce the geometry of the equilibrium conformation of the dimer with a polarized triple-zeta basis set are M06-L, PWB6K, and MPWB1K, each having a mean unsigned relative error across the two experimentally verifiable geometric parameters of only 8%. The success of M06-L is very encouraging because it is a local functional, which reduces the cost for large simulations. The M05-2X functional yields the most accurate binding energy of a conformation of the dimer for which a binding energy calculated at the CCSD(T) level of theory is available; M05-2X gives a binding energy for the system with a difference of merely 0.02 kcal/mol from that obtained by the CCSD(T) calculation. The M06 functional performs well in both categories by yielding a good representation of the geometry of the equilibrium structure and by giving a binding energy that is only 0.19 kcal/mol different from that calculated by CCSD(T). We conclude that the new generation of density functionals should be useful for a variety of problems in biochemistry and materials where aromatic functional groups can serve as hydrogen bond acceptors.     

Simulations of liquid ammonia based on the combined quantum mechanical/molecular mechanical (QM/MM) approach

Two combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, namely, HF/MM and B3LYP/MM, have been performed to investigate the local structure and dynamics of liquid ammonia. The most interesting region, a sphere containing a central reference molecule and all its nearest surrounding molecules (first coordination shell), was treated by the Hartree-Fock (HF) and hybrid density functional B3LYP methods, whereas the rest of the system was described by the classical pair potentials. On the basis of both HF and B3LYP methods, it is observed that the hydrogen bonding in this peculiar liquid is weak. The structure and dynamics of this liquid are suggested to be determined by the steric packing effects, rather than by the directional hydrogen bonding interactions. Compared to previous empirical as well as Car-Parrinello (CP) molecular dynamics studies, our QM/MM simulations provide detailed information that is in better agreement with experimental data.     

Liquid structures of water, methanol, and hydrogen fluoride at ambient conditions from first principles molecular dynamics simulations with a dispersion corrected density functional

Using first principles molecular dynamics simulations in the isobaric-isothermal ensemble (T = 300 K, p = 1 atm) with the Becke-Lee-Yang-Parr exchange/correlation functional and a dispersion correction due to Grimme, the hydrogen bonding networks of pure liquid water, methanol, and hydrogen fluoride are probed. Although an accurate density is found for water with this level of electronic structure theory, the average liquid densities for both hydrogen fluoride and methanol are overpredicted by 50 and 25%, respectively. The radial distribution functions indicate somewhat overstructured liquid phases for all three compounds. The number of hydrogen bonds per molecule in water is about twice as high as for methanol and hydrogen fluoride, though the ratio of cohesive energy over number of hydrogen bonds is lower for water. An analysis of the hydrogen-bonded aggregates revealed the presence of mostly linear chains in both hydrogen fluoride and methanol, with a few stable rings and chains spanning the simulation box in the case of hydrogen fluoride. Only an extremely small fraction of smaller clusters was found for water, indicating that its hydrogen bond network is significantly more extensive. A special form of water with on average about two hydrogen bonds per molecule yields a hydrogen-bonding environment significantly different from the other two compounds.     

Excitonic effects in a time-dependent density functional theory

Excited state properties of one-dimensional molecular materials are dominated by many-body interactions resulting in strongly bound confined excitons. These effects cannot be neglected or treated as a small perturbation and should be appropriately accounted for by electronic structure methodologies. We use adiabatic time-dependent density functional theory to investigate the electronic structure of one-dimensional organic semiconductors, conjugated polymers. Various commonly used functionals are applied to calculate the lowest singlet and triplet state energies and oscillator strengths of the poly(phenylenevinylene) and ladder-type (poly)(para-phenylene) oligomers. Local density approximations and gradient-corrected functionals cannot describe bound excitonic states due to lack of an effective attractive Coulomb interaction between photoexcited electrons and holes. In contrast, hybrid density functionals, which include long-range nonlocal and nonadiabatic corrections in a form of a fraction of Hartree-Fock exchange, are able to reproduce the excitonic effects. The resulting finite exciton sizes are strongly dependent on the amount of the orbital exchange included in the functional.     

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.     

A first principles simulation study of fluctuations of hydrogen bonds and vibrational frequencies of water at liquid-vapor interface

We present a theoretical study of the structure and dynamics of water-vapor interface by means of ab initio molecular dynamics simulations. The inhomogeneous density, hydrogen bond and orientational profiles, voids and vibrational frequency distributions are investigated. We have also studied various dynamical properties of the interface such as diffusion, orientational relaxation, hydrogen bond dynamics and vibrational frequency fluctuations. The diffusion and orientational relaxation of water molecules are found to be faster at the interface which can be correlated with the voids present in the system. The hydrogen bond dynamics, however, is found to be slightly slower at the interface than that in bulk water. The correlations of hydrogen bond relaxation with the dynamics of vibrational frequency fluctuations are also discussed.     

Efficient hybrid density functional calculations in solids: assessment of the Heyd-Scuseria-Ernzerhof screened Coulomb hybrid functional

The present work introduces an efficient screening technique to take advantage of the fast spatial decay of the short range Hartree-Fock (HF) exchange used in the Heyd-Scuseria-Ernzerhof (HSE) screened Coulomb hybrid density functional. The screened HF exchange decay properties and screening efficiency are compared with traditional hybrid functional calculations on solids. The HSE functional is then assessed using 21 metallic, semiconducting, and insulating solids. The examined properties include lattice constants, bulk moduli, and band gaps. The results obtained with HSE exhibit significantly smaller errors than pure density functional theory (DFT) calculations. For structural properties, the errors produced by HSE are up to 50% smaller than the errors of the local density approximation, PBE, and TPSS functionals used for comparison. When predicting band gaps of semiconductors, we found smaller errors with HSE, resulting in a mean absolute error of 0.2 eV (1.3 eV error for all pure DFT functionals). In addition, we present timing results which show the computational time requirements of HSE to be only a factor of 2-4 higher than pure DFT functionals. These results make HSE an attractive choice for calculations of all types of solids.     

Solvation structure of hydroxyl radical by Car-Parrinello molecular dynamics

Car-Parrinello molecular dynamics simulations of a hydroxyl radical in liquid water have been performed. Structural and dynamical properties of the solvated structure have been studied in details. The partial atom-atom radial distribution functions for the hydrated hydroxyl do not show drastic differences with the radial distribution functions for liquid water. The OH is found to be a more active hydrogen bond donor and acceptor than the water molecule, but the accepted hydrogen bonds are much weaker than for the hydroxide OH- ion. The first solvation shell of the OH is less structured than the water's one and contains a considerable fraction of water molecules that are not hydrogen bonded to the hydroxyl. Part of them are found to come closer to the solvated radical than the hydrogen bonded molecules do. The lifetime of the hydrogen bonds accepted by the hydroxyl is found to be shorter than the hydrogen bond lifetime in water. A hydrogen transfer between a water molecule and the OH radical has been observed, though it is a much rarer event than a proton transfer between water and an OH- ion. The velocity autocorrelation power spectrum of the hydroxyl hydrogen shows the properties both of the OH radical in clusters and of the OH- ion in liquid.     

Computational study of structural and dynamical properties of formamide-water mixtures

A molecular dynamics simulation study of structural and dynamical properties in liquid mixtures of formamide and water is presented. Site-site radial pair distribution functions, local mole fractions, pair energy distributions, and tetrahedral orientational order are the quantities analyzed to investigate the local structure in the simulated mixtures, along with a review of the intermolecular structure in terms of the distribution of hydrogen bonds. Our results indicate that there is a substitution of formamide molecules by water in the hydrogen bonds and a formation of a common hydrogen bond network. By analyzing the extent of tetrahedral order in the liquid as a function of composition, it is observed that whereas the tetrahedral network of liquid water is progressively lost by increasing the formamide concentration, the water structure within the first coordination shell is preserved and somewhat enhanced. The hydrogen-bond mean lifetimes were estimated by performing a time integration of the autocorrelation functions of bond occupation numbers. The lifetimes associated with hydrogen bonds between water, formamide, and interspecies pairs are found to increase with increasing formamide concentration. The lifetimes of the water hydrogen bonds show the largest variations, supporting the picture of an enhancement of the water structure among the nearest neighbors within the first coordination shell. We have used two different force field models for water, SPC/E [J. C. Berendsen et al., J. Phys. Chem. 91, 6269 (1987)] and TIP4P/2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005)]. Our results for structural and dynamical properties yield very small differences between those models, the TIP4P/2005 predicting a slightly more structured liquid and, consequently, exhibiting a slightly slower translational and librational dynamics.     

Hydrogen bonded structure and dynamics of liquid-vapor interface of water-ammonia mixture: an ab initio molecular dynamics study

We have carried out ab initio molecular dynamics simulations of a liquid-vapor interfacial system consisting of a mixture of water and ammonia molecules. We have made a detailed analysis of the structural and dynamical properties of the bulk and interfacial regions of the mixture. Among structural properties, we have looked at the inhomogeneous density profiles of water and ammonia molecules, hydrogen bond distributions, orientational profiles, and also vibrational frequency distributions of bulk and interfacial molecules. It is found that the interfacial molecules show preference for specific orientations so as to form water-ammonia hydrogen bonds at the interface with ammonia as the acceptor. The structure of the system is also investigated in terms of inter-atomic voids present in the system. Among the dynamical properties, we have calculated the diffusion, orientational relaxation, hydrogen bond dynamics, and vibrational spectral diffusion in bulk and interfacial regions. It is found that the diffusion and orientation relaxation of the interfacial molecules are faster than those of the bulk. However, the hydrogen bond lifetimes are longer at the interface which can be correlated with the time scales found from the decay of frequency time correlations.     

Hydrogen bond networks in water and methanol with varying interaction strengths

Metropolis Monte Carlo simulations of hydrogen-bonded liquids (water and methanol) were performed with the well tested effective pair potentials TIP5P and OPLS. The Coulomb contribution for the interaction potential was damped by a factor η varied from 1 to 0.49 for water and 1 to 0.15 for methanol. As a result, the networks formed by the hydrogen-bonded molecules presented interesting properties as a function of η, including small-world patterns and percolation transitions. These complex networks were analyzed by local (clustering coefficients, average degrees), semi-global (path lengths) and global (spectral densities) properties, and islands statistics. From these properties, small-world behavior was found for η in the range 0.60-0.75 for both liquids, interestingly independent of the molecular structure of the liquid. Phase transition behavior was observed for the average degrees and the clustering coefficient curves with critical values at 0.55 for water and 0.34 for methanol. Macroscopic properties such as mass density and vaporization enthalpy were also parametrically dependent on η and they presented phase transition behavior that coincides with the critical values obtained from the topological analysis. This is probably the first time that such phase transitions are observed for these quantities and shows a direct relation between macroscopic properties and topological features of hydrogen bond networks.     

First principles calculations on structure, bonding, and vibrational frequencies of SiP2

Pyrite type SiP(2) is reinvestigated by first principles calculations on various levels of functionals including local density approximation, generalized gradient approximation, Becke-Lee-Yang-Parr hybrid functional, and the Hartree-Fock method. SiP(2) is seen as a model compound with molecular [P-P] entities and [SiP(6)] octahedra. Structure and bonding are addressed by electronic structure calculations. Special attention is spent on P-P and Si-P bonds in terms of bond lengths and respective stretching modes from simulated Raman spectra. The electronic structure is analyzed in both direct and momentum space by the electron localization function and site projected density of states. The main goals of this work are to understand the nature of chemical bonding in SiP(2) and to compare and contrast the different methods of calculation.     

The AM05 density functional applied to solids

We show that the AM05 functional [Armiento and Mattsson, Phys. Rev. B 72, 085108 (2005)] has the same excellent performance for solids as the hybrid density functionals tested in Paier et al. [J. Chem. Phys. 124, 154709 (2006); 125, 249901 (2006)]. This confirms the original finding that AM05 performs exceptionally well for solids and surfaces. Hartree-Fock hybrid calculations are typically an order of magnitude slower than local or semilocal density functionals such as AM05, which is of a regular semilocal generalized gradient approximation form. The performance of AM05 is on average found to be superior to selecting the best of local density approximation and PBE for each solid. By comparing data from several different electronic-structure codes, we have determined that the numerical errors in this study are equal to or smaller than the corresponding experimental uncertainties.     

Databases for transition element bonding: metal-metal bond energies and bond lengths and their use to test hybrid, hybrid meta, and meta density functionals and generalized gradient approximations

We propose a data set of bond lengths for 8 selected transition metal dimers (Ag(2), Cr(2), Cu(2), CuAg, Mo(2), Ni(2), V(2), and Zr(2)) and another data set containing their atomization energies and the atomization energy of ZrV, and we use these for testing density functional theory. The molecules chosen for the test sets were selected on the basis of the expected reliability of the data and their ability to constitute a diverse and representative set of transition element bond types while the data sets are kept small enough to allow for efficient testing of a large number of computational methods against a very reliable subset of experimental data. In this paper we test 42 different functionals: 2 local spin density approximation (LSDA) functionals, 12 generalized gradient approximation (GGA) methods, 13 hybrid GGAs, 7 meta GGA methods, and 8 hybrid meta GGAs. We find that GGA density functionals are more accurate for the atomization energies of pure transition metal systems than are their meta, hybrid, or hybrid meta analogues. We find that the errors for atomization energies and bond lengths are not as large if we limit ourselves to dimers with small amounts of multireference character. We also demonstrate the effects of increasing the fraction of Hartree-Fock exchange in multireference systems by computing the potential energy curve for Cr(2) and Mo(2) with several functionals. We also find that BLYP is the most accurate functional for bond energies and is reasonably accurate for bond lengths. The methods that work well for transition metal bonds are found to be quite different from those that work well for organic and other main group chemistry.     

Embedded,graph-theoretically defined many-body approximations for wavefunction-in-DFT and DFT-in-DFT: Applications to gas- and condensed-phase ab initio molecular dynamics,and potential surfaces for quantum nuclear effects

We present a graph-theoretic approach to adaptively compute many-body approximations in an efficient manner to perform (a) accurate post-Hartree–Fock (HF) ab initio molecular dynamics (AIMD) at density functional theory (DFT) cost for medium- to large-sized molecular clusters, (b) hybrid DFT electronic structure calculations for condensed-phase simulations at the cost of pure density functionals, (c) reduced-cost on-the-fly basis extrapolation for gas-phase AIMD and condensed phase studies, and (d) accurate post-HF-level potential energy surfaces at DFT cost for quantum nuclear effects. The salient features of our approach are ONIOM-like in that (a) the full system (cluster or condensed phase) calculation is performed at a lower level of theory (pure DFT for condensed phase or hybrid DFT for molecular systems), and (b) this approximation is improved through a correction term that captures all many-body interactions up to any given order within a higher level of theory (hybrid DFT for condensed phase; CCSD or MP2 for cluster), combined through graph-theoretic methods. Specifically, a region of chemical interest is coarse-grained into a set of nodes and these nodes are then connected to form edges based on a given definition of local envelope (or threshold) of interactions. The nodes and edges together define a graph, which forms the basis for developing the many-body expansion. The methods are demonstrated through (a) ab initio dynamics studies on protonated water clusters and polypeptide fragments, (b) potential energy surface calculations on one-dimensional water chains such as those found in ion channels, and (c) conformational stabilization and lattice energy studies on homogeneous and heterogeneous surfaces of water with organic adsorbates using two-dimensional periodic boundary conditions.     

Investigating pressure effects on structural and dynamical properties of liquid methanol with many-body interactions

Molecular-dynamics simulations utilizing a many-body potential was used to study the pressure dependence of structural and dynamical properties for liquid methanol. The liquid density as a function of pressure agreed well with experiment, and a combination of radial and angular distribution functions were used to analyze molecular structure. From these distribution functions, it was observed that hydrogen bond strength increased with increasing pressure. This observation coincided with an increase in the molecular dipole as a function of pressure, having a significant effect on the observed increased hydrogen bond strength. Also, methanols were found to more strongly favor exactly two hydrogen bonds, with fewer methanols of zero, one, or three hydrogen bonds present at higher pressures. Furthermore, a majority of the compression with increased pressure was found to occur in regions perpendicular to the methanol hydrogen-oxygen bond vector. This was the case despite hydrogen-oxygen nonbonded distances between hydrogen bonding species being shorter, but their stiffer oxygen-hydrogen-(nonbonded) oxygen angle offsets this, resulting in their oxygen-oxygen distances being relatively unaffected. The methanol translational diffusion decreased significantly with increased pressure, while the rotational diffusion decreased at a similar magnitude around the oxygen-hydrogen and oxygen-carbon bond vectors, despite having very different overall diffusion. Finally, the hydrogen bond lifetime increased significantly with pressure, owing to the increased hydrogen bond strength, and the slower translational and rotational dynamics.