Stabilization of Ab Initio Molecular Dynamics Simulations at Large Time Steps

The Verlet method is still widely used to integrate the equations of motion in ab initio molecular dynamics simulations. We show that the stability limit of the Verlet method may be significantly increased by setting an upper limit on the kinetic energy of each atom with only a small loss in accuracy. The validity of this approach is demonstrated for molten lithium fluoride.Comment: 9 pages, 3 figure

Kinetic energy of protons in ice Ih and water: a path integral study

The kinetic energy of H and O nuclei has been studied by path integral molecular dynamics simulations of ice Ih and water at ambient pressure. The simulations were performed by using the q-TIP4P/F model, a point charge empirical potential that includes molecular flexibility and anharmonicity in the OH stretch of the water molecule. Ice Ih was studied in a temperature range between 210-290 K, and water between 230-320 K. Simulations of an isolated water molecule were performed in the range 210-320 K to estimate the contribution of the intramolecular vibrational modes to the kinetic energy. Our results for the proton kinetic energy, K_H, in water and ice Ih show both agreement and discrepancies with different published data based on deep inelastic neutron scattering experiments. Agreement is found for water at the experimental melting point and in the range 290-300 K. Discrepancies arise because data derived from the scattering experiments predict in water two maxima of K_H around 270 K and 277 K, and that K_H is lower in ice than in water at 269 K. As a check of the validity of the employed water potential, we show that our simulations are consistent with other experimental thermodynamic properties related to K_H, as the temperature dependence of the liquid density, the heat capacity of water and ice at constant pressure, and the isotopic shift in the melting temperature of ice upon isotopic substitution of either H or O atoms. Moreover, the temperature dependence of K_H predicted by the q-TIP4P/F model for ice Ih is found to be in good agreement to results of path integral simulations using ab initio density functional theory.Comment: 11 pages, 6 figures, 2 table

Path-integral molecular dynamics simulation of 3C-SiC

Molecular dynamics simulations of 3C-SiC have been performed as a function of pressure and temperature. These simulations treat both electrons and atomic nuclei by quantum mechanical methods. While the electronic structure of the solid is described by an efficient tight-binding Hamiltonian, the nuclei dynamics is treated by the path integral formulation of statistical mechanics. To assess the relevance of nuclear quantum effects, the results of quantum simulations are compared to others where either the Si nuclei, the C nuclei or both atomic nuclei are treated as classical particles. We find that the experimental thermal expansion of 3C-SiC is realistically reproduced by our simulations. The calculated bulk modulus of 3C-SiC and its pressure derivative at room temperature show also good agreement with the available experimental data. The effect of the electron-phonon interaction on the direct electronic gap of 3C-SiC has been calculated as a function of temperature and related to results obtained for bulk diamond and Si. Comparison to available experimental data shows satisfactory agreement, although we observe that the employed tight-binding model tends to overestimate the magnitude of the electron-phonon interaction. The effect of treating the atomic nuclei as classical particles on the direct gap of 3C-SiC has been assessed. We find that non-linear quantum effects related to the atomic masses are particularly relevant at temperatures below 250 K.Comment: 14 pages, 15 figure

Holonomic constraints : an analytical result

Systems subjected to holonomic constraints follow quite complicated dynamics that could not be described easily with Hamiltonian or Lagrangian dynamics. The influence of holonomic constraints in equations of motions is taken into account by using Lagrange multipliers. Finding the value of the Lagrange multipliers allows to compute the forces induced by the constraints and therefore, to integrate the equations of motions of the system. Computing analytically the Lagrange multipliers for a constrained system may be a difficult task that is depending on the complexity of systems. For complex systems, it is most of the time impossible to achieve. In computer simulations, some algorithms using iterative procedures estimate numerically Lagrange multipliers or constraint forces by correcting the unconstrained trajectory. In this work, we provide an analytical computation of the Lagrange multipliers for a set of linear holonomic constraints with an arbitrary number of bonds of constant length. In the appendix of the paper, one would find explicit formulas for Lagrange multipliers for systems having 1, 2, 3, 4 and 5 bonds of constant length, linearly connected.Comment: 13 pages, no figures. To appear in J. Phys. A : Math. The

Classical-path integral adaptive resolution in molecular simulation: towards a smooth quantum-classical coupling

Simulations that couple different classical molecular models in an adaptive way by changing the number of degrees of freedom on the fly, are available within reasonably consistent theoretical frameworks. The same does not occur when it comes to classical-quantum adaptivity. The main reason for this is the difficulty in describing a continuous transition between the two different kind of physical principles: probabilistic for the quantum and deterministic for the classical. Here we report the basic principles of an algorithm that allows for a continuous and smooth transition by employing the path integral description of atoms.Comment: 8 pages 4 figure

Path Integral Molecular Dynamics within the Grand Canonical-like Adaptive Resolution Technique: Simulation of Liquid Water

Quantum effects due to the spatial delocalization of light atoms are treated in molecular simulation via the path integral technique. Among several methods, Path Integral (PI) Molecular Dynamics (MD) is nowadays a powerful tool to investigate properties induced by spatial delocalization of atoms; however computationally this technique is very demanding. The abovementioned limitation implies the restriction of PIMD applications to relatively small systems and short time scales. One possible solution to overcome size and time limitation is to introduce PIMD algorithms into the Adaptive Resolution Simulation Scheme (AdResS). AdResS requires a relatively small region treated at path integral level and embeds it into a large molecular reservoir consisting of generic spherical coarse grained molecules. It was previously shown that the realization of the idea above, at a simple level, produced reasonable results for toy systems or simple/test systems like liquid parahydrogen. Encouraged by previous results, in this paper we show the simulation of liquid water at room conditions where AdResS, in its latest and more accurate Grand-Canonical-like version (GC-AdResS), is merged with two of the most relevant PIMD techniques available in literature. The comparison of our results with those reported in literature and/or with those obtained from full PIMD simulations shows a highly satisfactory agreement

Hydrogen and muonium in diamond: A path-integral molecular dynamics simulation

Isolated hydrogen, deuterium, and muonium in diamond have been studied by path-integral molecular dynamics simulations in the canonical ensemble. Finite-temperature properties of these point defects were analyzed in the range from 100 to 800 K. Interatomic interactions were modeled by a tight-binding potential fitted to density-functional calculations. The most stable position for these hydrogenic impurities is found at the C-C bond center. Vibrational frequencies have been obtained from a linear-response approach, based on correlations of atom displacements at finite temperatures. The results show a large anharmonic effect in impurity vibrations at the bond center site, which hardens the vibrational modes with respect to a harmonic approximation. Zero-point motion causes an appreciable shift of the defect level in the electronic gap, as a consequence of electron-phonon interaction. This defect level goes down by 70 meV when replacing hydrogen by muonium.Comment: 11 pages, 8 figure

On the Geometry and Entropy of Non-Hamiltonian Phase Space

We analyze the equilibrium statistical mechanics of canonical, non-canonical and non-Hamiltonian equations of motion by throwing light into the peculiar geometric structure of phase space. Some fundamental issues regarding time translation and phase space measure are clarified. In particular, we emphasize that a phase space measure should be defined by means of the Jacobian of the transformation between different types of coordinates since such a determinant is different from zero in the non-canonical case even if the phase space compressibility is null. Instead, the Jacobian determinant associated with phase space flows is unity whenever non-canonical coordinates lead to a vanishing compressibility, so that its use in order to define a measure may not be always correct. To better illustrate this point, we derive a mathematical condition for defining non-Hamiltonian phase space flows with zero compressibility. The Jacobian determinant associated with time evolution in phase space is altogether useful for analyzing time translation invariance. The proper definition of a phase space measure is particularly important when defining the entropy functional in the canonical, non-canonical, and non-Hamiltonian cases. We show how the use of relative entropies can circumvent some subtle problems that are encountered when dealing with continuous probability distributions and phase space measures. Finally, a maximum (relative) entropy principle is formulated for non-canonical and non-Hamiltonian phase space flows.Comment: revised introductio

Alchemical normal modes unify chemical space

In silico design of new molecules and materials with desirable quantum properties by high-throughput screening is a major challenge due to the high dimensionality of chemical space. To facilitate its navigation, we present a unification of coordinate and composition space in terms of alchemical normal modes (ANMs) which result from second order perturbation theory. ANMs assume a predominantly smooth nature of chemical space and form a basis in which new compounds can be expanded and identified. We showcase the use of ANMs for the energetics of the iso-electronic series of diatomics with 14 electrons, BN doped benzene derivatives (C$_{6-2x}$(BN)$_{x}$H$_6$ with $x = 0, 1, 2, 3$), predictions for over 1.8 million BN doped coronene derivatives, and genetic energy optimizations in the entire BN doped coronene space. Using Ge lattice scans as reference, the applicability ANMs across the periodic table is demonstrated for III-V and IV-IV-semiconductors Si, Sn, SiGe, SnGe, SiSn, as well as AlP, AlAs, AlSb, GaP, GaAs, GaSb, InP, InAs, and InSb. Analysis of our results indicates simple qualitative structure property rules for estimating energetic rankings among isomers. Useful quantitative estimates can also be obtained when few atoms are changed to neighboring or lower lying elements in the periodic table. The quality of the predictions often increases with the symmetry of system chosen as reference due to cancellation of odd order terms. Rooted in perturbation theory the ANM approach promises to generally enable unbiased compound exploration campaigns at reduced computational cost