Project topics and materials
Get Complete Project Material File(s) Now! »
Thermodynamics and phase transitions
Generally speaking, physical systems have several possible states that they can reach for a given pressure and temperature. For instance, at ambient pressure and summer temperature, water can be liquid, solid or gaseous, but if we take a piece of ice, it will melt into liquid. We say that the liquid state is stable, while the ice state is metastable, and that our piece of ice has performed a phase transition.
To describes this, historically physicists introduced the concept of free energy.
It is a function of thermodynamic quantities, like pressure P and temperature T , and of some order parameter. Without entering too much into the details, an order parameter is a function that describe in which state our system is, telling us for instance if it is liquid or solid. At fixed temperature and pressure, we speak about the Gibbs free energy and will note it G. Metastable states correspond to local minima of G, while the stable state is the absolute minimum of G ; different locally stable states are separated by free energy barriers.
In general phase transitions can be put into two group, depending on the conti-nuous properties of G. We speak about first order phase transition when first deriva-tives of G are discontinuous, like latent heat or density. We speak about second order transition when second derivatives of G are discontinuous, like heat capacity. Despite this relatively abstract definition, one important consequence of the discontinuous nature of first order transitions is that during transition, the system will have some of its parts that have completed the transition and others that haven’t, forming an interface between the two that the system will tend to minimize due to its energetic cost. Like for instance when water boils and bubbles of vapor form in it.
The ensemble of equilibrium stability information is generally grouped in a phase diagram, where we show the preferred physical states of matter at different thermo-dynamic variables. In this diagram we have open spaces where a single state is stable, separated by lines where phase transition will occur, that we call phase boundaries. In other terms, the open space corresponds to free energy with one single global mi-nimum,whereas along the lines we have two or three global minima with the same stability. The most common phase diagram uses pressure and temperature as ther-modynamic variables.
If the free energy information contained in a phase diagram is sufficient to des-cribe stability of the states available by a system and the nature of the transition bet-ween two states, it has some limitations. As it is an equilibrium quantity, it does not inform us about the out of equilibrium dynamics, like transition pathways and their kinetics. When studying kinetics, one generally wants to reduce the complex high-dimensional phase space dynamics of the system described by the Liouville equation to a more handy and human readable description. This can be done by coarse grai-ning into a few relevant order parameters (or into a discrete set of microstates) and using stochastic equations that represent the average behavior of exact trajectories, to represent how the system passes from one state to another. In this scheme, each minimum of G will effectively trap the system for a time long (from microseconds to millions years) in comparison to typical bond vibrations (few femtoseconds). Rarely, a state will rapidly jump over the barrier into another metastable state.
Studying the kinetics is a more subtle and complex problem than the study of free energy, but a complete reconstruction of kinetics contains more information and is essential to understand the experimental world, where very often the physical systems do not have the time to reach the equilibrium state. A feature that is exploited by organisms to maintain their living state, and also by humans for the synthesis of complex materials for technological applications. In this thesis, we will address both the thermodynamics and kinetics of transformation processes in water.
Numerical study of transformations of matter
As briefly stated previously, transformation processes are hard to grasp expe-rimentally, due to the small time and length scale involved. In theory, molecular dynamics simulations would be ideal to study such phenomena at microscopic scale, since long trajectories of the system would sample fluctuations within metastable states and all their transitions among them. In practice however, there are two main limitations.
The first one is linked to the way we describe our material numerically. A quan-tum approach is a priori more precise as it relies on first principles quantum mecha-nics, but it drastically limits the length scale accessible during a simulation due to computational cost. Thus it is in general of limited interest when we study phase transitions, as if size effect are important for the transition (and they generally are), the limited size of the system under study would spoil the meaning of the extracted data. Thus classical approaches are generally used when studying transformations of matter, as they give access to length scales thousands of times larger than in quantum simulation. The problem is that now we need to represent our material with some predefined force fields, generally tuned to reproduce at best properties found with either quantum simulation or experimental data. As you may guess, this approach can only give approximate results, even if some models reproduce quite well selected properties of the material under study.
The second limitation come from the current time scale available by compu-tation, as it is orders of magnitude smaller than what we would need for proper sampling of transitions in condensed matter. Quantum simulations are limited to the sub-nanosecond time scale, and classical atomistic simulations are limited to the sub-milliseconds time scale. (For instance ice nucleation might require from millise-cond to millions years depending on conditions). This limits drastically the range of transformations that one can study. To solve this problem and also to get insight not automatically provided by a simple long trajectory, several methods were developed in the last decades, called enhanced sampling methods [46].
Enhanced sampling methods
Despite being grouped under one terminology, enhanced sampling methods tend to solve two different problems : accelerated exploration of the available metastable states and/or precise sampling of the thermodynamic and kinetic properties. The latter consists in large accumulation of samples in relevant regions of the configuration space, to construct an estimate of equilibrium probability distributions (and possibly kinetic properties). In the present thesis, we used enhanced sampling methods that exploited the two following ideas :
1. add artificial, external biasing force to the natural forces of the system, in a way that enhances population of barrier regions, compared to its negligible equilibrium value.
2. generate several trajectories starting from specific configurations, only kee-ping those that fulfill some clever requirements.
Even if the methods are general and do not restrain themselves to the study of a specific material, what limits their use is the need to define some variables that will allow us to track the transformations studied. Such variables are generally called structural descriptors, collective variables, order parameters, or reaction coordinates by different scientific communities. We will talk about reaction coordinates or order parameters when their aim is to describe a specific reaction mechanism.
Reaction coordinates
Defining a reaction coordinate is not a difficult task, but it can be extremely chal-lenging to define a good reaction coordinate. Often this relies on specific knowledge about the system under study and tedious trial-and-error iterative process. Postpo-ning the discussion about methods to quantitatively assess the quality of a variable, it is important to note that efficiency of enhanced sampling methods is directly lin-ked to the quality of the reaction coordinates. Moreover, finding optimal reaction coordinates can lead to deeper understanding of the transformation processes.
A recently developed approach to define good reaction coordinates, regardless of the system under study, starts by describing transformations of matter as changes in a matrix that contains the inter-atomic bond network information [47, 7, 8]. One aim of this thesis is to apply, improve and further validate this approach by applying it to the study of water in the low temperature regime.
Complexity of water
Despite its molecular simplicity, water reveals to be very complex from the view-point of physico-chemical properties. These include numerous triple points, at least one critical point, more than 10 stable solid phases, no less than 3 metastable amor-phous phases and several theoretical liquid states [48], without speaking of its 74 anomalies compared to other liquids ! Among them, one of the most famous is that water expands upon freezing, increasing its volume by 9% under atmospheric pressure, implying that ice floats on water. Another one is its high density, with a maximum at 4◦, meaning that upon cooling or heating liquid water its density will decrease. This is why for instance the bottom of fresh water lakes keeps the same temperature regardless of the external temperature variation, as water at 4◦ will sink due to its higher density. This temperature-driven shift in density is also the origin of ocean currents such as the gulf’s stream, which have huge impact on land climate.
Just to have a glimpse of its complexity, figure 2.4.1 presents the stable phase diagram of water [48]. Every solid line represents a phase boundary. Along them, two phases will stably coexist in any relative proportions. When three such lines join, we have a triple point with three stably coexisting phases. As already mentioned, it is important to note that such diagram indicates the equilibrium properties of water, and does not tell us about how the phases are kinetically connected, or by which microscopic mechanisms water transforms from one phase to another when crossing a line.
The complexity of collective behavior of water is linked to its polarized molecule H2O. This polarization is due to the difference in electronegativity between oxygen and hydrogen nucleus, which leads to the hydrogen’s electron being more attracted by the oxygen atom, resulting in a slightly negatively charged oxygen atom, while hydrogen is slightly positively charged. As a result water molecules will form hydrogen bonds, where one of the hydrogen atom is linked to the oxygen atom of another molecule. Compared to other types of chemical bonds, the hydrogen bond is neither strong nor weak, meaning that it can be easily broken but will generally survive to thermal fluctuations [49].
At moderate pressures, each water molecule can form 4 hydrogen bonds, 2 in-volving its oxygen atom and 2 involving its hydrogen atoms. This 5 water molecules bonded together will optimally arrange themselves in tetrahedral shapes as shown in figure 2.4.2. In solid phases this local tetrahedral arrangement will extend to the whole system and produce crystalline structures. In liquid phases thermal energy will break, stretch or bend the hydrogen bonds, leading to only local clusters of tetrahe-dral structure, even if large chains of hydrogen-bonded molecules are present. At high pressure this tetrahedral structure will be enriched by supplementary water molecules and the soft hydrogen bond will evolved toward a strong shared proton due to nuclear quantum effect [50]. Properly describing water with classical force fields is not an easy task. As a result, until quite recently they was no model able to reproduce correctly water phy-sical properties, with the correct phase diagram and the correct structural properties of each phase [5].
Here, we will restrict ourselves to the study of water at low temperature, from 140 to 237 K, and moderate pressure, from 1 to 5000 atmospheres (0.1 MPa to 0.5 GPa). In this region, two transformations are of main interest : the supposed liquid-liquid transition of water and the homogeneous ice nucleation of water.
The “liquid-liquid transition”
If we look at water at very low temperature (less than 140 K), we see that water possesses at least three amorphous metastable phases : a low density amor-phous (LDA), a high density amorphous (HDA) and a very high density amorphous (VHDA). This property of a solid to exist under different amorphous form is called polyamorphism. Evidence of these three forms cames from experiments [17, 16, 18]. The low and high density amorphous phases are connected by a reversible first-order transition [19].
To explain the existence of these states and some of the water anomalies, several theories were formulated in the last three decades. Among them, one postulates the existence of two structurally different liquid, with a first order transition between them [20]. This hypothesis is hard to be verified experimentally in pure water. Indeed, its supposed location is right in an area, the so called no man’s land, where liquid water spontaneously transforms into ice, preventing any experimental measure on the underlying metastable liquid states. Its numerical study is also not trivial to pursue, as one needs to distinguish two similar liquid structures and to perform simulation at low temperature with slow dynamics, that renders very costly the statistically meaningful sampling of the corresponding states.
In fact, an uninterrupted series of experimental works have tried to probe the existence or absence of this transition, using salty water or confined systems to prevent freezing [51, 52, 53, 54, 55]. Results show that in salty water, despite existence of a LDA-HDA first order transition, a liquid-liquid transition does not exist [22, 21]. But for pure water the nature of the no man’s land remains unobserved.
In parallel, the computational scientific community has not been at rest, with an ongoing debate about the very existence of the liquid-liquid transition. Most of the debate crystallized around the sampling methods and the model used to describe water during the simulation. After long controversies, a study based on a sub-optimal model of water found clear and indubitable evidence of a first order liquid-liquid transition [23]. Whilst other studies found indirect evidence that this transition might not exist using a more recent and precise model of water [56, 57].
Most likely the debate in the scientific community will continue until experiments bring it to an end. Still, one of the aims of this thesis is to clarify the situation with robust evidence of the existence or absence of the transition with the use of a state-of-the-art water model and of enhanced sampling methods.
Homogeneous ice nucleation
Everyone knows that water freezes when it comes down to 273,15 K (0◦ C), but few know that in pure water the kinetic barrier is so large that freezing would require more than the age of the universe [36]. In practical situations, the freezing is accelerated by the presence of impurities or interfaces, that favor the apparition of small nuclei of ice subsequently growing until filling the whole system once they have attained some critical size. In the case of perfectly pure water we speak about homogeneous ice nucleation, otherwise we speak about heterogeneous ice nucleation. To study the mechanism of homogeneous ice nucleation, experiments are limited by the difficulty to prepare pure water samples and by the short time and length scale involved in the microscopic nucleation processes. One may ask what are the structural properties of the initial nucleus of ice, but for instance at 230 K, its size is only of a few nanometers and its growth can take place in few hundreds of nanoseconds. Note that if the growth of the nucleus is fast, its apparition is a rare event, so one may need to wait milliseconds or seconds to see one appear and propagate to the whole system.
We want to stress that despite the experimental difficulties in realizing homoge-neous nucleation of ice, the process remains of fundamental interest to understand crystallization and to address the more complex case of heterogeneous ice nucleation. The latter is more readily observed in real life, however it has additional difficulties related to the precise experimental characterization of the nucleation sites (defects, impurities, etc.). Also from a numerical viewpoint, further complications arise from the definition of interactions between different molecular or atomic species. It is a very challenging domain, where theoretical predictions are still scarce and often in strong disagreement with experiments. This is also linked to the limitation of the simplified models used, like classical nucleation theory.
If the time and length scales of detailed nucleation mechanisms are too small for experiments, they are on the reverse too big for simple numerical study. In fact, increasing the number of molecules increases rapidly the cost of the simulation, so either we can use less precise simulations or we can resort to smaller number of molecules. But even by limiting the number of molecules, milliseconds or seconds are just too much for nowadays supercomputer. So we need to use enhanced sampling methods, with all the complications around the definition of a reaction coordinate implied.
One of the aims of this thesis is to study nucleation with state-of-the-art molecu-lar models and enhanced sampling methods, coupled with the recent general approach to compute reaction coordinates based on the adjacency matrix.
Open challenges
As discussed above, on the one hand we need to develop general tools and me-thods to study transformations of matter. On the other hand water is a material that despite its molecular simplicity gives rise to highly complex behaviors, and requires special care to be investigated, especially in the supercooled region.
Can we successfully develop general methods to study transformations of matter ? And can we apply them to the difficult study of both the liquid-liquid transition and homogeneous ice nucleation of water ?
For the liquid-liquid transition, we need to study water at several low temperature conditions, to see how the stability of the hypothesized high and low density liquid states evolves and in which way they are mixing or coexisting together. In particular, a recurring question in the community concerns the reconstruction of free energy landscapes using accurate inter-atomic potentials, which is one of the aims of our work.
For homogeneous ice nucleation, we need to study the critical nuclei and whe-ther and how they spontaneously evolve from a “pure” crystalline state with unique symmetry to a disordered stack of hexagonal and cubic ice layers. We remark that the mechanism is believed to be far from trivial, with a complex interplay of different crystalline phases, and to defy classical nucleation theory.
To achieve these targets, we realised massive classical molecular dynamics si-mulations (more than 10 millions of cpu hours) coupled with a range of different enhanced sampling methods, to systematically reconstruct the mechanisms, the free energy landscapes and the kinetics of several transformation processes in water.
We anticipate that for the first time we reconstructed accurate free energy land-scapes in no man’s land to study the liquid-liquid transition with the TIP4P/2005 potential, observing the lack of free energy barriers and of a discontinuous transition, similarly to mW and contrary to the ST2 potential. We also provide robust evidence of the nucleation mechanism, starting from hexagonal or cubic nuclei that sponta-neously evolve toward stacking disordered nuclei, and eventually ice, applying for the first time rigorous aimless shooting techniques to the accurate TIP4P/ice potential. Our results provide also a solid background and database for further investigations of nucleation free energy landscapes and kinetic rates.
This thesis is structured in two parts. The first one presents the various methods and numerical tools we used, with their underlying physical and mathematical prin-ciples. The second one presents results obtained during the study of the liquid-liquid transition, or of the homogeneous nucleation of ice.
Table of contents :
1 Résumé en français
2 Introduction
2.1 Context
2.2 Thermodynamics and phase transitions
2.3 Numerical study of transformations of matter
2.3.1 Enhanced sampling methods
2.3.2 Reaction coordinates
2.4 Complexity of water
2.4.1 The “liquid-liquid transition”
2.4.2 Homogeneous ice nucleation
2.5 Open challenges
I Methods
3 Simulations of Condensed Matter
3.1 Classical molecular dynamics
3.1.1 Short introduction to statistical principle
3.1.2 Time evolution
3.1.3 Temperature coupling
3.1.4 Pressure coupling
3.1.5 Periodic boundary conditions
3.2 Force fields
3.2.1 The mW model
3.2.2 The ST2 model
3.2.3 The TIP4P family
3.3 About the choice of water model
4 Describing our systems with collective variables
4.1 Local collective variables
4.1.1 Coordination number
4.1.2 The Steinhardt parameters
4.1.3 The Chill+ algorithm
4.2 Global collective variable
4.2.1 The largest nucleus size
4.2.2 Cubicity
4.2.3 Permutation Invariant Vector (PIV) and PIV distance
4.2.4 Path collective variables
4.2.5 Commitment probability
4.3 Quality of reaction coordinates
4.3.1 Maximum likelihood optimization
4.4 Development of PIV in plumed
4.4.1 User interface
4.4.2 Counting sort algorithm
4.4.3 Computing the derivatives
5 Enhanced sampling methods
5.1 Free energy exploration and reconstruction
5.1.1 Umbrella sampling
5.1.2 Metadynamics
5.2 Transition path sampling
5.2.1 Seeding
5.2.2 Aimless shooting algorithm
5.2.3 Shooting range algorithm
II Results
6 Study of the Liquid-Liquid Phase Transition
6.1 Introduction
6.1.1 The two-states model
6.1.2 Experimental studies of the liquid-liquid transition
6.1.3 Numerical studies of the liquid-liquid transition
6.2 Simulation methods
6.2.1 Molecular dynamics parameters
6.2.2 States preparation
6.2.3 Order parameter definition
6.2.4 Free energy calculations
6.3 Exploration of the (P, T) diagram
6.3.1 First exploration using Berendsen barostat
6.3.2 Proper exploration with ensemble consistent barostat
6.3.3 Convergence assessment and error estimation
6.3.4 Schematic no man’s land phase diagram
6.4 Kinetic properties of the explored (P, T) conditions
6.5 Structural properties
6.5.1 Coordination number
6.5.2 Cluster analysis
6.6 Discussion
6.7 Conclusions and outlook
7 Study of the Homogeneous Ice Nucleation
7.1 Introduction
7.1.1 Classical nucleation theory
7.1.2 The ice I polymorphs
7.1.3 A bit of crystallography
7.1.4 Numerical study of nucleation of ice I
7.2 Generation of initial reactive trajectories
7.2.1 Freezing and melting with metadynamics
7.2.2 Exploration of the transition state with seeding
7.3 Order parameter quality
7.3.1 Choice of the order parameter
7.3.2 committor analysis
7.3.3 Maximum likelihood optimization
7.4 Sampling of the transition path ensemble
7.4.1 Standard aimless shooting
7.4.2 Aimless shooting within a range
7.4.3 Difference of efficiency between the two techniques
7.4.4 Critical nuclei evolution
7.5 Conclusions and outlook
8 Conclusion
References
