- Email: [email protected]
Location of phase equilibria by temperature-quench molecular dynamics simulations
Fluid Phase Equilibria 203 (2002) 1–14
Location of phase equilibria by temperature-quench molecular dynamics simulations Lev D. Gelb a,1 , Erich A. Müller b,∗ b
a Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL 32306-4390, USA Departamento de Termodinámica y Fenómenos de Transferencia, Universidad Simón Bol´ıvar, Caracas 1080, Venezuela
Received 2 January 2002; accepted 14 May 2002
Abstract We present a method to locate phase coexistence points using molecular dynamics (MD) simulations. By equilibrating a single phase in a relatively large simulation cell under NVT conditions and then quenching the system into a two-phase region, it spontaneously separates into two coexisting phases. The different phases are identified by collecting histograms of the local density and/or local composition of the pure phase regions, which can be done soon after the quench. The method can be used to locate vapor–liquid, liquid–liquid or solid–fluid equilibria. MD is a robust simulation technique that can be applied without difficulty to very complex molecules and mixtures, and can be implemented on large parallel computers. This approach, thus, has the potential to be both faster (in wall clock time) and more generally applicable than any other single technique. The method is demonstrated on test systems of single component and binary Lennard–Jones (LJ) fluids. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Molecular simulation; Method of calculation; Vapor–liquid equilibrium; Solid–fluid equilibria
1. Introduction The most popular method for locating phase coexistence in molecular simulation is still Gibbs Ensemble Monte Carlo [1–3], though it is being slowly supplanted by histogram-reweighting (HR) techniques, which offer substantially more precise location of critical points [4]. Monte Carlo methods have three well-known deficiencies. When simulating dense phases, typical of strongly interacting liquids and solids, equilibration is difficult to achieve because of the poor statistics associated with the insertion/deletion steps. Secondly, they are difficult to apply to systems containing very complex molecules without substantial system-specific modifications (such as Configurational-Bias Monte Carlo [5–7]). ∗
Corresponding author. E-mail addresses: [email protected] (L.D. Gelb), [email protected] (E.A. Müller). 1 Present address: Department of Chemistry, Washington University , St. Louis, MO 63130-1134, USA. 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 2 ) 0 0 1 7 4 - 7
2
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Thirdly, Monte Carlo methods are difficult to implement on parallel computers and few efforts have been made to improve this situation [8,9]. Kofke [10,11] suggested “Gibbs–Duhem integration”, an indirect method based on the integration of the Clausius–Clapeyron equation. However, to apply Gibbs–Duhem integration, one must start with a known phase coexistence point. The “NPT + test particle” method [12] may also be used to study phase equilibria, but suffers from the same limitations as other Monte Carlo simulations, since it relies on a particle insertion technique. In contrast, the molecular dynamics (MD) simulation method is applicable to both dilute and dense phases, easily parallelized, and widely used for simulating very large and complex molecules. Several parallelized MD codes are available, including AMBER [13], CHARMM [14], DL POLY [15], Gromos96 [16], Moldy [17], and Gromacs [18]. The application of MD to study phase equilibria is not new, and in fact predates GEMC. A typical implementation is that of Holcomb et al. [19,20], in which previously equilibrated fractions of a bulk liquid and bulk vapor phase are placed in a simulation box in the form of a slab of liquid surrounded by a vapor phase. The system is then allowed to evolve under NVT conditions until equilibrium is reached through diffusive mass transport. Goujon et al. [21] have applied a similar method, using MC instead of MD. This approach is still used to study interfacial properties [22–24], and while it can provide VLE or LLE coexistence data, it is not considered efficient because long times are required for complete equilibration of a two-phase system. In such systems, the last stage of equilibration consists of the minimization of the (already small) surface area between the two phases, which provides only a relatively weak driving force. What follows is the description of a method that allows the location of phase coexistence through a constant density simulation in which the temperature is changed in a single time-step (quenching) in order to place the system in a thermodynamically and mechanically unstable state. At this point, the system spontaneously separates into domains of stable, locally equilibrated phases. These phases are then detected and characterized by dividing the system into small blocks and analyzing the densities and/or the compositions in the blocks. Phase coexistence data is, therefore, gathered from the locally equilibrated system, and it is not necessary to continue the simulation until global equilibrium is reached, resulting in a dramatic savings of computational effort. A similar subsystem block density distribution analysis was used by Rovere et al. [25,26] to find the critical properties of a 2D Lennard–Jones (LJ) fluid, though was deemed to be less efficient than GEMC in this respect [27]. Our implementation differs from that of Rovere et al. in two aspects; we employ MD in place of Monte Carlo, and we attempt to avoid gathering data in the interfacial region by microstructural analysis of the fluid. We demonstrate the method on 3D off-lattice fluids and solids, exploring several types of phase coexistence. While the speed of single processor computers continues to grow exponentially, for the past decade the power of parallel computers has grown at a greater rate and the cost of parallel machines has dropped quickly with the proliferation of “Beowulf” type cluster systems [28]. Unlike Monte Carlo methods for molecular simulation, MD algorithms may be efficiently parallelized in several ways [8,29]. Calculations using spatial decompositions [30], such as the ones discussed in this work, are applicable to short-ranged potentials and complete in O(N/P) wall clock time, where N is the number of molecules in the system, and P is the number of processors used. Thus, for the rapid location of phase equilibria, MD methods have an advantage in that P can be increased in order to decrease the time before a solution is obtained. This fact has recently been exploited in several molecular simulation strategies to study thermophysical properties of fluids [8].
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
3
Fig. 1. Temperature T vs. density ρ diagram for a pure fluid. Point (a) corresponds to the initial one-phase system. During the simulation, the temperature is dropped to a lower value, point (b), placing the system in an unstable state. The fluid spontaneously separates into two phases in coexistence with each other; a dense liquid and a vapor.
2. Temperature quench molecular dynamics (TQMD) method 2.1. Method The basis of the method is summarized in Fig. 1. A single component liquid–vapor system is considered in the Figure, but the method is entirely general. One starts with a one-phase system (point “a”) with a given number of particles N, volume V and temperature T. The temperature is controlled by means of a “thermostat” affecting the equations of motion. Keeping both N and V constant, the temperature of the system is lowered abruptly by adjusting the thermostat. The target temperature (point “b”) must be such that the resulting state point lies within the spinodal envelope, e.g. at conditions that are both mechanically and thermodynamically unstable. As the system relaxes at the new state point, domains of liquid and vapor form. The connectivity and morphology of these phases will depend on their volume fractions, given by the usual lever rule, and the spinodal decomposition process has been studied elsewhere [31–34]. In Fig. 2, we present a sequence of snapshots that illustrate this evolution. After quenching, the system is far from equilibrium, and the driving force for diffusive transport is maximized. At short times the surface area between the two phases is very large, and the combination of these two factors ensures that local densities and concentrations stabilize at their equilibrium values very quickly. At later times, the reduction of surface tension between the two phases is the driving force for the observed “coarsening”, which is thus rather slow. We also note that no data is gathered at the initial point “a” in simulations of this type, and therefore careful equilibration at this temperature is not necessary. The MD method amounts to a numerical integration of Newton’s equations of motion for each particle, or in the case of constant temperature and/or pressure simulations, the equations of motion for a system of particles coupled to external heat and/or volume “baths”. In the calculations performed here, the integration is accomplished with a fifth-order Gear predictor–corrector algorithm [35]. The temperature
4
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 2. Snapshots from a simulation of 8000 molecules of a LJ fluid with CS potential at rcut = 5σ . Initial state (corresponding to point (a) of Fig. 1) is at a temperature T ∗ = 4, ρ ∗ = 0.33. The system is quenched to T ∗ = 0.9. Snapshots (b–h) are taken at 5000 step intervals (t ∗ = 0.004), showing the evolution of the system. The last snapshot (i), corresponds to an equilibrated state at 80 000 time steps showing a dense liquid region (the top and bottom of the simulation box) which coexists with an equilibrium vapor phase (a center slab of the simulation box).
is controlled with the Nosé–Hoover [36,37] thermostat, which ensures a canonical distribution of momenta under equilibrium conditions. For isothermal–isobaric (NPT) simulations, this approach is extended to also couple the volume of the simulation cell to a second external degree of freedom, which drives the system towards an imposed constant pressure [37]. 2.2. Parallel implementation In order to test the TQMD method under the conditions envisioned for its “real-world” application, we have implemented MD simulations on clustered workstations using a spatial–decomposition approach, in which a rectilinear volume is assigned to each processor using a 3D decomposition scheme. Communications between processors are handled using the message passing interface (MPI) library. Each processor is responsible for the storage of data pertaining to particles within its assigned volume, for the integration of the equations of motion for these particles, for the computation of the forces these particles
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
5
exert on each other, and for computation of one-half of the interactions (forces) between its own particles and those in neighboring domains. This approach is well documented [8,36]. At the beginning of each time-step, each processor updates its own particles’ positions using the first half of the predictor–corrector integration algorithm. The positions of particles that are within potential function range of neighboring domains are then packed into arrays and sent to those processors; particles near to edges and corners of a domain may be sent to more than one other processor. Upon receipt of all positions, each processor then computes interactions between its own particles and those which it has just been sent. In the next phase, these “partial force vectors” are sent back to the processors which own the relevant particles, where they are summed to obtain the complete force on every particle in the system. The corrector step of the integration algorithm is then made, completing the time-step. The potential functions used here are relatively short-ranged, so communication is only necessary between processors with immediately adjacent domains. This arrangement theoretically leads to perfectly linear scaling, in that if the simulation cell size is increased at the same rate as processors are added, the total wall clock time per time-step does not change. In practice, the use of thermostatted and/or baristatted equations of motion interferes with this scaling, because these integration algorithms require the calculation of the kinetic energy and virial of the entire system. This requires an “all-to-all” type of communication, implemented by calling the MPI ALLREDUCE function of the MPI library, which collects information from every processor using a tree algorithm that scales as the logarithm of the number of processors. We have found that the wait time associated with this operation is entirely negligible for calculations on as many as 216 processors. In addition, the volume decomposition approach leads to excellent scaling in the use of memory, and at no time does our code require the storage of the complete simulation cell on a single processor; this allows for the simulation of millions of atoms even on systems where the per-node memory is small. There is a single supralinear scaling in the use of memory in our code. Link-cell algorithms are used by the force computation routines to cut down the number of particle–particle interactions which must be considered; in the current code implementation the link-cell arrays span the entire system rather than just a single domain, so that this use of memory on each processor scales with the system size. For modern computers and systems smaller than about 100 nm spatial dimension, this extra cost in memory is small and can be ignored. 2.3. Interface detection For the generation of probability distributions of density and mole fraction, the simulation cell is divided into cubic sub-cells, in which the mole fractions of each species and the total density are calculated and accumulated into histograms. These cells must be large enough that they can give a reasonable estimate of the density of a single phase; if the cells are too small, the density histogram will be “quantized” due to the small integer number of particles that can fit in each sub-cell. For relatively small simulation cells, and for systems in which phase separation occurs relatively slowly, these sub-cells may contain significant portions of two or more phases, which leads to “smearing out” of the histogram peaks corresponding to each phase. We propose an algorithm for detecting and avoiding this situation based on a microscopic analysis of the fluid configuration, so as to only collect histogram data in sub-cells which contain entirely one phase. In the liquid–vapor case, the algorithm is as follows. Particles which are “in” an interface between two phases have coordination numbers reflecting the interfacial region; e.g. they have fewer neighbors than particles in the liquid phase, and more than particles in the vapor phase. By counting the local coordination of every particle, these can be identified unambiguously. Sub-cells which contain
6
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
more than 30% “interfacial” particles are then (arbitrarily) excluded from the histogram count; this is a conservative estimate, and is found to not alter the shape or position of the resulting histogram peaks. In the case of discrimination by mole fraction, a similar algorithm can be applied to a single species; for every particle, one counts the total number of neighbors of a single species, and locates particles on the interface by values of this quantity that fall between those of the pure phases. For more complex systems, other order parameters may be developed. In all cases, it is necessary only to have a rough estimate of the density or concentration differences between the two phases, which is easily obtained through computer graphics visualizations of the quenched system. It is important to note that this detection procedure and the data collection in this method does not affect the MD trajectory in any way, and can be done after the simulation has completed. Furthermore, the cost of this analysis is low, and can be easily repeated many times using different detection methods, interface definitions, etc. without increasing the total cost of the calculation. In other studies of vapor–liquid coexistence, we have determined that this is a reliable algorithm for interface location, provided that sensible definitions of the coordination sphere are used. In Fig. 3, we show a typical result from a system described below. The simulation cell, containing 16 000 particles, is divided into (9)3 sub-cells and densities are calculated in each sub-cell. If one includes sub-cells containing interfaces, the interpretation of the resulting histogram is difficult. Once the interfacial detection algorithm is run, the “noise” is reduced and a clearer picture of the coexisting densities is obtained. In the real application of the method, several configurations can be analyzed, greatly improving the statistical quality of these histograms.
Fig. 3. Histogram of the occurrence of given densities within a simulation box after equilibration. The system is a pure LJ fluid with CS potential at rcut = 5σ at T ∗ = 1.15 with overall density of ρ ∗ = 0.328 using N = 16 000 particles. The vertical scale is the relative occurrence probability, f, of each density. Simple histogramming of the density in cubic sub-cells of four on each side, gives the “raw” data (dashed line), while the use of the interface detection algorithm gives a density distribution characteristic of two-phase coexistence (solid line). Vertical hairlines note the equilibrium densities for this particular configuration, corresponding to ρ ∗ = 0.08 ± 0.02 (vapor phase) and ρ ∗ = 0.60 ± 0.02 (liquid phase).
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
7
2.4. Intermolecular potentials We have tested this method on the LJ fluid and binary LJ mixtures. The LJ potential is described by 6 σij 12 σij Uij (rij ) = 4εij − (1) rij rij where rij is the center–center intermolecular distance and σ and ε are the size and energy parameters. In this case, the Lorentz–Berthelot mixing rules are employed for the unlike-pair interactions, namely √ εij = εi εj and σij = (σi +σj )/2 where i and j refer to the pure component values. The potential and other thermodynamic properties, such as temperature T, pressure P, and density ρ, may be made dimensionless by taking U ∗ = U /ε1 , r ∗ = r/σ 1 , T ∗ = kT/ε1 , P ∗ = P σ13 /ε1 and ρ ∗ = ρσ13 where the subscript 1 refers to one component of the mixture, or to the single component in the case of a simple fluid. For computer simulations, due to the finite size of the simulation cells, the full potential is often truncated at a given cutoff radius, rcut , i.e. artificially made equal to zero at distances rij > rcut . In MD simulations, the potential is truncated at the cutoff distance and then shifted with the purpose of eliminating an impulse force at the truncation radius. The so called “cut and shifted” (CS) potential takes the form of Uij (r) − Uij (rcut ), for r ≤ rcut CS Uij (rij ) = (2) 0, for r > rcut There is an implicit assumption here that the forces affecting a central particle due to the effect of the forces from the particles outside the cut-off radius average out, and therefore have no effect on the central particle’s trajectory. In an inhomogeneous fluid, or for particles close to an interface, this is no longer true and additional forces must be considered [38,39]. Unless otherwise noted, here we will use a CS potential with an rcut of 5σ . 2.5. Simulation details After the initial temperature √ quenching, local equilibration is attained rapidly, usually within 10 000 time steps of t ∗ = t/(σ m/ε) = 0.004, where m is the mass of a molecule. From this point on, results are averaged from at least five different configurations, each spaced by 1000 time steps. In an effort to improve the statistics without increasing the run-length, and to insure that sizeable regions of “bulk-like” liquid and vapor are present, relatively large system sizes are employed. Good results were obtained with N = 8000 for single component systems and N = 32 000 for multi-component systems. The higher the temperature, the larger the system size which is needed, since the liquid–vapor interface becomes broad at high temperature, reducing the available bulk phase information in the system. Even though we have used cubic simulation cells with periodic boundary conditions in all Cartesian dimensions, one may opt to use a cell with one side smaller than the other two. Since upon equilibration, the system attempts to minimize its free energy by reducing the interfacial area, one is guaranteed that the interface will be planar (given that the overall density lies well between both bulk densities) and parallel to the Cartesian directors of the system box, (c.f. Fig. 2). A “skinny” box, in which one dimension is significantly smaller than the other two, has the advantage of reducing the size of the interface and maximizing the volume of the bulk phases, thus improving the quality of the density distribution histograms. Close to the critical point,
8
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 4. Influence of system size (N) on the estimated saturated liquid phase density ρ ∗ for a pure LJ fluid with CS potential at rcut = 5σ at T ∗ = 1.15. Dashed line is the assumed result.
however, this apparent advantage turns out to be an unfavorable choice, since long-range correlations are related to the smallest dimension of the simulation cell. Fig. 4 shows the effect of the size of the system on the resulting liquid densities. Systems sizes of the order of O(N) = 105 may be needed in extreme cases, which still can be completed in a reasonable time if the computation is done with a parallel machine. The calculations here have all been performed on a cluster consisting of 20 dual processor Pentium 550 MHz machines connected via a 1 GB/s Myrinet switch. Most calculations were made using only four “worker” processors and one “master” processor for input/output. 3. Results Fig. 5 and Table 1 show the calculated liquid–vapor data of a LJ fluid with a cut and shifted potential at a distance of 5σ using N = 32 000 and 64 000 for the higher temperatures. Histograms are calculated with a 0.01 spacing in density, and the subsystems are taken to be cubic boxes of side 5σ and 6σ for the liquid and vapor data points. Certainly, the choice of sub-cell size, although arbitrary, has an effect on the result. The smaller the sub-cell, the larger the dispersion of the histogram. The larger the sub-cell, the more likely it is to include “interfacial” molecules and thus, the whole sub-cell (and its constituent molecules) is excluded from the analysis. We have used a sub-cell size which gives “quantization” of the density of order 0.01, which is reasonably fine-grained. The system is explored up to densities of ρ ∗ = 0.86 without difficulty. The vapor densities at low temperatures are obscured due to the non-zero discretization of the histogram. The results follow the
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
9
Fig. 5. Temperature–density diagram for the LJ fluid. Solid circles are TQMD results for a CS potential at rcut = 5σ , data in Table 1. Solid line is the equation of state of Johnson et al. [40] for the full potential (Eq. (1)).
expected equation of state of Johnson et al. [40] for the full potential (since the cutoff radius is relatively large). In Figs. 6 and 7 and Table 2, we show the results of the vapor–liquid equilibria of a binary system. Our test mixture, which has been studied by other methods, consisted of equal sized LJ spheres (σ 2 /σ1 = 1) with ε2 /ε1 = 0.75. The results are compared with GEMC simulations of Harismiadis et al. [41], the “NPT + test particle” simulations of Lofti et al. [42], and the equation of state of Johnson et al. [40]. In Table 1 Densities ρ ∗ of a cut and shifted (rcut = 5σ ) LJ fluid along the vapor–liquid equilibria curve T∗
ρ ∗ (vapor)
1.23 1.2 1.16 1.15 1.1 1.0 0.9 0.8 0.7 0.65
0.125 ± 0.03 0.110 ± 0.025 0.070 ± 0.01 0.080 ± 0.01 0.060 ± 0.01 0.025 ± 0.01 0.010 ± 0.01 0.005
ρ ∗ (liquid) 0.58 ± 0.03 0.61 ± 0.02 0.62 ± 0.02 0.64 ± 0.01 0.69 ± 0.01 0.75 ± 0.01 0.79 ± 0.01 0.83 ± 0.01 0.85 ± 0.01
10
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 6. Temperature–composition diagram for a LJ mixture with σ 2 /σ1 = 1 and ε 2 /ε1 = 0.75, data in Table 2. Open squares are GEMC simulations of Harismiadis et al. [41], open circles are “NPT + test particle” simulations of Lofti et al. [42], closed circles are TQMD. The solid line is the equation of state of Johnson et al. [40].
both the pressure–composition and pressure–density diagrams the TQMD results are in good agreement with previously published data and are of equivalent quality. It appears that TQMD is not limited to fluid phase equilibria. If the final temperature is below the triple point, solid phases can nucleate during the quenching process. Table 3 shows the results of the Table 2 Vapor–liquid equilibrium data for a LJ mixture with σ 2 /σ1 = 1 and ε2 /ε1 = 0.75 T∗
x1
y1
ρ ∗ (vapor)
ρ ∗ (liquid)
0.027 0.03 0.04 0.044 0.058 0.069 0.084
1.00 ± 0.005 0.87 ± 0.003 0.56 ± 0.002 0.45 ± 0.002 0.30 ± 0.003 0.21 ± 0.003 0.11 ± 0.002
0.94 ± 0.003 0.80 ± 0.002 0.72 ± 0.002 0.53 ± 0.003 0.38 ± 0.003 0.20 ± 0.002
0.04 ± 0.005 0.04 ± 0.003 0.04 ± 0.002 0.05 ± 0.002 0.08 ± 0.003 0.10 ± 0.003 0.13 ± 0.002
0.69 ± 0.005 0.67 ± 0.003 0.65 ± 0.002 0.64 ± 0.002 0.62 ± 0.003 0.56 ± 0.003 0.48 ± 0.002
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
11
Fig. 7. Temperature–density diagram for a LJ mixture. Conditions and symbols as in Fig. 6.
dense-phase equilibrium densities along the sublimation line for a LJ fluid CS potential at rcut = 5σ . The coexisting vapor densities were too low (ρ ∗ < 0.01) to be accurately calculated without the use of excessively large simulation cells. For a LJ fluid, either a polycrystalline FCC structure or an amorphous solid may be attained. This latter is a metastable supercooled liquid that can be formed during the quenching procedure. In Fig. 8, we show a resulting configuration in which sections of imperfect FCC crystals coexist with a vapor phase, which for a simulation cell of that size, is almost void of molecules. In order to encourage growth of stable crystals, some care must be taken in the quenching Table 3 Solid densities ρ ∗ of a cut and shifted (rcut = 5σ ) LJ fluid along the sublimation (solid–vapor equilibria) curve T∗
ρ∗
Aggregation state
0.6 0.55 0.5 0.475 0.45 0.45 0.35
0.88 0.88 0.92 0.94 0.96 1.02 1.08
Amorphous Amorphous Amorphous Amorphous Amorphous FCC FCC
The errors in density are of the order of ±0.015 in reduced units. The amorphous states correspond to metastable sub-cooled liquid states.
12
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 8. Snapshot of a locally equilibrated configuration for a pure LJ fluid with a CS potential at rcut = 5σ . Temperature T∗ is 0.35. The initial (before quench) configuration is similar to that of Fig. 2a.
process, particularly performing the quenching in several steps and over long equilibration times. More extensive studies of the efficacy of this approach for locating solid–fluid equilibria will be made in future work.
4. Conclusions We present a method for the location of phase coexistence points based on canonical MD simulations coupled with a novel post-simulation analysis method. The method relies on the propensity of a system in a mechanically unstable state to rapidly separate into two coexisting phases separated by an interface. This interface, of course, is exactly what previously developed methods, such as GEMC are designed to avoid; the use of relatively large simulation cells, accessible to modern computing equipment, makes it possible to avoid the difficulties presented by interfaces in early work. Because local equilibration of densities and compositions occurs quickly, only short simulations are necessary; the extremely fast quenches accessible to simulation provide a large driving force for phase separation. The use of MD for these calculations enables both the use of massively parallel computers, since the efficient parallelization of MD is a solved problem, and the application of the method to molecular species of great size or
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
13
complexity; parallelized MD simulations of complex fluids, polymers, proteins, and DNA are now considered routine. This approach is particularly suited to the determination of phase equilibria in systems about which little is known. A sequence of quenches at different densities and quench temperatures could be used to locate points of phase coexistence quickly and robustly, and quenches which do not result in phase separation used to bound coexistence lines and surfaces. The generation of high temperature configurations with which to start the method is not difficult, and in any case, the properties of the initial configuration are of no real importance. The reliance of this method on large simulation cells allows a close approach to critical states, though for the precise location of critical points or lines, HR methods [3,43] are more suitable. Data gathered via TQMD can be fit with correlations or equations of state for reasonable location of critical points, and, should critical parameters be required to high precision, HR methods could then be applied. We note that previous applications of HR methods in VLE have been in systems where the location of the critical point is already known approximately; this information was used to start the simulations near the critical point, working down in temperature along the liquid and vapor branches of the coexistence curve [3]. We have shown that this method gives correct results for pure and multi-component vapor–liquid equilibria. In a preliminary application to solid–vapor equilibria, we observed that the solid phase could be accessed, at least in the case of simple fluids. For the simulation of simple fluids, this method uses more computer time than methods such as GEMC, but is easily parallelized and may thus be made to return in shorter real time. More importantly, this approach can be applied in a straightforward way to more types of phase equilibrium and more complex molecular species than any single competing method. List of symbols f relative occurrence probability k Boltzmann’s constant P pressure r intermolecular distance rcut cut off radius T temperature U intermolecular potential x liquid phase composition y vapor phase composition Greek letters ε Lennard–Jones potential energy parameter ρ density σ Lennard–Jones potential size parameter Subscripts 1 first component in a mixture i, j corresponding to molecules i and j Superscript ∗ reduced units
14
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Acknowledgements E.A.M. wishes to acknowledge the financial support given by the Agenda Petroleo (project 97-003530) to this work. Some of the calculations reported here were performed by Mar´ıa Eugenia Suárez. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813–826. A.Z. Panagiotopoulos, Mol. Simulation 9 (1992) 1–23. A.Z. Panagiotopoulos, J. Phys.: Condens. Matter 12 (2000) R25–R52. A.M. Ferrenberg, A.M. Swendsen, Phys. Rev. Lett. 61 (1988) 2635–2638. J.I. Siepmann, D. Frenkel, Mol. Phys. 75 (1992) 59–70. D. Frenkel, G.C.A.M. Mooij, B. Smit, J. Phys.: Condens. Matter 4 (1992) 3053–3076. J.J. de Pablo, M. Laso, U.W. Suter, J. Chem. Phys. 96 (1992) 2395–2403. G.S. Heffelfinger, Comp. Phys. Commun. 128 (2000) 219–237. L.D.J.C. Loyens, B. Smit, K. Esselink, Mol. Phys. 86 (1995) 171–183. D.A. Kofke, Mol. Phys. 78 (1993) 1331–1336. D.A. Kofke, J. Chem. Phys. 98 (1993) 4149–4162. D. Möller, J. Fischer, Mol. Phys. 69 (1990) 463–473. http://www.amber.ucsf.edu/amber/amber.html. http://yuri.harvard.edu/. http://www.dl.ac.uk/TCSC/Software/DL POLY/main.html. http://igc.ethz.ch/gromos/gromos.html. http://www.earth.ox.ac.uk/∼keith/moldy.html. http://www.gromacs.org/. C.D. Holcomb, P. Clancy, S.M. Thompson, J.A. Zollweg, Fluid Phase Equilib. 75 (1992) 185–196. C.D. Holcomb, P. Clancy, J.A. Zollweg, Mol. Phys. 78 (1993) 437–459. F. Goujon, P. Malfreyt, A. Boutin, A.H. Fuchs, Mol. Simulation 27 (2001) 99–114. J.S. Rowlinson, B. Widom, Molecular Theory of Capillarity, Clarendon Press, Oxford, 1989. M. Mecke, J. Winkelmann, M. Fischer, J. Chem. Phys. 110 (1999) 1188–1194. L.X. Dang, J. Chem. Phys. 110 (1999) 10113–10122. M. Rovere, W.D. Heermann, K. Binder, J. Phys.: Condens. Matter 2 (1990) 7009–7032. M. Rovere, P. Nielaba, K. Binder, Z. Phys. B 90 (1993) 215–228. B. Smit, Ph. de Smedt, D. Frenkel, Mol. Phys. 68 (1989) 931–950. http://www.beowulf.org/. S. Plimpton, J. Comp. Phys. 117 (1995) 1–19. D. Fincham, Mol. Simulation 1 (1987) 1–46. J.W. Cahn, J. Chem. Phys. 42 (1965) 93–99. M. Laradji, O.G. Mouritsen, S. Toxvaerd, Phys. Rev. E 53 (1996) 3672–3681. M. Laradji, S. Toxvaerd, O.G. Mouritsen, Phys. Rev. Lett. 77 (1996) 2253–2256. L.D. Gelb, K.E. Gubbins, Phys. Rev. E 56 (1997) 3185–3196. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. S. Nosé, Mol. Phys. 52 (1984) 255–268. W.G. Hoover, Phys. Rev. A 31 (1985) 1695–1697. A. Lofti, J. Vrabec, J. Fischer, Mol. Simulation 5 (1990) 233. M. Guo, B.C.Y. Lu, J. Chem. Phys. 106 (1997) 3685–3688. J.K. Johnson, J.A. Zollweg, K.E. Gubbins, Mol. Phys. 78 (1993) 591–618. V.I. Harismiadis, N.K. Koutras, D.P. Tassios, A.Z. Panagiotopoulos, Fluid Phase Equilib. 65 (1991) 1–18. A. Lofti, J. Vrabec, J. Fischer, Mol. Phys. 112 (1995) 173–197. F.A. Escobedo, J. Phys. Chem. 113 (2000) 8444–8456.
Location of phase equilibria by temperature-quench molecular dynamics simulations Lev D. Gelb a,1 , Erich A. Müller b,∗ b
a Department of Chemistry and Biochemistry, Florida State University, Tallahassee, FL 32306-4390, USA Departamento de Termodinámica y Fenómenos de Transferencia, Universidad Simón Bol´ıvar, Caracas 1080, Venezuela
Received 2 January 2002; accepted 14 May 2002
Abstract We present a method to locate phase coexistence points using molecular dynamics (MD) simulations. By equilibrating a single phase in a relatively large simulation cell under NVT conditions and then quenching the system into a two-phase region, it spontaneously separates into two coexisting phases. The different phases are identified by collecting histograms of the local density and/or local composition of the pure phase regions, which can be done soon after the quench. The method can be used to locate vapor–liquid, liquid–liquid or solid–fluid equilibria. MD is a robust simulation technique that can be applied without difficulty to very complex molecules and mixtures, and can be implemented on large parallel computers. This approach, thus, has the potential to be both faster (in wall clock time) and more generally applicable than any other single technique. The method is demonstrated on test systems of single component and binary Lennard–Jones (LJ) fluids. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Molecular simulation; Method of calculation; Vapor–liquid equilibrium; Solid–fluid equilibria
1. Introduction The most popular method for locating phase coexistence in molecular simulation is still Gibbs Ensemble Monte Carlo [1–3], though it is being slowly supplanted by histogram-reweighting (HR) techniques, which offer substantially more precise location of critical points [4]. Monte Carlo methods have three well-known deficiencies. When simulating dense phases, typical of strongly interacting liquids and solids, equilibration is difficult to achieve because of the poor statistics associated with the insertion/deletion steps. Secondly, they are difficult to apply to systems containing very complex molecules without substantial system-specific modifications (such as Configurational-Bias Monte Carlo [5–7]). ∗
Corresponding author. E-mail addresses: [email protected] (L.D. Gelb), [email protected] (E.A. Müller). 1 Present address: Department of Chemistry, Washington University , St. Louis, MO 63130-1134, USA. 0378-3812/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 3 8 1 2 ( 0 2 ) 0 0 1 7 4 - 7
2
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Thirdly, Monte Carlo methods are difficult to implement on parallel computers and few efforts have been made to improve this situation [8,9]. Kofke [10,11] suggested “Gibbs–Duhem integration”, an indirect method based on the integration of the Clausius–Clapeyron equation. However, to apply Gibbs–Duhem integration, one must start with a known phase coexistence point. The “NPT + test particle” method [12] may also be used to study phase equilibria, but suffers from the same limitations as other Monte Carlo simulations, since it relies on a particle insertion technique. In contrast, the molecular dynamics (MD) simulation method is applicable to both dilute and dense phases, easily parallelized, and widely used for simulating very large and complex molecules. Several parallelized MD codes are available, including AMBER [13], CHARMM [14], DL POLY [15], Gromos96 [16], Moldy [17], and Gromacs [18]. The application of MD to study phase equilibria is not new, and in fact predates GEMC. A typical implementation is that of Holcomb et al. [19,20], in which previously equilibrated fractions of a bulk liquid and bulk vapor phase are placed in a simulation box in the form of a slab of liquid surrounded by a vapor phase. The system is then allowed to evolve under NVT conditions until equilibrium is reached through diffusive mass transport. Goujon et al. [21] have applied a similar method, using MC instead of MD. This approach is still used to study interfacial properties [22–24], and while it can provide VLE or LLE coexistence data, it is not considered efficient because long times are required for complete equilibration of a two-phase system. In such systems, the last stage of equilibration consists of the minimization of the (already small) surface area between the two phases, which provides only a relatively weak driving force. What follows is the description of a method that allows the location of phase coexistence through a constant density simulation in which the temperature is changed in a single time-step (quenching) in order to place the system in a thermodynamically and mechanically unstable state. At this point, the system spontaneously separates into domains of stable, locally equilibrated phases. These phases are then detected and characterized by dividing the system into small blocks and analyzing the densities and/or the compositions in the blocks. Phase coexistence data is, therefore, gathered from the locally equilibrated system, and it is not necessary to continue the simulation until global equilibrium is reached, resulting in a dramatic savings of computational effort. A similar subsystem block density distribution analysis was used by Rovere et al. [25,26] to find the critical properties of a 2D Lennard–Jones (LJ) fluid, though was deemed to be less efficient than GEMC in this respect [27]. Our implementation differs from that of Rovere et al. in two aspects; we employ MD in place of Monte Carlo, and we attempt to avoid gathering data in the interfacial region by microstructural analysis of the fluid. We demonstrate the method on 3D off-lattice fluids and solids, exploring several types of phase coexistence. While the speed of single processor computers continues to grow exponentially, for the past decade the power of parallel computers has grown at a greater rate and the cost of parallel machines has dropped quickly with the proliferation of “Beowulf” type cluster systems [28]. Unlike Monte Carlo methods for molecular simulation, MD algorithms may be efficiently parallelized in several ways [8,29]. Calculations using spatial decompositions [30], such as the ones discussed in this work, are applicable to short-ranged potentials and complete in O(N/P) wall clock time, where N is the number of molecules in the system, and P is the number of processors used. Thus, for the rapid location of phase equilibria, MD methods have an advantage in that P can be increased in order to decrease the time before a solution is obtained. This fact has recently been exploited in several molecular simulation strategies to study thermophysical properties of fluids [8].
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
3
Fig. 1. Temperature T vs. density ρ diagram for a pure fluid. Point (a) corresponds to the initial one-phase system. During the simulation, the temperature is dropped to a lower value, point (b), placing the system in an unstable state. The fluid spontaneously separates into two phases in coexistence with each other; a dense liquid and a vapor.
2. Temperature quench molecular dynamics (TQMD) method 2.1. Method The basis of the method is summarized in Fig. 1. A single component liquid–vapor system is considered in the Figure, but the method is entirely general. One starts with a one-phase system (point “a”) with a given number of particles N, volume V and temperature T. The temperature is controlled by means of a “thermostat” affecting the equations of motion. Keeping both N and V constant, the temperature of the system is lowered abruptly by adjusting the thermostat. The target temperature (point “b”) must be such that the resulting state point lies within the spinodal envelope, e.g. at conditions that are both mechanically and thermodynamically unstable. As the system relaxes at the new state point, domains of liquid and vapor form. The connectivity and morphology of these phases will depend on their volume fractions, given by the usual lever rule, and the spinodal decomposition process has been studied elsewhere [31–34]. In Fig. 2, we present a sequence of snapshots that illustrate this evolution. After quenching, the system is far from equilibrium, and the driving force for diffusive transport is maximized. At short times the surface area between the two phases is very large, and the combination of these two factors ensures that local densities and concentrations stabilize at their equilibrium values very quickly. At later times, the reduction of surface tension between the two phases is the driving force for the observed “coarsening”, which is thus rather slow. We also note that no data is gathered at the initial point “a” in simulations of this type, and therefore careful equilibration at this temperature is not necessary. The MD method amounts to a numerical integration of Newton’s equations of motion for each particle, or in the case of constant temperature and/or pressure simulations, the equations of motion for a system of particles coupled to external heat and/or volume “baths”. In the calculations performed here, the integration is accomplished with a fifth-order Gear predictor–corrector algorithm [35]. The temperature
4
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 2. Snapshots from a simulation of 8000 molecules of a LJ fluid with CS potential at rcut = 5σ . Initial state (corresponding to point (a) of Fig. 1) is at a temperature T ∗ = 4, ρ ∗ = 0.33. The system is quenched to T ∗ = 0.9. Snapshots (b–h) are taken at 5000 step intervals (t ∗ = 0.004), showing the evolution of the system. The last snapshot (i), corresponds to an equilibrated state at 80 000 time steps showing a dense liquid region (the top and bottom of the simulation box) which coexists with an equilibrium vapor phase (a center slab of the simulation box).
is controlled with the Nosé–Hoover [36,37] thermostat, which ensures a canonical distribution of momenta under equilibrium conditions. For isothermal–isobaric (NPT) simulations, this approach is extended to also couple the volume of the simulation cell to a second external degree of freedom, which drives the system towards an imposed constant pressure [37]. 2.2. Parallel implementation In order to test the TQMD method under the conditions envisioned for its “real-world” application, we have implemented MD simulations on clustered workstations using a spatial–decomposition approach, in which a rectilinear volume is assigned to each processor using a 3D decomposition scheme. Communications between processors are handled using the message passing interface (MPI) library. Each processor is responsible for the storage of data pertaining to particles within its assigned volume, for the integration of the equations of motion for these particles, for the computation of the forces these particles
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
5
exert on each other, and for computation of one-half of the interactions (forces) between its own particles and those in neighboring domains. This approach is well documented [8,36]. At the beginning of each time-step, each processor updates its own particles’ positions using the first half of the predictor–corrector integration algorithm. The positions of particles that are within potential function range of neighboring domains are then packed into arrays and sent to those processors; particles near to edges and corners of a domain may be sent to more than one other processor. Upon receipt of all positions, each processor then computes interactions between its own particles and those which it has just been sent. In the next phase, these “partial force vectors” are sent back to the processors which own the relevant particles, where they are summed to obtain the complete force on every particle in the system. The corrector step of the integration algorithm is then made, completing the time-step. The potential functions used here are relatively short-ranged, so communication is only necessary between processors with immediately adjacent domains. This arrangement theoretically leads to perfectly linear scaling, in that if the simulation cell size is increased at the same rate as processors are added, the total wall clock time per time-step does not change. In practice, the use of thermostatted and/or baristatted equations of motion interferes with this scaling, because these integration algorithms require the calculation of the kinetic energy and virial of the entire system. This requires an “all-to-all” type of communication, implemented by calling the MPI ALLREDUCE function of the MPI library, which collects information from every processor using a tree algorithm that scales as the logarithm of the number of processors. We have found that the wait time associated with this operation is entirely negligible for calculations on as many as 216 processors. In addition, the volume decomposition approach leads to excellent scaling in the use of memory, and at no time does our code require the storage of the complete simulation cell on a single processor; this allows for the simulation of millions of atoms even on systems where the per-node memory is small. There is a single supralinear scaling in the use of memory in our code. Link-cell algorithms are used by the force computation routines to cut down the number of particle–particle interactions which must be considered; in the current code implementation the link-cell arrays span the entire system rather than just a single domain, so that this use of memory on each processor scales with the system size. For modern computers and systems smaller than about 100 nm spatial dimension, this extra cost in memory is small and can be ignored. 2.3. Interface detection For the generation of probability distributions of density and mole fraction, the simulation cell is divided into cubic sub-cells, in which the mole fractions of each species and the total density are calculated and accumulated into histograms. These cells must be large enough that they can give a reasonable estimate of the density of a single phase; if the cells are too small, the density histogram will be “quantized” due to the small integer number of particles that can fit in each sub-cell. For relatively small simulation cells, and for systems in which phase separation occurs relatively slowly, these sub-cells may contain significant portions of two or more phases, which leads to “smearing out” of the histogram peaks corresponding to each phase. We propose an algorithm for detecting and avoiding this situation based on a microscopic analysis of the fluid configuration, so as to only collect histogram data in sub-cells which contain entirely one phase. In the liquid–vapor case, the algorithm is as follows. Particles which are “in” an interface between two phases have coordination numbers reflecting the interfacial region; e.g. they have fewer neighbors than particles in the liquid phase, and more than particles in the vapor phase. By counting the local coordination of every particle, these can be identified unambiguously. Sub-cells which contain
6
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
more than 30% “interfacial” particles are then (arbitrarily) excluded from the histogram count; this is a conservative estimate, and is found to not alter the shape or position of the resulting histogram peaks. In the case of discrimination by mole fraction, a similar algorithm can be applied to a single species; for every particle, one counts the total number of neighbors of a single species, and locates particles on the interface by values of this quantity that fall between those of the pure phases. For more complex systems, other order parameters may be developed. In all cases, it is necessary only to have a rough estimate of the density or concentration differences between the two phases, which is easily obtained through computer graphics visualizations of the quenched system. It is important to note that this detection procedure and the data collection in this method does not affect the MD trajectory in any way, and can be done after the simulation has completed. Furthermore, the cost of this analysis is low, and can be easily repeated many times using different detection methods, interface definitions, etc. without increasing the total cost of the calculation. In other studies of vapor–liquid coexistence, we have determined that this is a reliable algorithm for interface location, provided that sensible definitions of the coordination sphere are used. In Fig. 3, we show a typical result from a system described below. The simulation cell, containing 16 000 particles, is divided into (9)3 sub-cells and densities are calculated in each sub-cell. If one includes sub-cells containing interfaces, the interpretation of the resulting histogram is difficult. Once the interfacial detection algorithm is run, the “noise” is reduced and a clearer picture of the coexisting densities is obtained. In the real application of the method, several configurations can be analyzed, greatly improving the statistical quality of these histograms.
Fig. 3. Histogram of the occurrence of given densities within a simulation box after equilibration. The system is a pure LJ fluid with CS potential at rcut = 5σ at T ∗ = 1.15 with overall density of ρ ∗ = 0.328 using N = 16 000 particles. The vertical scale is the relative occurrence probability, f, of each density. Simple histogramming of the density in cubic sub-cells of four on each side, gives the “raw” data (dashed line), while the use of the interface detection algorithm gives a density distribution characteristic of two-phase coexistence (solid line). Vertical hairlines note the equilibrium densities for this particular configuration, corresponding to ρ ∗ = 0.08 ± 0.02 (vapor phase) and ρ ∗ = 0.60 ± 0.02 (liquid phase).
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
7
2.4. Intermolecular potentials We have tested this method on the LJ fluid and binary LJ mixtures. The LJ potential is described by 6 σij 12 σij Uij (rij ) = 4εij − (1) rij rij where rij is the center–center intermolecular distance and σ and ε are the size and energy parameters. In this case, the Lorentz–Berthelot mixing rules are employed for the unlike-pair interactions, namely √ εij = εi εj and σij = (σi +σj )/2 where i and j refer to the pure component values. The potential and other thermodynamic properties, such as temperature T, pressure P, and density ρ, may be made dimensionless by taking U ∗ = U /ε1 , r ∗ = r/σ 1 , T ∗ = kT/ε1 , P ∗ = P σ13 /ε1 and ρ ∗ = ρσ13 where the subscript 1 refers to one component of the mixture, or to the single component in the case of a simple fluid. For computer simulations, due to the finite size of the simulation cells, the full potential is often truncated at a given cutoff radius, rcut , i.e. artificially made equal to zero at distances rij > rcut . In MD simulations, the potential is truncated at the cutoff distance and then shifted with the purpose of eliminating an impulse force at the truncation radius. The so called “cut and shifted” (CS) potential takes the form of Uij (r) − Uij (rcut ), for r ≤ rcut CS Uij (rij ) = (2) 0, for r > rcut There is an implicit assumption here that the forces affecting a central particle due to the effect of the forces from the particles outside the cut-off radius average out, and therefore have no effect on the central particle’s trajectory. In an inhomogeneous fluid, or for particles close to an interface, this is no longer true and additional forces must be considered [38,39]. Unless otherwise noted, here we will use a CS potential with an rcut of 5σ . 2.5. Simulation details After the initial temperature √ quenching, local equilibration is attained rapidly, usually within 10 000 time steps of t ∗ = t/(σ m/ε) = 0.004, where m is the mass of a molecule. From this point on, results are averaged from at least five different configurations, each spaced by 1000 time steps. In an effort to improve the statistics without increasing the run-length, and to insure that sizeable regions of “bulk-like” liquid and vapor are present, relatively large system sizes are employed. Good results were obtained with N = 8000 for single component systems and N = 32 000 for multi-component systems. The higher the temperature, the larger the system size which is needed, since the liquid–vapor interface becomes broad at high temperature, reducing the available bulk phase information in the system. Even though we have used cubic simulation cells with periodic boundary conditions in all Cartesian dimensions, one may opt to use a cell with one side smaller than the other two. Since upon equilibration, the system attempts to minimize its free energy by reducing the interfacial area, one is guaranteed that the interface will be planar (given that the overall density lies well between both bulk densities) and parallel to the Cartesian directors of the system box, (c.f. Fig. 2). A “skinny” box, in which one dimension is significantly smaller than the other two, has the advantage of reducing the size of the interface and maximizing the volume of the bulk phases, thus improving the quality of the density distribution histograms. Close to the critical point,
8
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 4. Influence of system size (N) on the estimated saturated liquid phase density ρ ∗ for a pure LJ fluid with CS potential at rcut = 5σ at T ∗ = 1.15. Dashed line is the assumed result.
however, this apparent advantage turns out to be an unfavorable choice, since long-range correlations are related to the smallest dimension of the simulation cell. Fig. 4 shows the effect of the size of the system on the resulting liquid densities. Systems sizes of the order of O(N) = 105 may be needed in extreme cases, which still can be completed in a reasonable time if the computation is done with a parallel machine. The calculations here have all been performed on a cluster consisting of 20 dual processor Pentium 550 MHz machines connected via a 1 GB/s Myrinet switch. Most calculations were made using only four “worker” processors and one “master” processor for input/output. 3. Results Fig. 5 and Table 1 show the calculated liquid–vapor data of a LJ fluid with a cut and shifted potential at a distance of 5σ using N = 32 000 and 64 000 for the higher temperatures. Histograms are calculated with a 0.01 spacing in density, and the subsystems are taken to be cubic boxes of side 5σ and 6σ for the liquid and vapor data points. Certainly, the choice of sub-cell size, although arbitrary, has an effect on the result. The smaller the sub-cell, the larger the dispersion of the histogram. The larger the sub-cell, the more likely it is to include “interfacial” molecules and thus, the whole sub-cell (and its constituent molecules) is excluded from the analysis. We have used a sub-cell size which gives “quantization” of the density of order 0.01, which is reasonably fine-grained. The system is explored up to densities of ρ ∗ = 0.86 without difficulty. The vapor densities at low temperatures are obscured due to the non-zero discretization of the histogram. The results follow the
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
9
Fig. 5. Temperature–density diagram for the LJ fluid. Solid circles are TQMD results for a CS potential at rcut = 5σ , data in Table 1. Solid line is the equation of state of Johnson et al. [40] for the full potential (Eq. (1)).
expected equation of state of Johnson et al. [40] for the full potential (since the cutoff radius is relatively large). In Figs. 6 and 7 and Table 2, we show the results of the vapor–liquid equilibria of a binary system. Our test mixture, which has been studied by other methods, consisted of equal sized LJ spheres (σ 2 /σ1 = 1) with ε2 /ε1 = 0.75. The results are compared with GEMC simulations of Harismiadis et al. [41], the “NPT + test particle” simulations of Lofti et al. [42], and the equation of state of Johnson et al. [40]. In Table 1 Densities ρ ∗ of a cut and shifted (rcut = 5σ ) LJ fluid along the vapor–liquid equilibria curve T∗
ρ ∗ (vapor)
1.23 1.2 1.16 1.15 1.1 1.0 0.9 0.8 0.7 0.65
0.125 ± 0.03 0.110 ± 0.025 0.070 ± 0.01 0.080 ± 0.01 0.060 ± 0.01 0.025 ± 0.01 0.010 ± 0.01 0.005
ρ ∗ (liquid) 0.58 ± 0.03 0.61 ± 0.02 0.62 ± 0.02 0.64 ± 0.01 0.69 ± 0.01 0.75 ± 0.01 0.79 ± 0.01 0.83 ± 0.01 0.85 ± 0.01
10
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 6. Temperature–composition diagram for a LJ mixture with σ 2 /σ1 = 1 and ε 2 /ε1 = 0.75, data in Table 2. Open squares are GEMC simulations of Harismiadis et al. [41], open circles are “NPT + test particle” simulations of Lofti et al. [42], closed circles are TQMD. The solid line is the equation of state of Johnson et al. [40].
both the pressure–composition and pressure–density diagrams the TQMD results are in good agreement with previously published data and are of equivalent quality. It appears that TQMD is not limited to fluid phase equilibria. If the final temperature is below the triple point, solid phases can nucleate during the quenching process. Table 3 shows the results of the Table 2 Vapor–liquid equilibrium data for a LJ mixture with σ 2 /σ1 = 1 and ε2 /ε1 = 0.75 T∗
x1
y1
ρ ∗ (vapor)
ρ ∗ (liquid)
0.027 0.03 0.04 0.044 0.058 0.069 0.084
1.00 ± 0.005 0.87 ± 0.003 0.56 ± 0.002 0.45 ± 0.002 0.30 ± 0.003 0.21 ± 0.003 0.11 ± 0.002
0.94 ± 0.003 0.80 ± 0.002 0.72 ± 0.002 0.53 ± 0.003 0.38 ± 0.003 0.20 ± 0.002
0.04 ± 0.005 0.04 ± 0.003 0.04 ± 0.002 0.05 ± 0.002 0.08 ± 0.003 0.10 ± 0.003 0.13 ± 0.002
0.69 ± 0.005 0.67 ± 0.003 0.65 ± 0.002 0.64 ± 0.002 0.62 ± 0.003 0.56 ± 0.003 0.48 ± 0.002
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
11
Fig. 7. Temperature–density diagram for a LJ mixture. Conditions and symbols as in Fig. 6.
dense-phase equilibrium densities along the sublimation line for a LJ fluid CS potential at rcut = 5σ . The coexisting vapor densities were too low (ρ ∗ < 0.01) to be accurately calculated without the use of excessively large simulation cells. For a LJ fluid, either a polycrystalline FCC structure or an amorphous solid may be attained. This latter is a metastable supercooled liquid that can be formed during the quenching procedure. In Fig. 8, we show a resulting configuration in which sections of imperfect FCC crystals coexist with a vapor phase, which for a simulation cell of that size, is almost void of molecules. In order to encourage growth of stable crystals, some care must be taken in the quenching Table 3 Solid densities ρ ∗ of a cut and shifted (rcut = 5σ ) LJ fluid along the sublimation (solid–vapor equilibria) curve T∗
ρ∗
Aggregation state
0.6 0.55 0.5 0.475 0.45 0.45 0.35
0.88 0.88 0.92 0.94 0.96 1.02 1.08
Amorphous Amorphous Amorphous Amorphous Amorphous FCC FCC
The errors in density are of the order of ±0.015 in reduced units. The amorphous states correspond to metastable sub-cooled liquid states.
12
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Fig. 8. Snapshot of a locally equilibrated configuration for a pure LJ fluid with a CS potential at rcut = 5σ . Temperature T∗ is 0.35. The initial (before quench) configuration is similar to that of Fig. 2a.
process, particularly performing the quenching in several steps and over long equilibration times. More extensive studies of the efficacy of this approach for locating solid–fluid equilibria will be made in future work.
4. Conclusions We present a method for the location of phase coexistence points based on canonical MD simulations coupled with a novel post-simulation analysis method. The method relies on the propensity of a system in a mechanically unstable state to rapidly separate into two coexisting phases separated by an interface. This interface, of course, is exactly what previously developed methods, such as GEMC are designed to avoid; the use of relatively large simulation cells, accessible to modern computing equipment, makes it possible to avoid the difficulties presented by interfaces in early work. Because local equilibration of densities and compositions occurs quickly, only short simulations are necessary; the extremely fast quenches accessible to simulation provide a large driving force for phase separation. The use of MD for these calculations enables both the use of massively parallel computers, since the efficient parallelization of MD is a solved problem, and the application of the method to molecular species of great size or
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
13
complexity; parallelized MD simulations of complex fluids, polymers, proteins, and DNA are now considered routine. This approach is particularly suited to the determination of phase equilibria in systems about which little is known. A sequence of quenches at different densities and quench temperatures could be used to locate points of phase coexistence quickly and robustly, and quenches which do not result in phase separation used to bound coexistence lines and surfaces. The generation of high temperature configurations with which to start the method is not difficult, and in any case, the properties of the initial configuration are of no real importance. The reliance of this method on large simulation cells allows a close approach to critical states, though for the precise location of critical points or lines, HR methods [3,43] are more suitable. Data gathered via TQMD can be fit with correlations or equations of state for reasonable location of critical points, and, should critical parameters be required to high precision, HR methods could then be applied. We note that previous applications of HR methods in VLE have been in systems where the location of the critical point is already known approximately; this information was used to start the simulations near the critical point, working down in temperature along the liquid and vapor branches of the coexistence curve [3]. We have shown that this method gives correct results for pure and multi-component vapor–liquid equilibria. In a preliminary application to solid–vapor equilibria, we observed that the solid phase could be accessed, at least in the case of simple fluids. For the simulation of simple fluids, this method uses more computer time than methods such as GEMC, but is easily parallelized and may thus be made to return in shorter real time. More importantly, this approach can be applied in a straightforward way to more types of phase equilibrium and more complex molecular species than any single competing method. List of symbols f relative occurrence probability k Boltzmann’s constant P pressure r intermolecular distance rcut cut off radius T temperature U intermolecular potential x liquid phase composition y vapor phase composition Greek letters ε Lennard–Jones potential energy parameter ρ density σ Lennard–Jones potential size parameter Subscripts 1 first component in a mixture i, j corresponding to molecules i and j Superscript ∗ reduced units
14
L.D. Gelb, E.A. Müller / Fluid Phase Equilibria 203 (2002) 1–14
Acknowledgements E.A.M. wishes to acknowledge the financial support given by the Agenda Petroleo (project 97-003530) to this work. Some of the calculations reported here were performed by Mar´ıa Eugenia Suárez. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813–826. A.Z. Panagiotopoulos, Mol. Simulation 9 (1992) 1–23. A.Z. Panagiotopoulos, J. Phys.: Condens. Matter 12 (2000) R25–R52. A.M. Ferrenberg, A.M. Swendsen, Phys. Rev. Lett. 61 (1988) 2635–2638. J.I. Siepmann, D. Frenkel, Mol. Phys. 75 (1992) 59–70. D. Frenkel, G.C.A.M. Mooij, B. Smit, J. Phys.: Condens. Matter 4 (1992) 3053–3076. J.J. de Pablo, M. Laso, U.W. Suter, J. Chem. Phys. 96 (1992) 2395–2403. G.S. Heffelfinger, Comp. Phys. Commun. 128 (2000) 219–237. L.D.J.C. Loyens, B. Smit, K. Esselink, Mol. Phys. 86 (1995) 171–183. D.A. Kofke, Mol. Phys. 78 (1993) 1331–1336. D.A. Kofke, J. Chem. Phys. 98 (1993) 4149–4162. D. Möller, J. Fischer, Mol. Phys. 69 (1990) 463–473. http://www.amber.ucsf.edu/amber/amber.html. http://yuri.harvard.edu/. http://www.dl.ac.uk/TCSC/Software/DL POLY/main.html. http://igc.ethz.ch/gromos/gromos.html. http://www.earth.ox.ac.uk/∼keith/moldy.html. http://www.gromacs.org/. C.D. Holcomb, P. Clancy, S.M. Thompson, J.A. Zollweg, Fluid Phase Equilib. 75 (1992) 185–196. C.D. Holcomb, P. Clancy, J.A. Zollweg, Mol. Phys. 78 (1993) 437–459. F. Goujon, P. Malfreyt, A. Boutin, A.H. Fuchs, Mol. Simulation 27 (2001) 99–114. J.S. Rowlinson, B. Widom, Molecular Theory of Capillarity, Clarendon Press, Oxford, 1989. M. Mecke, J. Winkelmann, M. Fischer, J. Chem. Phys. 110 (1999) 1188–1194. L.X. Dang, J. Chem. Phys. 110 (1999) 10113–10122. M. Rovere, W.D. Heermann, K. Binder, J. Phys.: Condens. Matter 2 (1990) 7009–7032. M. Rovere, P. Nielaba, K. Binder, Z. Phys. B 90 (1993) 215–228. B. Smit, Ph. de Smedt, D. Frenkel, Mol. Phys. 68 (1989) 931–950. http://www.beowulf.org/. S. Plimpton, J. Comp. Phys. 117 (1995) 1–19. D. Fincham, Mol. Simulation 1 (1987) 1–46. J.W. Cahn, J. Chem. Phys. 42 (1965) 93–99. M. Laradji, O.G. Mouritsen, S. Toxvaerd, Phys. Rev. E 53 (1996) 3672–3681. M. Laradji, S. Toxvaerd, O.G. Mouritsen, Phys. Rev. Lett. 77 (1996) 2253–2256. L.D. Gelb, K.E. Gubbins, Phys. Rev. E 56 (1997) 3185–3196. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. S. Nosé, Mol. Phys. 52 (1984) 255–268. W.G. Hoover, Phys. Rev. A 31 (1985) 1695–1697. A. Lofti, J. Vrabec, J. Fischer, Mol. Simulation 5 (1990) 233. M. Guo, B.C.Y. Lu, J. Chem. Phys. 106 (1997) 3685–3688. J.K. Johnson, J.A. Zollweg, K.E. Gubbins, Mol. Phys. 78 (1993) 591–618. V.I. Harismiadis, N.K. Koutras, D.P. Tassios, A.Z. Panagiotopoulos, Fluid Phase Equilib. 65 (1991) 1–18. A. Lofti, J. Vrabec, J. Fischer, Mol. Phys. 112 (1995) 173–197. F.A. Escobedo, J. Phys. Chem. 113 (2000) 8444–8456.