Determining the phase diagram of water from direct coexistence simulations: The phase diagram of the TIP4P/2005 model revisited

Abstract

Computing phase diagrams of model systems is an essential part of computational condensed matter physics. In this paper we discuss in detail the interface pinning (IP) method for calculation of the Gibbs free energy difference between a solid and a liquid. This is done in a single equilibrium simulation by applying a harmonic field that biases the system towards two-phase configurations. The Gibbs free energy difference between the phases is determined from the average force that the applied field exerts on the system. As a test system we study the Lennard-Jones model. It is shown that the coexistence line can be computed efficiently to a high precision when the IP method is combined with the Newton-Raphson method for finding roots. Statistical and systematic errors are investigated. Advantages and drawbacks of the IP method are discussed. The high pressure part of the temperature-density coexistence region is outlined by isomorphs. I. INTRODUCTIONAn important aspect of computational condensed matter physics is to compute phase diagrams of model systems. The naive approach is to preform a long-time simulation at a selected state point and hope that the system by itself finds its preferred phase, i.e. the phase with the lowest Gibbs free energy. In most cases this strategy is not viable since first-order transitions are associated with hysteresis. Thus the system is likely to be stuck in a metastable phase. The origin of this hysteresis effect is the formation of an interface between two phases. The surface tension of the interface gives rise to a free energy barrier that the system has to overcome to transform from one phase to the other [1]. A conceptual appealing and widely used strategy to overcome the hysteresis problem is to preform simulations starting from an initial configuration with two phases in a periodic box [2][3][4][5][6][7][8][9][10][11][12][13][14]. When a steady state situation is reached the stable phase will grow at the expense of the other phase. The disadvantages of this approach are that: i) it relates to a non-equilibrium computation that cannot be done ad infinitum; ii) a sufficiently large system is needed to reach the steady state growth; iii) thermal fluctuation will result in some probability that the system will move towards the disfavored phase. In a recent paper [15] these problems were resolved by introducing the so-called "interface pinning" (IP) method. In short, the idea of this method is to compute the average force needed to keep the system in the two-phases state. This is done by connecting the system to a harmonic field which couples to an order-parameter that distinguished between the two phases of interest. The Gibbs free energy difference between the phases is determined by the average force that the applied field exerts on the system. Thus, a standard * ulf.pedersen@tuwien.ac.at ad infinitum equilibrium simulation gives the information needed to computed the Gibbs free energy difference between the two phases of interest.The purpose of this paper is to give a detailed ...

Searching for crystal-ice domains in amorphous ices

“…1 calculated using, as a reference, the unit cells of ice (a) α =II, (b) α =III and (c) α =IV for LDA (black) and for HDA Ih (red) at p= 0.01 GPa and at p= 3.0 GPa, respectively. The green distrubutions represent the histograms for the corresponding bulk ices computed at the following thermodynamic conditions [44]: T = 210 K and p = 4 GPa for ice II, and T = 250 K and p = 3 GPa for ice III and IV. The green line emphasizes the value S = 0.8 used as a cutoff to identify ice-like environments.…”

Competition between ices Ih and Ic in homogeneous water freezing

“…59 This potential provides a reasonably good global description of real water. 59,68 Although TIP4P/2005 is probably the rigid non-polarizable model that best describes properties of liquid water, 63 it fails in predicting the melting temperature, underestimating it by more than 20 K. 59 By contrast, the TIP4P/Ice model 60 has been tailored to closely reproduce the coexistence temperature between ice Ih and liquid water at normal pressure. In this work we use both TIP4P/2005 and TIP4P/Ice to study the competition between ice Ih and ice Ic in the nucleation of ice from supercooled water at normal pressure.…”

High precision determination of the melting points of water TIP4P/2005 and water TIP4P/Ice models by the direct coexistence technique

“…In recent years, the technique of direct coexistence has been implemented successfully to calculate the equilibrium conditions, and thus melting temperatures, in different real systems such as solid water phases, [51][52][53][54] ionic compounds, 55,56 clathrate hydrates, [57][58][59] or methanol. 60 This method was originally proposed by Ladd and Woodcock 61,62 for a Lennard-Jones system.…”