Comparing van der Waals DFT methods for water on NaCl(001) and MgO(001).

In this work, a range of van der Waals type density functionals are applied to the H2O/NaCl(001) and H2O/MgO(001) interface systems to explore the effect of an explicit dispersion treatment. The functionals we use are the self-consistent vdW functionals vdW-DF, vdW-DF2, optPBE-vdW, optB88-vdW, optB86b-vdW, and vdW-DF-cx, as well as the dispersion-corrected PBE-TS and PBE-D2 methods; they are all compared with the standard PBE functional. For both NaCl(001) and MgO(001), we find that the dispersion-flavoured functionals stabilize the water-surface interface by approximately 20%-40% compared to the PBE results. For NaCl(001), where the water molecules remain intact for all overlayers, the dominant contribution to the adsorption energy from "density functional theory dispersion" stems from the water-surface interactions rather than the water-water interactions. The optPBE-vdW and vdW-DF-cx functionals yield adsorption energies in good agreement with available experimental values for both NaCl and MgO. To probe the strengths of the perturbations of the adsorbed water molecules, we also calculated water dipole moments and found an increase up to 85% for water at the MgO(001) surface and 70% at the NaCl(001) surface, compared to the gas-phase dipole moment.

Treatment of delocalized electron transfer in periodic and embedded cluster DFT calculations: The case of Cu on ZnO (10(1)0).

We assess the consequences of the interface model-embedded-cluster or periodic-slab model-on the ability of DFT calculations to describe charge transfer (CT) in a particularly challenging case where periodic-slab calculations indicate a delocalized charge-transfer state. Our example is Cu atom adsorption on ZnO(10(1)0), and in fact the periodic slab calculations indicate three types of CT depending on the adsorption site: full CT, partial CT, and no CT. Interestingly, when full CT occurs in the periodic calculations, the calculated Cu atom adsorption energy depends on the underlying ZnO substrate supercell size, since when the electron enters the ZnO it delocalizes over as many atoms as possible. In the embedded-cluster calculations, the electron transferred to the ZnO delocalizes over the entire cluster region, and as a result the calculated Cu atom adsorption energy does not agree with the value obtained using a large periodic supercell, but instead to the adsorption energy obtained for a periodic supercell of roughly the same size as the embedded cluster. Different density functionals (of GGA and hybrid types) and basis sets (local atom-centered and plane-waves) were assessed, and we show that embedded clusters can be used to model Cu adsorption on ZnO(10(1)0), as long as care is taken to account for the effects of CT.

Calculation of the stability of nonperiodic solids using classical force fields and the method of increments: N2O as an example.

Combining classical force fields for the Hartree-Fock (HF) part and the method of increments for post-HF contributions, we calculate the cohesive energy of the ordered and randomly disordered nitrous oxide (N2 O) solid. At 0 K, ordered N2 O is most favorable with a cohesive energy of -27.7 kJ/mol. At temperatures above 60 K, more disordered structures become compatible and a phase transition to completely disordered N2 O is predicted. Comparison with experiment in literature suggests that experimentally prepared N2 O crystals are mainly disordered due to a prohibitively high activation energy of ordering processes. © 2015 Wiley Periodicals, Inc.

An efficient and extensible format, library, and API for binary trajectory data from molecular simulations.

Molecular dynamics simulations is an important application in theoretical chemistry, and with the large high-performance computing resources available today the programs also generate huge amounts of output data. In particular in life sciences, with complex biomolecules such as proteins, simulation projects regularly deal with several terabytes of data. Apart from the need for more cost-efficient storage, it is increasingly important to be able to archive data, secure the integrity against disk or file transfer errors, to provide rapid access, and facilitate exchange of data through open interfaces. There is already a whole range of different formats used, but few if any of them (including our previous ones) fulfill all these goals. To address these shortcomings, we present "Trajectory Next Generation" (TNG)--a flexible but highly optimized and efficient file format designed with interoperability in mind. TNG both provides state-of-the-art multiframe compression as well as a container framework that will make it possible to extend it with new compression algorithms without modifications in programs using it. TNG will be the new file format in the next major release of the GROMACS package, but it has been implemented as a separate library and API with liberal licensing to enable wide adoption both in academic and commercial codes.

Band-Filling Correction Method for Accurate Adsorption Energy Calculations: A Cu/ZnO Case Study.

We present a simple method, the "band-filling correction", to calculate accurate adsorption energies (Eads) in the low coverage limit from finite-size supercell slab calculations using DFT. We show that it is necessary to use such a correction if charge transfer takes place between the adsorbate and the substrate, resulting in the substrate bands either filling up or becoming depleted. With this correction scheme, we calculate Eads of an isolated Cu atom adsorbed on the ZnO(101̅0) surface. Without the correction, the calculated Eads is highly coverage-dependent, even for surface supercells that would typically be considered very large (in the range from 1 nm × 1 nm to 2.5 nm × 2.5 nm). The correction scheme works very well for semilocal functionals, where the corrected Eads is converged within 0.01 eV for all coverages. The correction scheme also works well for hybrid functionals if a large supercell is used and the exact exchange interaction is screened.

Fitting properties from density functional theory based molecular dynamics simulations to parameterize a rigid water force field.

In the quest towards coarse-grained potentials and new water models, we present an extension of the force matching technique to parameterize an all-atom force field for rigid water. The methodology presented here allows to improve the matching procedure by first optimizing the weighting exponents present in the objective function. A new gauge for unambiguously evaluating the quality of the fit has been introduced; it is based on the root mean square difference of the distributions of target properties between reference data and fitted potentials. Four rigid water models have been parameterized; the matching procedure has been used to assess the role of the ghost atom in TIP4P-like models and of electrostatic damping. In the former case, burying the negative charge inside the molecule allows to fit better the torques. In the latter, since short-range interactions are damped, a better fit of the forces is obtained. Overall, the best performing model is the one with a ghost atom and with electrostatic damping. The approach shown in this paper is of general validity and could be applied to any matching algorithm and to any level of coarse graining, also for non-rigid molecules.

Cations strongly reduce electron-hopping rates in aqueous solutions.

We study how the ultrafast intermolecular hopping of electrons excited from the water O1s core level into unoccupied orbitals depends on the local molecular environment in liquid water. Our probe is the resonant Auger decay of the water O1s core hole (lifetime ∼3.6 fs), by which we show that the electron-hopping rate can be significantly reduced when a first-shell water molecule is replaced by an atomic ion. Decays resulting from excitations at the O1s post-edge feature (∼540 eV) of 6 m LiBr and 3 m MgBr(2) aqueous solutions reveal electron-hopping times of ∼1.5 and 1.9 fs, respectively; the latter represents a 4-fold increase compared to the corresponding value in neat water. The slower electron-hopping in electrolytes, which shows a strong dependence on the charge of the cations, can be explained by ion-induced reduction of water-water orbital mixing. Density functional theory electronic structure calculations of solvation geometries obtained from molecular dynamics simulations reveal that this phenomenon largely arises from electrostatic perturbations of the solvating water molecules by the solvated ions. Our results demonstrate that it is possible to deliberately manipulate the rate of charge transfer via electron-hopping in aqueous media.

Structure and dynamics of water dangling OH bonds in hydrophobic hydration shells. Comparison of simulation and experiment.

Molecular dynamics and electric field strength simulations are performed in order to quantify the structural, dynamic, and vibrational properties of non-H-bonded (dangling) OH groups in the hydration shell of neopentane, as well as in bulk water. The results are found to be in good agreement with the experimentally observed high-frequency (∼3660 cm(-1)) OH band arising from the hydration shell of neopentanol dissolved in HOD/D(2)O, obtained by analyzing variable concentration Raman spectra using multivariate curve resolution (Raman-MCR). The simulation results further indicate that hydration shell dangling OH groups preferentially point toward the central carbon atom of neopentane to a degree that increases with the lifetime of the dangling OH.

On the origins of core-electron chemical shifts of small biomolecules in aqueous solution: insights from photoemission and ab initio calculations of glycine(aq).

The local electronic structure of glycine in neutral, basic, and acidic aqueous solution is studied experimentally by X-ray photoelectron spectroscopy and theoretically by molecular dynamics simulations accompanied by first-principle electronic structure and spectrum calculations. Measured and computed nitrogen and carbon 1s binding energies are assigned to different local atomic environments, which are shown to be sensitive to the protonation/deprotonation of the amino and carboxyl functional groups at different pH values. We report the first accurate computation of core-level chemical shifts of an aqueous solute in various protonation states and explicitly show how the distributions of photoelectron binding energies (core-level peak widths) are related to the details of the hydrogen bond configurations, i.e. the geometries of the water solvation shell and the associated electronic screening. The comparison between the experiments and calculations further enables the separation of protonation-induced (covalent) and solvent-induced (electrostatic) screening contributions to the chemical shifts in the aqueous phase. The present core-level line shape analysis facilitates an accurate interpretation of photoelectron spectra from larger biomolecular solutes than glycine.

Trajectory NG: portable, compressed, general molecular dynamics trajectories.

We present general algorithms for the compression of molecular dynamics trajectories. The standard ways to store MD trajectories as text or as raw binary floating point numbers result in very large files when efficient simulation programs are used on supercomputers. Our algorithms are based on the observation that differences in atomic coordinates/velocities, in either time or space, are generally smaller than the absolute values of the coordinates/velocities. Also, it is often possible to store values at a lower precision. We apply several compression schemes to compress the resulting differences further. The most efficient algorithms developed here use a block sorting algorithm in combination with Huffman coding. Depending on the frequency of storage of frames in the trajectory, either space, time, or combinations of space and time differences are usually the most efficient. We compare the efficiency of our algorithms with each other and with other algorithms present in the literature for various systems: liquid argon, water, a virus capsid solvated in 15 mM aqueous NaCl, and solid magnesium oxide. We perform tests to determine how much precision is necessary to obtain accurate structural and dynamic properties, as well as benchmark a parallelized implementation of the algorithms. We obtain compression ratios (compared to single precision floating point) of 1:3.3-1:35 depending on the frequency of storage of frames and the system studied.

Al3+, Ca2+, Mg2+, and Li+ in aqueous solution: calculated first-shell anharmonic OH vibrations at 300 K.

The anharmonic OH stretching vibrational frequencies, ν(OH), for the first-shell water molecules around the Li(+), Ca(2+), Mg(2+), and Al(3+) ions in dilute aqueous solutions have been calculated based on classical molecular dynamics (MD) simulations and quantum-mechanical (QM) calculations. For Li(+)(aq), Ca(2+)(aq), Mg(2+)(aq), and Al(3+)(aq), our calculated IR frequency shifts, Δν(OH), with respect to the gas-phase water frequency, are about -300, -350, -450, and -750 cm(-1), compared to -290, -290, -420, and -830 cm(-1) from experimental infrared (IR) studies. The agreement is thus quite good, except for the order between Li(+) and Ca(2+). Given that the polarizing field from the Ca(2+) ion ought to be larger than that from Li(+)(aq), our calculated result seems reasonable. Also the absolute OH frequencies agree well with experiment. The method we used is a sequential four-step procedure: QM(electronic) to make a force field+MD simulation+QM(electronic) for point-charge-embedded M(n+) (H(2)O)(y) (second shell) (H(2)O)(z) (third shell) clusters+QM(vibrational) to yield the OH spectrum. The many-body Ca(2+)-water force-field presented in this paper is new. IR intensity-weighting of the density-of-states frequency distributions was carried out by means of the squared dipole moment derivatives.

Sequential penning ionization: harvesting energy with ions.

We report the observation of the ejection of electrons caused by collisions of excited atoms with ions, rather than neutrals, leading to the production of doubly charged ions. Doping superfluid He droplets with methyl iodide and exposing them to electrons enhances the formation of doubly charged iodine atoms at the threshold for the production of two metastable He atoms. These observations point toward a novel ionization process where doubly charged ions are produced by sequential Penning ionization. In some cases, depending on the neutral target, the process also leads to a subsequent Coulomb explosion of the dopant.

OD vibrations and hydration structure in an Al3+ (aq) solution from a Car-Parrinello molecular-dynamics simulation.

The optimized geometry, energetics, and vibrational properties of Al(D2O)n3+ clusters, with n = 1,2,4, and 6, have been studied using plane waves, different local basis sets, different methodologies [density-functional theory, MP2, CCSD(T)], and different functionals (BLYP, PBE). Moreover, Car-Parrinello molecular-dynamics (MD) simulations using the BLYP functional, plane waves, and the Vanderbilt ultrasoft pseudopotentials have been performed for an aqueous Al3+ solution with 1 ion and 32 D2O molecules in a periodic box at room temperature, studied for 10 ps. The cluster calculations were performed to pinpoint possible shortcomings of the electronic structure description used in the Car-Parinello MD (CPMD) simulation. For the clusters, the hydration structure and interaction energies calculated with the 'BLYP/plane-wave' approach agree well with high-level ab initio methods but the exchange-correlation functional introduces errors in the OD stretching frequencies (both in the absolute values and in the ion-induced shifts). For the aqueous solution, the CPMD simulation yields structural properties in good agreement with experimental data. The CPMD-simulated OD stretching vibrational band for the first-shell water molecules around Al3+ is strongly downshifted by the influence of the ion and is compared with experimental data from the literature. To make such a comparison meaningful, the influences of a number of systematic effects have been addressed, such as the exchange-correlation functional, the fictitious electron mass, anharmonicity effects, and the small box size in the simulation. Each of these factors (except the last one) is found to affect the OD frequency by 100 cm(-1) or more. The final "corrected" frequencies agree with experiment within approximately 30 cm(-1) for bulk water but are too little downshifted for the first-shell Al3+ (aq) water molecules (by approximately 200 cm(-1)).

Distorted five-fold coordination of Cu2+(aq) from a Car-Parrinello molecular dynamics simulation.

The solvation shell structure and dynamics of a single Cu2+ ion in a periodic box with 32 water molecules under ambient conditions has been investigated using Car-Parrinello molecular dynamics simulations in a time-window of 18 ps. Five-fold coordination with four equidistant equatorial water molecules at 2.00 A and one axial water molecule at 2.45 A from the Cu2+ ion is found. A "hole" without water molecules is found on the opposite side of the axial water. The ion-water bonding character for the equatorial water molecules is different from that of the axial water molecules, as shown by a localized orbital analysis of the electronic structure. Moreover, the calculated OD stretching vibrational band for the equatorial water molecules lies ca. 175 cm-1 below the axial-water band, in good agreement with experimental data. The equatorial-water band lies below, and the axial-water band above, the pure liquid D2O band, also in agreement with experimental data.

Using MD snapshots in ab initio and DFT calculations: OH vibrations in the first hydration shell around Li+(aq).

The average OH stretching vibrational frequency for the water molecules in the first hydration shell around a Li(+) ion in a dilute aqueous solution was calculated by a hybrid molecular dynamics + quantum-mechanical ("MD + QM") approach. Using geometry configurations from a series of snapshots from an MD simulation, the anharmonic, uncoupled OH stretching frequencies were calculated for 100 first-shell OH oscillators at the B3LYP and HF/6-31G(d,p) levels of theory, explicitly including the first shell and the relevant second shell water molecules into charge-embedded supermolecular QM calculations. Infrared intensity-weighting of the density-of-states (DOS) distributions by means of the squared dipole moment derivatives (which vary by a factor of 20 over the OH stretching frequency band at the B3LYP level), changes the downshift from approximately -205 to -275 cm(-1) at the B3LYP level. Explicit inclusion of relevant third-shell water molecules in the supermolecular cluster leads to a further downshift by approximately -30 cm(-1). Our final estimated average downshift is approximately -305 cm(-1). The experimental value lies somewhere in the range between -290 and -420 cm(-1). Also, the absolute nu(OH) frequency is well reproduced in our calculations. "In-liquid" instantaneous correlation curves between nu(OH) and various typical H-bond strength parameters such as R(O...O), R(H...O), the intramolecular OH bond length, and the IR intensity are presented. Some of these correlations are robust and persist also for the rather distorted instantaneous geometries in the liquid; others are less so.

Car-Parrinello molecular dynamics simulation of Fe 3+ (aq).

The optimized geometry and energetic properties of Fe(D2O)n 3+ clusters, with n = 4 and 6, have been studied with density-functional theory calculations and the BLYP functional, and the hydration of a single Fe 3+ ion in a periodic box with 32 water molecules at room temperature has been studied with Car-Parrinello molecular dynamics and the same functional. We have compared the results from the CPMD simulation with classical MD simulations, using a flexible SPC-based water model and the same number of water molecules, to evaluate the relative strengths and weaknesses of the two MD methods. The classical MD simulations and the CPMD simulations both give Fe-water distances in good agreement with experiment, but for the intramolecular vibrations, the classical MD yields considerably better absolute frequencies and ion-induced frequency shifts. On the other hand, the CPMD method performs considerably better than the classical MD in describing the intramolecular geometry of the water molecule in the first hydration shell and the average first shell...second shell hydrogen-bond distance. Differences between the two methods are also found with respect to the second-shell water orientations. The effect of the small box size (32 vs 512 water molecules) was evaluated by comparing results from classical simulations using different box sizes; non-negligible effects are found for the ion-water distance and the tilt angles of the water molecules in the second hydration shell and for the O-D stretching vibrational frequencies of the water molecules in the first hydration shell.

Many-body potentials for aqueous Li(+), Na(+), Mg(2+), and Al(3+): comparison of effective three-body potentials and polarizable models.

Many-body potentials for the aqueous Li(+), Na(+), Mg(2+), and Al(3+) ions have been constructed from ab initio cluster calculations. Pure pair, effective pair, effective three-body, and effective polarizable models were created and used in subsequent molecular dynamics simulations. The structures of the first and second solvation shells were studied using radial distribution functions and angular-radial distribution functions. The effective three-body and polarizable potentials yield similar first-shell structures, while the contraction of the O-O distances between the first and second solvation shells is more pronounced with the polarizable potentials. The definition of the tilt angle of the water molecules around the ions is discussed. When a proper definition is used, it is found that for Li(+), Mg(2+), and Al(3+) the water molecules prefer a trigonal orientation, but for Na(+) a tetrahedral orientation (ion in lone-pair direction) is preferred. The self-diffusion coefficients for the water molecules and the ions were calculated; the ionic values follow the order obtained from experiment, although the simulated absolute values are smaller than experiment for Mg(2+) and Al(3+).