Elucidating the vibrational spectra of hydrogen-bonded aggregates in solution: electronic structure calculations with implicit solvent and first-principles molecular dynamics simulations with explicit solvent for 1-hexanol in n-hexane

Fourier transform infrared spectroscopy is a popular method for the experimental investigation of hydrogen-bonded aggregates, but linking spectral information to microscopic information on aggregate size distribution and aggregate architecture is an arduous task. Static electronic structure calculations with an implicit solvent model, Car-Parrinello molecular dynamics (CPMD) using the Becke-Lee-Yang-Parr (BLYP) exchange and correlation energy functionals and classical molecular dynamics simulations for the all-atom version of the optimized parameters for liquid simulations (OPLS-AA) force field were carried out for an ensemble of 1-hexanol aggregates solvated in n-hexane. The initial configurations for these calculations were size-selected from a distribution of aggregates obtained from a large-scale Monte Carlo simulation. The vibrational spectra computed from the static electronic structure calculations for monomers and dimers and from the CPMD simulations for aggregates up to pentamers demonstrate the extent of the contribution of dangling or nondonating hydroxyl groups found in linear and branched aggregates to the "monomeric" peak. Furthermore, the computed spectra show that there is no simple relationship between peak shift and aggregate size nor architecture, but the effect of hydrogen-bond cooperativity is shown to differentiate polymer-like (cooperative) and dimer-like (noncooperative) hydrogen bonds in the vibrational spectrum. In contrast to the static electronic structure calculations and the CPMD simulations, the classical molecular dynamics simulations greatly underestimate the vibrational peak shift due to hydrogen bonding.     

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).     

Ion-orbital coupling in Car-Parrinello calculations of hydrogen-bond vibrational dynamics: case study with the NH3-HCl dimer

We have performed Car-Parrinello molecular dynamics (CPMD) calculations of the hydrogen-bonded NH(3)-HCl dimer. Our main aim is to establish how ionic-orbital coupling in CPMD affects the vibrational dynamics in hydrogen-bonded systems by characterizing the dependence of the calculated vibrational frequencies upon the orbital mass in the adiabatic limit of Car-Parrinello calculations. We use the example of the NH(3)-HCl dimer because of interest in its vibrational spectrum, in particular the magnitude of the frequency shift of the H-Cl stretch due to the anharmonic interactions when the hydrogen bond is formed. We find that an orbital mass of about 100 a.u. or smaller is required in order for the ion-orbital coupling to be linear in orbital mass, and the results for which can be accurately extrapolated to the adiabatic limit of zero orbital mass. We argue that this is general for hydrogen-bonded systems, suggesting that typical orbital mass values used in CPMD are too high to accurately describe vibrational dynamics in hydrogen-bonded systems. Our results also show that the usual application of a scaling factor to the CPMD frequencies to correct for the effects of orbital mass is not valid. For the dynamics of the dimer, we find that the H-Cl stretch and the N-H-Cl bend are significantly coupled, suggesting that it is important to include the latter degree of freedom in quantum dynamical calculations. Results from our calculations with deuterium-substitution show that both these degrees of freedom have significant anharmonic interactions. Our calculated frequency for the H-Cl stretch using the Becke-exchange Lee-Yang-Parr correlation functional compares reasonably well with a previous second-order M?ller-Plesset calculation with anharmonic corrections, although it is low compared to the experimental value for the dimer trapped in a neon-matrix.     

Using force-matching to reveal essential differences between density functionals in ab initio molecular dynamics simulations

The exchange-correlation (XC) functional and value of the electronic fictitious mass μ can be two major sources of systematic errors in ab initio Car-Parrinello Molecular Dynamics (CPMD) simulations, and have a significant impact on the structural and dynamic properties of condensed-phase systems. In this work, an attempt is made to identify the origin of differences in liquid water properties generated from CPMD simulations run with the BLYP and HCTH∕120 XC functionals and two different values of μ (representative of "small" and "large" limits) by analyzing the effective pairwise atom-atom interactions. The force-matching (FM) algorithm is used to map CPMD interactions into non-polarizable, empirical potentials defined by bonded interactions, pairwise short-ranged interactions in numerical form, and Coulombic interactions via atomic partial charges. The effective interaction models are derived for the BLYP XC functional with μ=340 a.u. and μ=1100 a.u. (BLYP-340 and BLYP-1100 simulations) and the HCTH∕120 XC functional with μ=340 a.u. (HCTH-340 simulation). The BLYP-340 simulation results in overstructured water with slow dynamics. In contrast, the BLYP-1100 and HCTH-340 simulations both produce radial distribution functions (indicative of structure) that are in reasonably good agreement with experiment. It is shown that the main difference between the BLYP-340 and HCTH-340 effective potentials arises in the short-ranged nonbonded interactions (in hydrogen bonding regions), while the difference between the BLYP-340 and BLYP-1100 interactions is mainly in the long-ranged electrostatic components. Collectively, these results demonstrate how the FM method can be used to further characterize various simulation ensembles (e.g., density-functional theory via CPMD). An analytical representation of each effective interaction water model, which is easy to implement, is presented.     

Theoretical calculations: Can Gibbs free energy for intermolecular complexes be predicted efficiently and accurately?

The theoretical study has been performed to refine the procedure for calculations of Gibbs free energy with a relative accuracy of less than 1 kcal/mol. Three benchmark intermolecular complexes are examined via several quantum-chemical methods, including the second-order Moller-Plesset perturbation (MP2), coupled cluster (CCSD(T)), and density functional (BLYP, B3LYP) theories augmented by Dunnings correlation-consistent basis sets. The effects of electron correlation, basis set size, and anharmonicity are systematically analyzed, and the results are compared with available experimental data. The results of the calculations suggest that experimental accuracy can be reached only by extrapolation of MP2 and CCSD(T) total energies to the complete basis set. The contribution of anharmonicity to the zero point energy and TDeltaSint values is fairly small. The new, economic way to reach chemical accuracy in the calculations of the thermodynamic parameters of intermolecular interactions is proposed. In addition, interaction energy (De) and free energy change (DeltaA) for considered species have been evaluated by Carr-Parrinello molecular dynamics (CPMD) simulations and static BLYP-plane wave calculations. The free energy change along the reaction paths were determined by the thermodynamic integration/"Blue Moon Ensemble" technique. Comparison between obtained values, and available experimental and conventional ab initio results has been made. We found that the accuracy of CPMD simulations is affected by several factors, including statistical uncertainty and convergence of constrained forces (TD integration), and the nature of DFT (density functional theory) functional. The results show that CPMD technique is capable of reproducing interaction and free energy with an accuracy of 1 kcal/mol and 2-3 kcal/mol respectively.     

Lagrangian formulation with dissipation of Born-Oppenheimer molecular dynamics using the density-functional tight-binding method

An important element determining the time requirements of Born-Oppenheimer molecular dynamics (BOMD) is the convergence rate of the self-consistent solution of Roothaan equations (SCF). We show here that improved convergence and dynamics stability can be achieved by use of a Lagrangian formalism of BOMD with dissipation (DXL-BOMD). In the DXL-BOMD algorithm, an auxiliary electronic variable (e.g., the electron density or Fock matrix) is propagated and a dissipative force is added in the propagation to maintain the stability of the dynamics. Implementation of the approach in the self-consistent charge density functional tight-binding method makes possible simulations that are several hundred picoseconds in lengths, in contrast to earlier DFT-based BOMD calculations, which have been limited to tens of picoseconds or less. The increase in the simulation time results in a more meaningful evaluation of the DXL-BOMD method. A comparison is made of the number of iterations (and time) required for convergence of the SCF with DXL-BOMD and a standard method (starting with a zero charge guess for all atoms at each step), which gives accurate propagation with reasonable SCF convergence criteria. From tests using NVE simulations of C(2)F(4) and 20 neutral amino acid molecules in the gas phase, it is found that DXL-BOMD can improve SCF convergence by up to a factor of two over the standard method. Corresponding results are obtained in simulations of 32 water molecules in a periodic box. Linear response theory is used to analyze the relationship between the energy drift and the correlation of geometry propagation errors.     

Identification of active sites of biomolecules. 1. Methyl-alpha-mannopyranoside and Fe(III)

Car-Parrinello molecular dynamics (CPMD) simulations, DFT chemical reactivity index calculations, and mass spectrometric measurements are combined in an integrated effort to elucidate the details of the coordination of a transition-metal ion to a carbohydrate. The impact of the interaction with the FeIII ion on the glycosidic linkage conformation of methyl-alpha-d-mannopyranoside is studied by classical molecular dynamics (MD) and CPMD simulations. This study shows that FeIII interacts with specific hydroxyl oxygen atoms of the carbohydrate, affecting the ground state carbohydrate conformation. These conformational details are discussed in terms of a set of supporting experiments involving electrospray ionization mass spectrometry, and CPMD simulations clearly indicate that the specific conformational preference is due to intramolecular hydrogen bonding. Classical MD simulations proved insensitive to these important chemical properties. Thus, we demonstrate the importance of chemical reactivity calculations and CPMD simulations in predicting the active sites of biological molecules toward metal cations.     

Ab initio molecular dynamics study of 45S5 bioactive silicate glass

Bioglass 45S5, the prototype of bioactive melt-quenched silicate glasses, was modeled by means of Car-Parrinello molecular dynamics (CPMD) simulations. Although long-range structural properties cannot be modeled by using this ab initio approach, the accuracy of CPMD simulations is exploited here to provide insight into the short-range structure and to analyze vibrational and electronic properties of this biomaterial. Detailed structural analysis in the short-range scale provided insight into the local environment of modifier Na and Ca ions: a possible key role of these cations in organizing the glass network by connecting different chains and fragments into specific, rather flexible geometries was proposed. The individual contributions of different species to the vibrational density of states were separated and discussed, allowing the identification of specific features in the vibrational spectrum, such as those related to phosphate groups. The components of the electronic density of states were also analyzed, enabling us to identify correlations between the electronic structure and the structural properties, such as the different bonding character of Si-O bonds involving bridging or nonbridging oxygen atoms.     

Direct dynamics study of ultrafast vibrational energy relaxation in ice Ih

Car-Parrinello molecular dynamics (CPMD) and a previously developed wave packet model are used to study ultrafast relaxation in water clusters. Water clusters of 15 water molecules are used to represent ice Ih. The relaxation is studied by exciting a symmetric or an asymmetric stretch mode of the central water molecule. The CPMD results suggest that relaxation occurs within 100 fs. This is in agreement with experimental work by Woutersen and Bakker and the earlier wave packet calculations. The CPMD results further indicate that the excitation energy is transferred both intramolecularly and intermolecularly on roughly the same time scale. The intramolecular energy transfer occurs predominantly between the symmetric and asymmetric modes while the bend mode is largely left unexcited on the short time scale studied here.     

Optimization of auxiliary basis sets for the LEDO expansion and a projection technique for LEDO-DFT

We present a systematic procedure for the optimization of the expansion basis for the limited expansion of diatomic overlap density functional theory (LEDO-DFT) and report on optimized auxiliary orbitals for the Ahlrichs split valence plus polarization basis set (SVP) for the elements H, Li--F, and Na--Cl. A new method to deal with near-linear dependences in the LEDO expansion basis is introduced, which greatly reduces the computational effort of LEDO-DFT calculations. Numerical results for a test set of small molecules demonstrate the accuracy of electronic energies, structural parameters, dipole moments, and harmonic frequencies. For larger molecular systems the numerical errors introduced by the LEDO approximation can lead to an uncontrollable behavior of the self-consistent field (SCF) process. A projection technique suggested by L?wdin is presented in the framework of LEDO-DFT, which guarantees for SCF convergence. Numerical results on some critical test molecules suggest the general applicability of the auxiliary orbitals presented in combination with this projection technique. Timing results indicate that LEDO-DFT is competitive with conventional density fitting methods.     

Excited state geometries within time-dependent and restricted open-shell density functional theories

Singlet excited state geometries of a set of medium sized molecules with different characteristic lowest excitations are studied. Geometry optimizations of excited states are performed with two closely related restricted open-shell Kohn–Sham methods and within linear response to time-dependent density functional theory. The results are compared to wave-function based methods. Excitation energies (vertical and adiabatic) calculated from the open-shell methods show systematic errors depending on the type of excitation. However, for all states accessible by the restricted methods a good agreement for the geometries with time-dependent density functional theory and wave-function based methods is found. An analysis of the energy with respect to the mixing angle for the singly occupied orbitals reveals that some states (mostly [n→π^*]) are stable when symmetry constraints are relaxed and others (mostly [π→π^*]) are instable. This has major implications on the applicability of the restricted open-shell methods in molecular dynamics simulations.     

Valence-bond calculations on ZNO and HGO using integrals computed through the semiempirical AM1 method

Valence-bond calculations have been carried out on ZnO and HgO using a basis set of Slatertype atomic orbitals and the one- and two-electron integrals as computed in the semiempirical AM 1 molecular orbital method. The zero differential overlap approximation has been used to calculate integrals between atomic orbital Slater determinants using the rules for matrix elements between determinants formed by orthogonal orbitals. Diabatic and adiabatic curves have been analyzed for the two systems, and results compared with molecular orbital AM 1 results. © 1992 John Wiley & Sons, Inc.     

An overlap criterion for selection of core orbitals

 An overlap criterion is defined that connects the identification of core orbitals in a molecular system, which can be problematic, to that in isolated atoms, which is well defined. This approach has been tested on a variety of troublesome systems that have been identified in the literature, including molecules containing third-row main-group elements, and is shown to remove errors of up to 100 kcal/mol arising from an inconsistent treatment of core orbitals at different locations on a potential-energy surface. For some systems and choices of core orbitals, errors as large as 19 kcal/mol can be introduced even when consistent sets of orbitals are frozen, and the new method is shown to identify these cases of substantial core–valence mixing. Finally, even when there is limited core–valence mixing, the frozen-core approximation can introduce errors of more than 5 kcal/mol, which is much larger than the presumed accuracy of models such as G2 and CBS-QB3. The source of these errors includes interatomic core–core and core–valence dispersion forces. Received: 31 August 2001 / Accepted: 11 October 2001 / Published online: 9 January 2002     

Preferred conformation of the glycosidic linkage of methyl-beta-mannose

The conformational preference of the glycosidic linkage of methyl-beta-mannose was studied in the gas phase and in aqueous solution by ab initio calculations, and by molecular dynamics (MD) and Car-Parrinello molecular dynamics (CPMD) simulations. MD simulations were performed with various water potential functions to study the impact of the chosen water potential on the predicted conformational preference of the glycosidic linkage of the carbohydrate in solution. This study shows that the trans (t) orientation of the glycosidic linkage of methyl-beta-mannose is preferred over its gauche clockwise (g+) orientation in solution. CPMD simulations clearly indicate that this preference is due to intermolecular hydrogen bonding with surrounding water molecules, whereas no such information could be demonstrated by MD simulations. This study demonstrates the importance of ab initio molecular dynamics simulations in studying the structural properties of carbohydrate-water interactions.     

How does ammonium dynamically interact with benzene in aqueous media? A first principle study using the Car-Parrinello molecular dynamics method

The Car-Parrinello molecular dynamics (CPMD) method was used to study the dynamic characteristics of the cation-pi interaction between ammonium and benzene in gaseous and aqueous media. The results obtained from the CPMD calculation on the cation-pi complex in the gaseous state were very similar to those calculated from the Gaussian98 program with DFT and MP2 algorithms, demonstrating that CPMD is a valid approach for studying this system. Unlike the interaction in the gaseous state, our 12-ps CPMD simulation showed that the geometry of the complex in aqueous solution changes frequently in terms of the interaction angles and distances. Furthermore, the simulation revealed that the ammonium is constantly oscillating above the benzene plane in an aqueous environment and interacts with benzene mostly through three of its hydrogen atoms. In contrast, the interaction of the cation with the aromatic molecule in the gaseous state involves two hydrogen atoms. In addition, the free energy profile in aqueous solution was studied using constrained CPMD simulations, resulting in a calculated binding free energy of -5.75 kcal/mol at an optimum interaction distance of approximately 3.25 A, indicating that the cation-pi interaction between ammonium and benzene is stable even in aqueous solution. Thus, this CPMD study suggested that the cation-pi interaction between an ammonium (group) and an aromatic structure could take place even on surfaces of protein or nucleic acids in solution.     

Rupture force analysis and the associated systematic errors in force spectroscopy by AFM

Force spectroscopy is a new and valuable tool in physical chemistry and biophysics. However, data analysis has yet to be standardized, hindering the advancement of the technique. In this article, treatment of the rupture forces is described in the framework of the Bell-Evans model, and the systematic errors associated with the tether effect for approaches that utilize the most probable, the median, and the mean rupture forces are compared. It is shown that significant systematic errors in the dissociation rate can result from nonlinear loading with polymeric tethers even if the apparent loading rate is used in the analysis. Analytical expressions for the systematic errors are provided for the most probable and median forces. The use of these expressions to correct the associated systematic errors is illustrated by the analysis of the measured rupture forces between single hexadecane molecules in water. It is noted that the measured distributions of rupture forces often contain high forces that are unaccounted for by theoretical models. Experimental data indicate that the most significant effect of the high forces "tail" is on the dissociation rate obtained from the median force analysis whereas the barrier width appears to be unaffected.     

Mechanisms of proton transfer in Nafion: elementary reactions at the sulfonic acid groups

Proton transfer reactions at the sulfonic acid groups in Nafion were theoretically studied, using complexes formed from triflic acid (CF3SO3H), H3O+ and H2O, as model systems. The investigations began with searching for potential precursors and transition states at low hydration levels, using the test-particle model (T-model), density functional theory (DFT) and ab initio calculations. They were employed as starting configurations in Born-Oppenheimer molecular dynamics (BOMD) simulations at 298 K, from which elementary reactions were analyzed and categorized. For the H3O+-H2O complexes, BOMD simulations suggested that a quasi-dynamic equilibrium could be established between the Eigen and Zundel complexes, and that was considered to be one of the most important elementary reactions in the proton transfer process. The average lifetime of H3O+ obtained from BOMD simulations is close to the lowest limit, estimated from low-frequency vibrational spectroscopy. It was demonstrated that proton transfer reactions at -SO3H are not concerted, due to the thermal energy fluctuation and the existence of various quasi-dynamic equilibria, and -SO3H could directly and indirectly mediate proton transfer reactions through the formation of proton defects, as well as the -SO3- and -SO3H2+ transition states.     

Ab initio molecular dynamics of liquid 1,3-dimethylimidazolium chloride

Density-functional-based Car-Parrinello molecular dynamics (CPMD) simulations have been performed for the ionic liquid 1,3-dimethylimidazolium chloride, [dmim]Cl, at 438 K. The local structure of the liquid is described in terms of various partial radial distribution functions and anisotropic spatial distributions, which reveal a significant extent of hydrogen bonding. The cation-anion distribution simulated with the BP86 functional is in qualitative agreement with the structural model derived from neutron diffraction data for the liquid, whereas the theoretical cation-cation distribution shows less satisfactory accord. Population analyses indicate noticeable charge transfer from anions to cations, and specific CH...Cl hydrogen bonds are characterized in terms of donor-acceptor interactions between lone pairs on Cl and antibonding sigma(CH) orbitals.     

Charge controlled changes in the cluster and spin dynamics of Sc3N@C80(CF3)2: the flexible spin density distribution and its impact on ESR spectra

Cluster and spin dynamics of a Sc(3)N@C(80)(CF(3))(2) derivative are studied by DFT in different charge states, from -3 to +1. For the neutral Sc(3)N@C(80)(CF(3))(2), static DFT computations of many cluster conformers as well as Born-Oppenheimer molecular dynamics (BOMD) show that addition of two CF(3) groups to Sc(3)N@C(80) significantly changes dynamics of the Sc(3)N cluster: instead of free rotation as in Sc(3)N@C(80), the cluster in Sc(3)N@C(80)(CF(3))(2) exhibits only hindered motions. Similar cluster dynamics is found in the mono- and trianions of Sc(3)N@C(80)(CF(3))(2), while free rotation of the cluster is found in the cation. In the radical species, motions of the cluster dramatically change spin-density distribution. Spin populations of the metal atoms and the carbon cage are followed along the BOMD trajectories to reveal the details of the spin-flow. (45)Sc ESR hyperfine coupling constants integrated over BOMD trajectories are found to be substantially different from the results of static DFT computations, which emphasizes that cluster dynamics should be taken into account for reliable predictions of spectroscopic properties.