Acceleration schemes for ab initio molecular-dynamics simulations and electronic-structure calculations

First principles molecular dynamics calculations of the mechanical properties of endofullerenes containing noble gas atoms or small molecules

The mechanical properties of endofullerenes have been investigated by performing compression tests using finite temperature first principles molecular dynamics calculations. We considered various X@C₆₀ systems, with X a single noble gas atom (He, Ne, Ar, Kr, or Xe), small molecules (H₂O, CH₄), or small helium clusters. In the absence of compression, it is observed that there is no or at best a negligible effect of X on the properties of C₆₀. The compression simulations revealed several original findings. First, the influence of X on the stiffness of X@C₆₀ can be quantified, although it is at most 12% for the studied cases. Next, both energy and contact force variations as a function of strain are demonstrated to depend on X. However, this is not the case for the yield strain and for the failure mechanism of the C₆₀ shell. Finally, it is shown that the X@C₆₀ compression could bring X to be in a high stress state. In the specific cases of H₂O and CH₄ molecules, a mechanism of stress assisted dissociation is observed.

Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory-Based Ab Initio Molecular Dynamics II: Extensions to the Isobaric-Isoenthalpic and Isobaric-Isothermal Ensembles

In the previous paper of this series [Ko, H.-Y. et al. J. Chem. Theory Comput. 2020, 16, 3757-3785], we presented a theoretical and algorithmic framework based on a localized representation of the occupied space that exploits the inherent sparsity in the real-space evaluation of the exact exchange (EXX) interaction in finite-gap systems. This was accompanied by a detailed description of exx, a massively parallel hybrid message-passing interface MPI/OpenMP implementation of this approach in Quantum ESPRESSO (QE) that enables linear scaling hybrid density functional theory (DFT)-based ab initio molecular dynamics (AIMD) in the microcanonical/canonical (NVE/NVT) ensembles of condensed-phase systems containing 500-1000 atoms (in fixed orthorhombic cells) with a wall time cost comparable to semi-local DFT. In this work, we extend the current capabilities of exx to enable hybrid DFT-based AIMD simulations of large-scale condensed-phase systems with general and fluctuating cells in the isobaric-isoenthalpic/isobaric-isothermal (NpH/NpT) ensembles. The theoretical extensions to this approach include an analytical derivation of the EXX contribution to the stress tensor for systems in general simulation cells with a computational complexity that scales linearly with system size. The corresponding algorithmic extensions to exx include optimized routines that (i) handle both static and fluctuating simulation cells with non-orthogonal lattice symmetries, (ii) solve Poisson's equation in general/non-orthogonal cells via an automated selection of the auxiliary grid directions in the Natan-Kronik representation of the discrete Laplacian operator, and (iii) evaluate the EXX contribution to the stress tensor. Using this approach, we perform a case study on a variety of condensed-phase systems (including liquid water, a benzene molecular crystal polymorph, and semi-conducting crystalline silicon) and demonstrate that the EXX contributions to the energy and stress tensor simultaneously converge with an appropriate choice of exx parameters. This is followed by a critical assessment of the computational performance of the extended exx module across several different high-performance computing architectures via case studies on (i) the computational complexity due to lattice symmetry during NpT simulations of three different ice polymorphs (i.e., ice Ih, II, and III) and (ii) the strong/weak parallel scaling during large-scale NpT simulations of liquid water. We demonstrate that the robust and highly scalable implementation of this approach in the extended exx module is capable of evaluating the EXX contribution to the stress tensor with negligible cost (<1%) as well as all other EXX-related quantities needed during NpT simulations of liquid water (with a very tight 150 Ry planewave cutoff) in ≈5.2 s ((H₂O)₁₂₈) and ≈6.8 s ((H₂O)₂₅₆) per AIMD step. As such, the extended exx module presented in this work brings us another step closer to routinely performing hybrid DFT-based AIMD simulations of sufficient duration for large-scale condensed-phase systems across a wide range of thermodynamic conditions.

Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory Based Ab Initio Molecular Dynamics. 1. Theory, Algorithm, and Performance

By including a fraction of exact exchange (EXX), hybrid functionals reduce the self-interaction error in semilocal density functional theory (DFT) and thereby furnish a more accurate and reliable description of the underlying electronic structure in systems throughout biology, chemistry, physics, and materials science. However, the high computational cost associated with the evaluation of all required EXX quantities has limited the applicability of hybrid DFT in the treatment of large molecules and complex condensed-phase materials. To overcome this limitation, we describe a linear-scaling approach that utilizes a local representation of the occupied orbitals (e.g., maximally localized Wannier functions (MLWFs)) to exploit the sparsity in the real-space evaluation of the quantum mechanical exchange interaction in finite-gap systems. In this work, we present a detailed description of the theoretical and algorithmic advances required to perform MLWF-based ab initio molecular dynamics (AIMD) simulations of large-scale condensed-phase systems of interest at the hybrid DFT level. We focus our theoretical discussion on the integration of this approach into the framework of Car-Parrinello AIMD, and highlight the central role played by the MLWF-product potential (i.e., the solution of Poisson's equation for each corresponding MLWF-product density) in the evaluation of the EXX energy and wave function forces. We then provide a comprehensive description of the exx algorithm implemented in the open-source Quantum ESPRESSO program, which employs a hybrid MPI/OpenMP parallelization scheme to efficiently utilize the high-performance computing (HPC) resources available on current- and next-generation supercomputer architectures. This is followed by a critical assessment of the accuracy and parallel performance (e.g., strong and weak scaling) of this approach when AIMD simulations of liquid water are performed in the canonical (NVT) ensemble. With access to HPC resources, we demonstrate that exx enables hybrid DFT-based AIMD simulations of condensed-phase systems containing 500-1000 atoms (e.g., (H₂O)₂₅₆) with a wall time cost that is comparable to that of semilocal DFT. In doing so, exx takes us one step closer to routinely performing AIMD simulations of complex and large-scale condensed-phase systems for sufficiently long time scales at the hybrid DFT level of theory.

Electronic structure, optoelectronic properties and enhanced photocatalytic response of GaN–GeC van der Waals heterostructures: a first principles study

In this work, we systematically studied the electronic structure and optical characteristics of van der Waals (vdW) heterostructure composed of a single layer of GaN and GeC using first principles calculations. The GaN–GeC vdW heterostructure exhibits indirect band gap semiconductor properties and possesses type-II energy band arrangement, which will help the separation of photogenerated carriers and extend their lifetime. In addition, the band edge positions of the GaN–GeC heterostructure meet both the requirements of water oxidation and reduction energy, indicating that the photocatalysts have the potential for water decomposition. The GaN–GeC heterostructure shows obvious absorption peaks in the visible region, leading to the efficient use of solar energy. Tensile and compressive strains of up to 10% are also proposed. Tensile strain leads to an increase in the blue shift of optical absorption, whereas a red shift is observed in the case of the compressive strain. These fascinating characteristics make the GaN–GeC vdW heterostructure a highly effective photocatalyst for water splitting.

In this work, we systematically studied the electronic structure and optical characteristics of van der Waals (vdW) heterostructure composed of a single layer of GaN and GeC using first principles calculations.

Importance of van der Waals effects on the hydration of metal ions from the Hofmeister series

The van der Waals (vdW) interaction plays a crucial role in the description of liquid water. Based on ab initio molecular dynamics simulations, including the non-local and fully self-consistent density-dependent implementation of the Tkatchenko-Scheffler dispersion correction, we systematically studied the aqueous solutions of metal ions (K⁺, Na⁺, and Ca²⁺) from the Hofmeister series. Similar to liquid water, the vdW interactions strengthen the attractions among water molecules in the long-range, leading to the hydrogen bond networks softened in all the ion solutions. However, the degree that the hydration structure is revised by the vdW interactions is distinct for different ions, depending on the strength of short-range interactions between the hydrated ion and surrounding water molecules. Such revisions by the vdW interactions are important for the understanding of biological functionalities of ion channels.

Direct inversion of the iterative subspace with contracted planewave basis functions

Ways to reduce the computational cost of periodic electronic structure calculations by using basis functions corresponding to linear combinations of planewaves have been examined recently. These contracted planewave (CPW) basis functions correspond to Fourier series representations of atom-centered basis functions, and thus provide access to some beneficial properties of planewave (PW) and localized basis functions. This study reports the development and assessment of a direct inversion of the iterative subspace (DIIS) method that employs unique properties of CPW basis functions to efficiently converge electronic wavefunctions. This method relies on access to a PW-based representation of the electronic structure to provide a means of efficiently evaluating matrix-vector products involving the application of the Fock matrix to the occupied molecular orbitals. These matrix-vector products are transformed into a form permitting the use of direct diagonalization techniques and DIIS methods typically employed with atom-centered basis sets. The abilities of this method are assessed through periodic Hartree-Fock calculations of a range of molecules and solid-state systems. The results show that the method reported in this study is approximately five times faster than CPW-based calculations in which the entire Fock matrix is calculated. This method is also found to be weakly dependent upon the size of the basis set, thus permitting the use of larger CPW basis sets to increase variational flexibility with a minor impact on computational performance. © 2018 Wiley Periodicals, Inc.

Ab initio theory and modeling of water

Water is of the utmost importance for life and technology. However, a genuinely predictive ab initio model of water has eluded scientists. We demonstrate that a fully ab initio approach, relying on the strongly constrained and appropriately normed (SCAN) density functional, provides such a description of water. SCAN accurately describes the balance among covalent bonds, hydrogen bonds, and van der Waals interactions that dictates the structure and dynamics of liquid water. Notably, SCAN captures the density difference between water and ice Ih at ambient conditions, as well as many important structural, electronic, and dynamic properties of liquid water. These successful predictions of the versatile SCAN functional open the gates to study complex processes in aqueous phase chemistry and the interactions of water with other materials in an efficient, accurate, and predictive, ab initio manner.

Accelerating ab initio Molecular Dynamics and Probing the Weak Dispersive Forces in Dense Liquid Hydrogen

We propose an ab initio molecular dynamics method, capable of dramatically reducing the autocorrelation time required for the simulation of classical and quantum particles at finite temperatures. The method is based on an efficient implementation of a first order Langevin dynamics modified by means of a suitable, position dependent acceleration matrix S. Here, we apply this technique to both Lennard-Jones models, to demonstrate the accuracy and speeding-up of the sampling, and within a quantum Monte Carlo based wave function approach, for determining the phase diagram of high-pressure hydrogen with simulations much longer than the autocorrelation time. With the proposed method, we are able to equilibrate in a few hundred steps even close to the liquid-liquid phase transition (LLT). Within our approach, we find that the LLT transition is consistent with recent density functionals predicting a much larger transition pressure when the long range dispersive forces are taken into account.

The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water

In this work, we report the results of a series of density functional theory (DFT) based ab initio molecular dynamics (AIMD) simulations of ambient liquid water using a hierarchy of exchange-correlation (XC) functionals to investigate the individual and collective effects of exact exchange (Exx), via the PBE0 hybrid functional, non-local van der Waals/dispersion (vdW) interactions, via a fully self-consistent density-dependent dispersion correction, and an approximate treatment of nuclear quantum effects, via a 30 K increase in the simulation temperature, on the microscopic structure of liquid water. Based on these AIMD simulations, we found that the collective inclusion of Exx and vdW as resulting from a large-scale AIMD simulation of (H2O)128 significantly softens the structure of ambient liquid water and yields an oxygen-oxygen structure factor, SOO(Q), and corresponding oxygen-oxygen radial distribution function, gOO(r), that are now in quantitative agreement with the best available experimental data. This level of agreement between simulation and experiment demonstrated herein originates from an increase in the relative population of water molecules in the interstitial region between the first and second coordination shells, a collective reorganization in the liquid phase which is facilitated by a weakening of the hydrogen bond strength by the use of a hybrid XC functional, coupled with a relative stabilization of the resultant disordered liquid water configurations by the inclusion of non-local vdW/dispersion interactions. This increasingly more accurate description of the underlying hydrogen bond network in liquid water also yields higher-order correlation functions, such as the oxygen-oxygen-oxygen triplet angular distribution, POOO(θ), and therefore the degree of local tetrahedrality, as well as electrostatic properties, such as the effective molecular dipole moment, that are in much better agreement with experiment.

A combined experimental and theoretical study of ion solvation in liquid N-methylacetamide

Most current biomolecular simulations are based on potential energy functions that treat the electrostatic energy as a sum of pairwise Coulombic interactions between effective fixed atomic charges. This approximation, in which many-body induced polarization effects are included in an average way, is expected to be satisfactory for a wide range of systems, but less accurate for processes involving the transfer and partition of ions among heterogeneous environments. The limitations of these potential energy functions are perhaps most obvious in studies of ion permeation through membrane channels. In many cases, the pore is so narrow that the permeating ion must shed most of its surrounding water molecules and the large energetic loss due to dehydration must be compensated by coordination with protein atoms. Interactions of cations with protein backbone carbonyl oxygens, in particular, play a critical role in several important biological channels. As a first step toward meeting the challenge of developing an accurate explicit accounting for induced polarization effects, the present work combines experiments and computation to characterize the interactions of alkali and halide ions with N-methylacetamide chosen to represent the peptide bond. From solubility measurements, we extract the solvation free energies of KCl and NaCl in liquid N-methylacetamide. Polarizable models based on the Drude oscillator are then developed and compared with available experimental and ab initio data. The good agreement for a range of structural and thermodynamic properties in the gas and condensed phases suggests that the polarizable models provide an accurate representation of ion-amide interactions in biological systems.

Vibrational frequencies in Car-Parrinello molecular dynamics

Car-Parrinello molecular dynamics (CPMD) are widely used to investigate the dynamical properties of molecular systems. An important issue in such applications is the dependence of dynamical quantities such as molecular vibrational frequencies upon the fictitious orbital mass μ. Although it is known that the correct Born-Oppenheimer dynamics are recovered at zero μ, it is not clear how these dynamical quantities are to be rigorously extracted from CPMD calculations. Our work addresses this issue for vibrational frequencies. We show that when the system is sufficiently close to the ground state the calculated ionic vibrational frequencies are ω(M) = ω(0M)[1 -C(μ/M)] for small μ/M, where ω(0M) is the Born-Oppenheimer ionic frequency, M the ionic mass, and C a constant that depends upon the ion-orbital coupling force constants. Our analysis also provides a quantitative understanding of the orbital oscillation amplitudes, leading to a relationship between the adiabaticity of a system and the ion-orbital coupling constants. In particular, we show that there is a significant systematic dependence of calculated vibrational frequencies upon how close the CPMD trajectory is to the Born-Oppenheimer surface. We verify our analytical results with numerical simulations for N(2), Sn(2), and H/Si(100)-(2×1).

QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials

QUANTUM ESPRESSO is an integrated suite of computer codes for electronic-structure calculations and materials modeling, based on density-functional theory, plane waves, and pseudopotentials (norm-conserving, ultrasoft, and projector-augmented wave). The acronym ESPRESSO stands for opEn Source Package for Research in Electronic Structure, Simulation, and Optimization. It is freely available to researchers around the world under the terms of the GNU General Public License. QUANTUM ESPRESSO builds upon newly-restructured electronic-structure codes that have been developed and tested by some of the original authors of novel electronic-structure algorithms and applied in the last twenty years by some of the leading materials modeling groups worldwide. Innovation and efficiency are still its main focus, with special attention paid to massively parallel architectures, and a great effort being devoted to user friendliness. QUANTUM ESPRESSO is evolving towards a distribution of independent and interoperable codes in the spirit of an open-source project, where researchers active in the field of electronic-structure calculations are encouraged to participate in the project by contributing their own codes or by implementing their own ideas into existing codes.

Inflationary dynamics for matrix eigenvalue problems

Many fields of science and engineering require finding eigenvalues and eigenvectors of large matrices. The solutions can represent oscillatory modes of a bridge, a violin, the disposition of electrons around an atom or molecule, the acoustic modes of a concert hall, or hundreds of other physical quantities. Often only the few eigenpairs with the lowest or highest frequency (extremal solutions) are needed. Methods that have been developed over the past 60 years to solve such problems include the Lanczos algorithm, Jacobi-Davidson techniques, and the conjugate gradient method. Here, we present a way to solve the extremal eigenvalue/eigenvector problem, turning it into a nonlinear classical mechanical system with a modified Lagrangian constraint. The constraint induces exponential inflationary growth of the desired extremal solutions.

Understanding the dielectric properties of liquid amides from a polarizable force field

The role played by electronic polarization in the dielectric properties of liquid N-methyl acetamide (NMA) is examined using molecular dynamics simulations with a polarizable force field based on classical Drude oscillators. The model presented is the first force field shown to reproduce the anomalously large dielectric constant of liquid NMA. Details of the molecular polarizability are found to be important. For instance, all elements of the polarizability tensor, rather then just the trace, impact on the condensed phase properties. Two factors related to electronic polarizability are found to contribute to this large dielectric constant. First is the significant enhancement of the mean amide molecular dipole magnitude, which is 50% larger in the liquid than in the gas phase. Second is the consequent strong hydrogen bonding between molecular neighbors that enhances the orientational alignment of the molecular dipoles. Polarizable models of amide compounds that have two (acetamide) and zero (N,N-dimethyl acetamide) polar hydrogen-bond donor atoms are also investigated. Experimentally, the neat liquid dielectric constants at 373 K are 100 for NMA, 66 for acetamide and 26 for N,N-dimethyl acetamide. The polarizable models replicate this trend, predicting a dielectric constant of 92+/-5 for NMA, 66+/-3 for acetamide and 23+/-1 for N,N-dimethyl acetamide.

Theoretical Study of α-84 Phycocyanobilin Chromophore from the Thermophilic Cyanobacterium Synechococcus elongatus

Time-dependent density functional theory (TD-DFT) calculations were performed to obtain vertical excitation energies from the ground state to different low-lying singlet excited states of the protonated alpha-84 phycocyanobilin chromophore (alpha-84 PCBH(+)). It clearly emerges that three gradient-corrected approximation functionals (B3LYP, PBE0, and PBEPBE) show a similar description, confirming the proposed valence assignment of the strongest UV-vis absorption band at 618 nm. Moreover, our results show that there are not appreciable differences, in terms of excitation wavelength of the main peak, between the alpha-84 PCBH(+) chromophore and a model system in which the two propionic chains have not been taken into account. Finally, with the precise aim of investigating the effects of alpha-84 PCBH(+) conformational fluctuations on its electronic properties, vertical excitation energies obtained for the potential energy local minimum structure were also refined using a recently proposed TD-DFT/principal component analysis/Car-Parrinello molecular dynamics computational approach. Interestingly, and in line with previous results on another photosensitive complex, this study essentially suggests that interaction with the surrounding environment (protein matrix plus solvent molecules) coupled with the large amplitude fluctuation of the whole C-Phycocyanin (C-PC) pigment protein can affect the electronic properties of the alpha-84 PCB chromophore and therefore its biological activity.

Quasi-one-dimensional K-O chain in PTCDA thin films: evidence from first-principles calculations

Using density functional theory calculations we have found that K atoms in a PTCDA (3,4:9,10-perylenetetracarboxylic dianhydride) crystal form a quasi-one-dimensional (1D) K-O chain interacting with carboxylic oxygen of the terminal anhydride groups of PTCDA. The K-K distance in the chain (3.72 Angstrom) is commensurate to the periodicity of the organic semiconductor. We obtain that the K-O structure is stabilized by charge transfer from K to PTCDA molecules, forming prevalently ionic bonds: the electronic density of the chemistry induced gap states is essentially delocalized on the perylene core of PTCDA, while potassium appears spoiled of its charge. Band dispersion along the direction of molecular stack is evaluated to be 0.2 eV in pure PTCDA crystal and 0.5 eV in the K-doped system, confirming that the interaction occurs between different molecular planes.

A theoretical study of aqueous solvation of K comparing ab initio, polarizable, and fixed-charge models

The hydration of K(+) is studied using a hierarchy of theoretical approaches, including ab initio Born-Oppenheimer molecular dynamics and Car-Parrinello molecular dynamics, a polarizable force field model based on classical Drude oscillators, and a nonpolarizable fixed-charge potential based on the TIP3P water model. While models based more directly on quantum mechanics offer the possibility to account for complex electronic effects, polarizable and fixed-charges force fields allow for simulations of large systems and the calculation of thermodynamic observables with relatively modest computational costs. A particular emphasis is placed on investigating the sensitivity of the polarizable model to reproduce key aspects of aqueous K(+), such as the coordination structure, the bulk hydration free energy, and the self diffusion of K(+). It is generally found that, while the simple functional form of the polarizable Drude model imposes some restrictions on the range of properties that can simultaneously be fitted, the resulting hydration structure for aqueous K(+) agrees well with experiment and with more sophisticated computational models. A counterintuitive result, seen in Car-Parrinello molecular dynamics and in simulations with the Drude polarizable force field, is that the average induced molecular dipole of the water molecules within the first hydration shell around K(+) is slightly smaller than the corresponding value in the bulk. In final analysis, the perspective of K(+) hydration emerging from the various computational models is broadly consistent with experimental data, though at a finer level there remain a number of issues that should be resolved to further our ability in modeling ion hydration accurately.

First-principle molecular dynamics with ultrasoft pseudopotentials: parallel implementation and application to extended bioinorganic systems

We present a plane-wave ultrasoft pseudopotential implementation of first-principle molecular dynamics, which is well suited to model large molecular systems containing transition metal centers. We describe an efficient strategy for parallelization that includes special features to deal with the augmented charge in the contest of Vanderbilt's ultrasoft pseudopotentials. We also discuss a simple approach to model molecular systems with a net charge and/or large dipole/quadrupole moments. We present test applications to manganese and iron porphyrins representative of a large class of biologically relevant metalorganic systems. Our results show that accurate density-functional theory calculations on systems with several hundred atoms are feasible with access to moderate computational resources.

Ab initio molecular dynamics of liquid hydrogen chloride

We carried out an ab initio molecular dynamics simulation of liquid hydrogen chloride (l-HCl) at a temperature of 313 K. Comparison with inelastic neutron scattering data shows that the simulation achieves an overall good description of the structural correlations, improving significantly upon a description based on classical interaction potentials. Despite some minor differences between theory and experiment in the H-H partial structure factor, the simulation gives a description of the hydrogen bonding in impressive agreement with experiment, for both the amount and the bond-length distribution of the bonds. In the simulation, 40% of the molecules are nonbonded, while the hydrogen-bonded chains are short, principally consisting of dimers (25%) and trimers (15%). Neighboring molecules in the simulation are found to form L-shaped arrangements, like in the isolated (HCl)(2) dimer and in crystalline phases of HCl. The time correlation of the molecular-axis orientation is found to be characterized by a very short decay time (0.13 ps), consistent with the short length of the hydrogen-bonded chains. Other dynamical properties investigated in this work include the diffusion coefficient and the vibrational density of states. We evaluated the molecular dipole of the HCl molecule in the liquid using a definition based on the coupling of rotational modes to an external electric field. The average dipole moment (1.53 D) derived in this way is found to be considerably larger than for the isolated molecule (1.11 D). Our results show that the dipole moment in [script-l]-HCl undergoes large fluctuations, both in orientation and in modulus. Upon the onset of an external field, such dipole fluctuations concur to reduce the fluctuations of the dielectric response.

Ab Initio and Empirical Model MD Simulation Studies of Solvent Effects on the Properties of N-Methylacetamide along a cis−trans Isomerization Pathway

The properties of N-methylacetamide along a cis-trans isomerization pathway described by twisting about the C(O)-N bond are examined at finite temperature both in vacuo and in explicit water solvent. Two distinctly different theoretical descriptions, an ab initio (DFT-BLYP) and an empirical (CHARMM22) model, are studied in order to permit an assessment of the dominant forces active in the system. An analysis of the solvent structure at equilibrium and changes in solvation structure accompanying isomerization is, therefore, given for each model. Many-body polarization effects absent under CHARMM22 but present in the ab initio model are found to have a profound influence on the system. The electronic structure of the NMA molecule predicted by the ab initio method along the reaction coordinate is examined in order to shed further light on changes in peptide "partial-double" bond character [C(O)-N] as isomerization takes place. A new statistical-mechanical interpretation of the entropy change during a chemical reaction is presented to help interpret the thermochemistry of the simple reaction.

Structural Correlations and Motifs in Liquid Water at Selected Temperatures: Ab Initio and Empirical Model Predictions

To gain further insight into liquid water's structure and explore the role of different physical forces underlying the interaction between water molecules, the radial and angular structure of water is probed as a function of temperature for a carefully selected set of theoretical models. Simulations are performed with empirical rigid, empirical polarizable, empirical flexible with classical and quantum nuclei, and ab initio models with classical nuclei at 300 and 353 K and quantum nuclei at 300 K. The predicted radial distribution functions, spatial maps, and angular distributions of the neighboring water molecules are consistent with a model of liquid water in which water molecules are tetrahedrally coordinated. In addition, three-dimensional joint distribution functions are introduced and analyzed. By comparison of the functions obtained for hexagonal ice to those of liquid water, several thermally disordered, ice-like cluster structures are identified in the fluid.

Curvy-steps approach to constraint-free extended-Lagrangian ab initio molecular dynamics, using atom-centered basis functions: Convergence toward Born–Oppenheimer trajectories

A dynamical extension of the "curvy-steps" approach to linear-scaling self-consistent field calculations is presented, which yields an extended-Lagrangian formulation of ab initio molecular dynamics. An exponential parametrization of the one-electron density matrix, expressed in terms of atom-centered Gaussian basis functions, facilitates propagation along the manifold of density matrices in a geometrically correct fashion that automatically enforces idempotency constraints. The extended Lagrangian itself is constraint free, thus neither density matrix purification nor expensive, iterative solution for Lagrange multipliers is required. Propagation is highly efficient, and time steps compare favorably to those used in Car-Parrinello molecular dynamics simulations. The behavior of the method, especially with regard to the maintenance of adiabatic decoupling of nuclei and electrons, is examined for a sequence of diatomic molecules, and comparison is made to trajectories propagated on the converged Born-Oppenheimer surface. Certain claims to the contrary notwithstanding, our results demonstrate that vibrational frequencies may depend on the value of the fictitious mass parameter, even in an atom-centered basis. Light-atom stretching frequencies can be significantly redshifted, even when the nuclear and electronic energy scales are well separated. With a sufficiently small fictitious mass and a short time step, accurate frequencies can be obtained; we characterize appropriate values of these parameters for a wide range of vibrational frequencies.

Surface Reaction of Alkynes and Alkenes with H-Si(111): A Density Functional Theory Study

Recent experiments on the addition of alkene and alkyne molecules to H-terminated silicon surfaces have provided evidence for a surface chain reaction initiated at isolated Si dangling bonds and involving an intermediate carbon radical state, which, after abstraction of a hydrogen atom from a neighboring Si-H unit, transforms into a stable adsorbed species plus a new Si dangling bond. Using periodic density functional theory (DFT) calculations, together with an efficient method for determining reaction pathways, we have studied the initial steps of this chain reaction for a few different terminal alkynes and alkenes interacting with an isolated Si dangling bond on an otherwise H-saturated Si(111) surface. Calculated minimum energy pathways (MEPs) indicate that the chain mechanism is viable in the case of C(2)H(2), whereas for C(2)H(4) the stabilization of the intermediate state is so small and the barrier for H-abstraction so (relatively) large that the molecule is more likely to desorb than to form a stable adsorbed species. For phenylacetylene and styrene, stabilization of the intermediate state and decrease of the H-abstraction barrier take place. While a stable adsorbed species exists in both cases, the overall heat of adsorption is larger for the alkyne molecules.

First-principles string molecular dynamics: An efficient approach for finding chemical reaction pathways

A recently proposed approach, called "string method," allows us to find minimum energy pathways connecting two metastable states of a system [W. E et al., Phys. Rev. B 66, 052301 (2002)]. So far this approach has been only used with empirical force field parametrizations of the atomic potential energy surface or in the context of macroscopic continuum models. Here we show that the string method can be efficiently combined with first-principles molecular dynamics to provide an accurate description of chemical reaction pathways and barriers. We illustrate the first-principles string molecular dynamics by applying it to the study of a surface chemical reaction, for which extensive experimental and theoretical works are available, namely, the adsorption of H(2) on the reconstructed Si(100) surface.