Article pubs.acs.org/JCTC
Leveraging Gibbs Ensemble Molecular Dynamics and Hybrid Monte Carlo/Molecular Dynamics for Efficient Study of Phase Equilibria Thomas E. Gartner, III,† Thomas H. Epps, III,*,†,‡ and Arthi Jayaraman*,†,‡ †
Department of Chemical and Biomolecular Engineering, University of Delaware, 150 Academy Street, Newark Delaware 19716, United States ‡ Department of Materials Science and Engineering, University of Delaware, 201 DuPont Hall, Newark Delaware 19716 S Supporting Information *
ABSTRACT: We describe an extension of the Gibbs ensemble molecular dynamics (GEMD) method for studying phase equilibria. Our modifications to GEMD allow for direct control over particle transfer between phases and improve the method’s numerical stability. Additionally, we found that the modified GEMD approach had advantages in computational efficiency in comparison to a hybrid Monte Carlo (MC)/MD Gibbs ensemble scheme in the context of the single component Lennard-Jones fluid. We note that this increase in computational efficiency does not compromise the close agreement of phase equilibrium results between the two methods. However, numerical instabilities in the GEMD scheme hamper GEMD’s use near the critical point. We propose that the computationally efficient GEMD simulations can be used to map out the majority of the phase window, with hybrid MC/MD used as a follow up for conditions under which GEMD may be unstable (e.g., nearcritical behavior). In this manner, we can capitalize on the contrasting strengths of these two methods to enable the efficient study of phase equilibria for systems that present challenges for a purely stochastic GEMC method, such as dense or low temperature systems, and/or those with complex molecular topologies.
I. INTRODUCTION Fluid phase equilibria is of critical importance to a vast array of scientific and engineering applications, and interest in phase behavior has motivated the development of a multitude of computational methods to predict and understand phase coexistence.1−5 One such widely applied computational method is Gibbs ensemble Monte Carlo (GEMC),6,7 which has been used to study the phase behavior of small molecules,7−11 short chain hydrocarbons,12−14 polymers,15,16 and fluid mixtures,7,17 and has been adapted to membrane equilibria,7 fluids confined in porous media,18 and many other systems.19 In the GEMC method, the two phases in equilibrium are simulated in two physically separate simulation boxes, without modeling the interface between them explicitly. This strategy allows for the determination of coexistence properties with reduced system size effects and reduced simulation time in comparison to the large systems and long simulations necessary for the inclusion of both phases in the same simulation domain.20 Despite these advantages, GEMC suffers from prohibitively low computational efficiency in noteworthy cases (e.g., dense and macromolecular fluids) because of its dependence on stochastic particle insertion when transferring molecules between phases. Numerous strategies have been employed to improve the applicability of GEMC to dense and macromolecular systems, for which the probability of successful random particle insertion is vanishingly small. For example, configurational bias techniques facilitate higher efficiency insertions for chains by sequentially inserting each bead of a chain rather than the entire chain at once.12−14,16,21 Expanded ensembles further improve © 2016 American Chemical Society
efficiency by coupling the sequential insertion of beads in a chain with the sequential deletion of beads from a chain in the other phase.15 A similar technique uses continual fractional component MC moves to gradually insert and delete whole molecules.22 Alternatively, cavity biased insertion techniques improve efficiency by searching for voids in the insertion domain rather than choosing a location at random.23 These advanced methods help to address some of the limitations of GEMC; however, they introduce additional layers of complexity to any implementation of a Gibbs ensemble simulation scheme. Another approach to increasing the applicability of the Gibbs ensemble is to translate the Gibbs ensemble idea to a molecular dynamics (MD) formalism.24−27 Gibbs ensemble molecular dynamics (GEMD) works by defining a fourth degree of freedom for each bead/particle in the system, ξ, which governs whether a bead/particle is treated as existing in one phase (e.g., liquid phase) or the other (e.g., gas phase). Components transfer dynamically between the phases through motion along this ξ coordinate, analogous to motion along the three positional degrees of freedom in the normal Cartesian coordinates. This GEMD method allows for the simultaneous transfer of multiple beads/particles and, unlike GEMC, does not depend on random particle insertion moves. Additionally, GEMD harnesses some of the usual advantages of MD simulation, including its ability to equilibrate complex systems Received: June 4, 2016 Published: September 29, 2016 5501
DOI: 10.1021/acs.jctc.6b00575 J. Chem. Theory Comput. 2016, 12, 5501−5510
Article
Journal of Chemical Theory and Computation without the advanced stochastic moves needed in GEMC, and to calculate dynamic properties such as diffusion coefficients, which is not possible with the GEMC method. Hentschke and co-workers presented initial GEMD results for the LennardJones (LJ) fluid24 and expanded their work to n-hexane,25 adsorption in a solid pore,26 and the solvent swelling of a polymeric network28 in later publications. However, these past implementations of GEMD suffer from numerical instabilities in particle transfer (described in more detail in the Results and Discussion section) that need to be addressed to make the method more useful for phase equilibria studies. Herein, we present an extension and modification of the GEMD method to improve numerical stability and to control the particle transfer along the ξ coordinate. We applied a velocity rescaling thermostat to the ξ coordinate to control the exploration of the ξ domain and demonstrated the use of a “ξ annealing” procedure to remove dependence on the initial configuration of the simulation and achieve a stable equilibrium condition. The updated GEMD method for the single component LJ fluid achieved equivalent phase equilibrium results in comparison to those from a hybrid MC/MD method in which MD displacement runs are interspersed with Monte Carlo volume and particle exchanges. These two methods also agreed with prior simulation results by Panagiotopoulos (GEMC)7 and Hentschke (GEMD)26 for the same system. Notably, GEMD required a small fraction of the simulation time needed for hybrid MC/MD to achieve the same phase equilibrium results, and GEMD displayed an increasing advantage in computational efficiency at low temperatures. On the basis of these results, we conclude that the GEMD method represents an efficient simulation technique to study phase behavior of systems challenging for stochastic GEMC without resorting to advanced MC biasing moves.
Figure 1. Illustration of the GEMD scheme. Red represents vapor phase beads (near ξ = 0); blue represents liquid phase beads (near ξ = 1); and white represents transition state beads. The top snapshot contains all the components in the system (both liquid and gas phase) that share the same Cartesian space. The lower snapshots show the system separated along the ξ coordinate.
U=
i
