Mass-zero constrained dynamics for simulations based on orbital-free density functional theory.

A new algorithm for efficient and fully time-reversible integration of first-principles molecular dynamics based on orbital-free density functional theory (OFDFT) is presented. The algorithm adapts to this nontrivial case, the recently introduced Mass-Zero (MaZe) constrained dynamics. The formalism ensures that full adiabatic separation is enforced between nuclear and electronic degrees of freedom and, consequently, that the exact Born-Oppenheimer probability for the nuclei is sampled. Numerical integration of the MaZe dynamics combines standard molecular dynamics algorithms, e.g., Verlet or velocity Verlet, with the SHAKE method to impose the minimum conditions on the electronic degrees of freedom as a set of constraints. The developments presented in this work, which include a bespoke adaptation of the standard SHAKE algorithm, ensure that the quasilinear scaling of OFDFT is preserved by the new method for a broad range of kinetic and exchange-correlation functionals, including nonlocal ones. The efficiency and accuracy of the approach are demonstrated via calculations of static and dynamic properties of liquid sodium in the constant energy and constant temperature ensembles.

Probing anharmonic phonons by quantum correlators: A path integral approach.

We devise an efficient scheme to determine vibrational properties from Path Integral Molecular Dynamics (PIMD) simulations. The method is based on zero-time Kubo-transformed correlation functions and captures the anharmonicity of the potential due to both temperature and quantum effects. Using analytical derivations and numerical calculations on toy-model potentials, we show that two different estimators built upon PIMD correlation functions fully characterize the phonon spectra and the anharmonicity strength. The first estimator is associated with the force-force quantum correlators and, in the weak anharmonic regime, yields reliable zero-point motion frequencies and thermodynamic properties of the quantum system. The second one is instead connected to displacement-displacement correlators and accurately probes the lowest-energy phonon excitations, regardless of the anharmonicity strength of the system. We also prove that the use of generalized eigenvalue equations, in place of the standard normal mode equations, leads to a significant speed-up in the PIMD phonon calculations, both in terms of faster convergence rate and smaller time step bias. Within this framework, using ab initio PIMD simulations, we compute phonon dispersions of diamond and of the high-pressure I41/amd phase of atomic hydrogen. We find that in the latter case, the anharmonicity is stronger than previously estimated and yields a sizeable red-shift in the vibrational spectrum of atomic hydrogen.

Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems.

Classical molecular dynamics simulations have recently become a standard tool for the study of electrochemical systems. State-of-the-art approaches represent the electrodes as perfect conductors, modeling their responses to the charge distribution of electrolytes via the so-called fluctuating charge model. These fluctuating charges are additional degrees of freedom that, in a Born-Oppenheimer spirit, adapt instantaneously to changes in the environment to keep each electrode at a constant potential. Here, we show that this model can be treated in the framework of constrained molecular dynamics, leading to a symplectic and time-reversible algorithm for the evolution of all the degrees of freedom of the system. The computational cost and the accuracy of the new method are similar to current alternative implementations of the model. The advantage lies in the accuracy and long term stability guaranteed by the formal properties of the algorithm and in the possibility to systematically introduce additional kinematic conditions of arbitrary number and form. We illustrate the performance of the constrained dynamics approach by enforcing the electroneutrality of the electrodes in a simple capacitor consisting of two graphite electrodes separated by a slab of liquid water.

Initial metal-metal bond breakage detected by fs X-ray scattering in the photolysis of Ru3(CO)12 in cyclohexane at 400 nm.

Using femtosecond resolution X-ray solution scattering at a free electron laser we were able to directly observe metal-metal bond cleavage upon photolysis at 400 nm of Ru3(CO)12, a prototype for the photochemistry of transition metal carbonyls. This leads to the known single intermediate Ru3(CO)11(µ-CO)*, with a bridging ligand (µCO) and where the asterisk indicates an open Ru3-ring. This loses a CO ligand on a picosecond time scale yielding a newly observed triple bridge intermediate, Ru3(CO)8(µ-CO)3*. This loses another CO ligand to form the previously observed Ru3(CO)10, which returns to Ru3(CO)12via the known single-bridge Ru3(CO)10(µ-CO). These results indicate that contrary to long standing hypotheses, metal-metal bond breakage is the only chemical reaction immediately following the photolysis of Ru3(CO)12 at 400 nm. Combined with previous picosecond resolution X-ray scattering data and time resolved infrared spectroscopy these results yield a new mechanism for the photolysis of Ru3(CO)12.

Computing thermal Wigner densities with the phase integration method.

We discuss how the Phase Integration Method (PIM), recently developed to compute symmetrized time correlation functions [M. Monteferrante, S. Bonella, and G. Ciccotti, Mol. Phys. 109, 3015 (2011)], can be adapted to sampling/generating the thermal Wigner density, a key ingredient, for example, in many approximate schemes for simulating quantum time dependent properties. PIM combines a path integral representation of the density with a cumulant expansion to represent the Wigner function in a form calculable via existing Monte Carlo algorithms for sampling noisy probability densities. The method is able to capture highly non-classical effects such as correlation among the momenta and coordinates parts of the density, or correlations among the momenta themselves. By using alternatives to cumulants, it can also indicate the presence of negative parts of the Wigner density. Both properties are demonstrated by comparing PIM results to those of reference quantum calculations on a set of model problems.

Modeling proton-induced damage on 2-deoxy-D-ribose. Conformational analysis.

Modeling proton-induced damage in biological systems, in particular in DNA building blocks, is of major concern in studies on cancer proton therapy. This is indeed an extremely complex process and analysis of the mechanism at the molecular level is of crucial interest. Such collision reactions of protons on biological targets induce different reactions: excitation and ionization of the biomolecule, fragmentation of the ionized species, and charge transfer from the projectile ion toward the biomolecular target. In order to have an insight into such mechanisms, we have performed a theoretical approach of two of the most important steps, the fragmentation and the charge transfer processes. For that purpose, we have considered collision of protons with isolated 2-deoxy-D-ribose by means of ab-initio molecular dynamics and quantum chemistry molecular methods. The conformation of the sugar moiety has been analyzed and appears to induce important effects, in particular different fragmentation patterns have been pointed out with regard to the conformation, and significant variations of the charge transfer cross sections have been exhibited.

Nuclear Velocity Perturbation Theory of Vibrational Circular Dichroism.

We report the first implementation of vibrational circular dichroism (VCD) within density functional theory (DFT) using the nuclear velocity perturbation (NVP) theory. In order to support VCD calculations in large-scale systems such as solvated (bio)molecules and supramolecular assemblies, we have chosen a plane-wave electronic structure code (CPMD). This implementation allows the incorporation of fully anharmonic effects in VCD spectra on the basis of ab initio molecular dynamics simulations. On the conceptual level, we compare our NVP results for rigid molecules with an existing implementation based on the magnetic field perturbation (MFP) technique using a Gaussian basis set and find an excellent agreement. Regarding numerical aspects, we analyze our results for their correct origin dependence and gauge invariance of the physical observables. The correlation with experimental data is very satisfactory, with certain deviations mainly due to the level of electronic structure theory used.

Ultrafast nonadiabatic fragmentation dynamics of doubly charged uracil in a gas phase.

A combination of time-dependent density functional theory and Born-Oppenheimer molecular dynamics methods is used to investigate fragmentation of doubly charged gas-phase uracil in collisions with 100 keV protons. The results are in good agreement with ion-ion coincidence measurements. Orbitals of similar energy and/or localized in similar bonds lead to very different fragmentation patterns, thus showing the importance of intramolecular chemical environment. In general, the observed fragments do not correspond to the energetically most favorable dissociation path, which is due to dynamical effects occurring in the first few femtoseconds after electron removal.

Theoretical investigation of the ultrafast dissociation of ionised biomolecules immersed in water: direct and indirect effects.

Theoretical simulations are particularly well suited to investigate, at a molecular level, direct and indirect effects of ionising radiations in DNA, as in the particular case of irradiation by swift heavy ions such as those used in hadron therapy. In the past recent years, we have developed the modeling at the microscopic level of the early stages of the Coulomb explosion of DNA molecules immersed in liquid water that follows the irradiation by swift heavy ions. To that end, Time-Dependent Density Functional Theory molecular dynamics simulations (TD-DFT MD) have been developed where localised Wannier orbitals are propagated. This latter enables to separate molecular orbitals of each water molecule from the molecular orbitals of the biomolecule. Our main objective is to demonstrate that the double ionisation of one molecule of the liquid sample, either one water molecule from the solvent or the biomolecule, may be in some cases responsible for the formation of an atomic oxygen as a direct consequence of the molecule Coulomb explosion. Our hypothesis is that the molecular double ionisation arising from irradiation by swift heavy ions (about 10% of ionisation events by ions whose velocity is about the third of speed of light), as a primary event, though maybe less probable than other events resulting from the electronic cascading (for instance, electronic excitations, electron attachments), may be systematically more damageable (and more lethal), as supported by experiments that have been carried out in our group in the 1990s (in studies of damages created by K holes in DNA). The chemical reactivity of the produced atomic oxygen with other radicals present in the medium will ultimately lead to chemical products that are harmful to DNA. In the present paper, we review our theoretical methodology in an attempt that the community be familiar with our new approach. Results on the production of atomic oxygen as a result of the double ionisation of water or as a result of the double ionisation of the Uracil RNA base will be presented.

Molecular dynamics study of the coordination sphere of trivalent lanthanum in a highly concentrated LiCl aqueous solution: a combined classical and ab initio approach.

The first coordination sphere of trivalent lanthanum in a highly concentrated (14 M) lithium chloride solution is studied with a combination of classical molecular dynamics and density functional theory based first principle molecular dynamics. This method enables us to obtain a solvation shell of La3+ containing 2 chloride ions and 6 water molecules. After refinement using first principle molecular dynamics, the resulting cation-water and cation-anion distances are in very good agreement with experiment. The 2Cl- and the 6 water molecules arrange in a square antiprism around La3+. Exchange of water molecules was also observed in the first-principle simulation, with an intermediate structure comprising 7 water molecules stable for 2.5 ps. Finally, evaluation of dipole moments using maximally localized Wannier functions shows a substantial polarization of the choride anions and the water molecules in the first solvation shell of trivalent lanthanum.

Ab Initio Molecular Dynamics Study of a Highly Concentrated LiCl Aqueous Solution.

The properties of a highly concentrated aqueous lithium chloride solution (|LiCl| = 14 mol L(-1)) are investigated using Car-Parrinello molecular dynamics. The coordination spheres of lithium ions, chloride ions, and water molecules are described successively. On the whole, our simulation provides results-distances and coordination numbers-in very good agreement with experimental data. The lithium solvation shell is found to exhibit a tetrahedral configuration on average, with three stable clusters observed during the simulation: Li(+)-4H2O, Li(+)(H2O)3Cl(-), and Li(+)(H2O)2(Cl(-))2. The chloride coordination sphere is logically formed by strong Cl-H hydrogen bonds with neighboring water molecules, for a mean coordination number of 4.4. The structuring of water molecules is strongly affected by the high concentration in LiCl. The hydrogen bond network is globally broken down, but little variation is calculated on water dipoles (µ = 3.07 D) because of the strong polarization from Li(+) and Cl(-) ions. We also point out some of the characteristic features of such a highly concentrated solution: water bridging between Li(+) and Cl(-) hydration spheres, Li(+)-Cl(-) ion-pairing, and intermediate behavior between dilute solutions and molten salts. Finally, the reliability of our simulation to describe ion-pairing is discussed.

Time-resolved observation of the Eigen cation in liquid water.

Experimental observation and time relaxation measurement of the hydrated proton Eigen form [H(3)O(+)(H(2)O)(3)] are presented here. Vibrational time-resolved spectroscopy is used with an original method of investigating the proton excess in water. The anharmonicity of the time-resolved spectra is characteristic of the Eigen-type proton geometry. Proton relaxation occurs in less than 200 fs. A calculation of the potential energy confirms the experimental result and the Eigen cation lifetime is in good agreement with previous molecular dynamics simulations.

Extracting effective normal modes from equilibrium dynamics at finite temperature.

A general method for obtaining effective normal modes of a molecular system from molecular dynamics simulations is presented. The method is based on a localization criterion for the Fourier transformed velocity time-correlation functions of the effective modes. For a given choice of the localization function used, the method becomes equivalent to the principal mode analysis (PMA) based on covariance matrix diagonalization. On the other hand, a proper choice of the localization function leads to a novel method with a strong analogy with the usual normal mode analysis of equilibrium structures, where the Hessian system at the minimum energy structure is replaced by the thermal averaged Hessian, although the Hessian itself is never actually calculated. This method does not introduce any extra numerical cost during the simulation and bears the same simplicity as PMA itself. It can thus be readily applied to ab initio molecular dynamics simulations. Three such examples are provided here. First we recover effective normal modes of an isolated formaldehyde molecule computed at 20 K in very good agreement with the results of a normal mode analysis performed at its equilibrium structure. We then illustrate the applicability of the method for liquid phase studies. The effective normal modes of a water molecule in liquid water and of a uracil molecule in aqueous solution can be extracted from ab initio molecular dynamics simulations of these two systems at 300 K.

Ab initio simulation of carbon clustering on an Ni(111) surface: a model of the poisoning of nickel-based catalysts.

In this work, we studied the poisoning of a nickel surface due to carbon. Performing ab initio simulations, within the framework of density functional theory, we computed the surface energy of the nickel (111) surface as a function of carbon coverage. On the basis of these results, we can assert that the stable state of the nickel/carbon surface is either a clean nickel surface or a fully carbon-covered nickel surface, which has a graphitic configuration. The relative stability of the two states depends on the temperature and partial pressure of the carbon gas. At fixed nominal carbon coverage, the most stable configurations are those forming carbon clusters. However, the nickel sites hosting these clusters change from hexagonal close packed/face centered cubic (hcp/fcc) sites to on-top sites when going toward larger clusters. This indicates that poisoning due to graphitic patches occurs on on-top sites.

Recombination of photodissociated iodine: a time-resolved x-ray-diffraction study.

A time-resolved x-ray-diffraction experiment is presented that aims to study the recombination of laser-dissociated iodine molecules dissolved in CCl4. This process is monitored over an extended time interval from pico- to microseconds. The variations of atom-atom distances are probed with a milliangstrom resolution. A recent theory of time-resolved x-ray diffraction is used to analyze the experimental data; it employs the correlation function approach of statistical mechanics. The most striking outcome of this study is the experimental determination of time-dependent I-I atom-atom distribution functions. The structure of the CCl4 solvent changes simultaneously; the solvent thus appears as a reaction partner rather than an inert medium hosting it. Thermal expansion of the system is nonuniform in time, an effect due to the presence of the acoustic horizon. One concludes that a time-resolved x-ray diffraction permits real-time visualization of solvent and solute motions during a chemical reaction.

Infrared Spectroscopy of N-Methylacetamide Revisited by ab Initio Molecular Dynamics Simulations.

The density functional theory based molecular dynamics simulation method ("Car-Parrinello") was applied in a numerical study of the electronic properties, hydrogen bonding, and infrared spectroscopy of the trans and cis isomer of N-methylacetamide in aqueous solution. A detailed analysis of the electronic structure of the solvated molecules, in terms of localized Wannier functions and Born atomic charges, is presented. Two schemes for the computation of the solute infrared absorption spectrum are investigated: In the first method the spectrum is determined by Fourier transforming the time correlation function of the solute dipole as determined from the Wannier function analysis. The second method uses instead the molecular current-current correlation function computed from the Born charges and atomic velocities. The resulting spectral properties of trans- and cis-NMA are carefully compared to each other and to experimental results. We find that the two solvated isomers can be clearly distinguished by their infrared spectral profile in the 1000-2000 cm(-)(1) range.