2274
J. Phys. Chem. B 2007, 111, 2274-2287
Coarse-Grained Molecular Dynamics Modeling of Associating Fluids: Thermodynamics, Liquid Structure, and Dynamics in the Limit of Zero Association Strength Ting Li*,†,‡ and Erik Nies*,†,§ Polymer Research DiVision, Department of Chemistry, Katholieke UniVersiteit LeuVen, Celestijnenlaan 200F, B-3001 HeVerlee, Belgium, State Key Laboratory of Polymer Physics and Chemistry, Center for Molecular Science, Institute of Chemistry, Chinese Academy of Sciences, Zhongguancun, Beijing 100080, People’s Republic of China, and Laboratory of Polymer Technology, EindhoVen UniVersity of Technology, P.O. Box 513, 5600MB EindhoVen, The Netherlands ReceiVed: September 20, 2006; In Final Form: NoVember 28, 2006
A continuous coarse-grained potential model for associating fluids, consisting of an off-center specific site bonded with a harmonic potential to a center particle, has been developed and used in canonical molecular dynamics simulations. The thermodynamic, structural, and dynamic properties of the limiting nonassociating reference coarse-grained fluid are investigated as functions of the mass distribution and bond strength between center and site particles. It is theoretically shown and confirmed by simulation that in this limit variations in these potential parameters do not alter the equation of state of the reference coarse-grained fluid but that they have profound influences on both the translational and the rotational dynamics. From the simulation results we arrive at some guidelines that should be kept in mind in the selection of appropriate values for the model parameters. This work provides the precursory knowledge for the study of coarse-grained associating fluids using the conventional molecular dynamics method.
1. Introduction Molecular association, such as that introduced by hydrogen bonding, plays an extremely important role in the properties of complex fluids. From pure small molecule systems to macromolecular mixtures, it has attracted considerable research interests of both experiment and theory in the past decades. For instance, a distinct liquid-liquid closed-loop miscibility gap in aqueous fluid mixtures can be observed,1-5 and the strong, localized, temperature-sensitive, and directional associating interaction accounts for the miscibility at lower temperatures and the existence of a lower critical solution temperature (LCST).6 The statistical associating fluid theory (SAFT)7-9 has been proven a successful theoretical method to deal with associating fluids as demonstrated by numerous comparisons with experimental and computer simulation results. On the basis of Wertheim’s first-order thermodynamic perturbation theory (TPT),10-14 SAFT incorporates the associating intermolecular interaction as a perturbation to a reference system accounting for the dispersive interactions. An equation of state (EOS) was developed and was applied to a wide range of associating fluids including water, small molecule organic compounds, and polymers. Nowadays, several SAFT-based theories are available for the study of the thermodynamical properties of various materials involving association interactions.15,16 As a molecular-based EOS, SAFT was extensively tested and improved using molecular simulations. Even in the first publica* Authors to whom correspondence should be addressed. Phone: +32 16 327614 (T.L.); +32 16 327481 or +32 16 327418 (E.N.). Fax: +32 16 327990 (T.L.); +32 16 327990 (E.N.). E-mail:
[email protected];
[email protected]. † Katholieke Universiteit Leuven. ‡ Institute of Chemistry, Chinese Academy of Sciences. § Eindhoven University of Technology.
tion that served as the precursor to what now is known as the SAFT theory, Chapman et al. had already used the Metropolis Monte Carlo simulation method to test their theory for binary associating liquid mixtures.17 The pair potential model used in their simulation and theory for the molecular interactions consisted of a hard sphere and an off-center point charge dipole. Excellent agreement between theory and simulation was found, including the consistency for the energy of mixing and the fraction of monomers, i.e., the molecules not involved in an association interaction. Later on, similar models were extensively used in the studies of associating fluids in both theory and simulation. A common character of these kinds of potential models is that the nonassociating interactions are modeled by, for example, a central hard sphere or a square-well potential for the dispersive interaction, and the specific interactions are accounted for by off-center attractive sites. In almost all of the previous studies, except for fully atomistic simulations, the specific interaction was also defined to be a square-well potential, adding some constraints on distance and angle between the interacting sites to ensure the saturation of the specific interaction. These potential models are applied in Monte Carlo simulations and in the SAFT theory without any difficulty, and good agreement between theory and simulations was observed. However, the discontinuous square-well potential for the specific interaction is not suitable to be applied in conventional molecular dynamics (MD) simulations, because the potential has no continuous first-order spatial derivative, making the forces required in MD simulations indeterminate. This problem hindered the study of transport properties of associating fluids, which certainly is also an important and practical aspect of complex fluids. Therefore, a special molecular dynamics algorithm is needed, and in fact, the original MD algorithm developed by Alder and Wainwright was such a method that could deal with discontinuous potentials.18 Later, this type of
10.1021/jp066162k CCC: $37.00 © 2007 American Chemical Society Published on Web 02/09/2007
Coarse-Grained MD Modeling of Fluids MD algorithm was further developed in the discontinuous molecular dynamics (DMD) simulation method. The essential idea of DMD is that the particles in the system retain linear motion until a collision happens, and collisions only happen at, for example, the repulsive and attractive boundary of the squarewell potential. This simplification endows DMD with great efficiency, because the frequent position or velocity updates, at the core of conventional MD, are not necessary and a large time step can be used. DMD was successfully employed in some previous thermodynamic studies. To make DMD applicable to polymers, Rapaport introduced covalent bond fluctuations between two adjacent beads, allowing them to vibrate between σ and σ(1 + δ), and thus bond-stretch collisions were introduced.19 Bellemans then further modified the bond fluctuation contribution by allowing bonded beads to slightly interpenetrate.20 Later, the method was extended and applied by Chapela et al.21 and by Denlinger and Hall22 to hard chain systems. However, before Gulati and Hall’s work,23 DMD could only be conducted in the microcanonical ensemble. On the basis of the conventional DMD, Gulati and Hall developed a new discontinuous canonical molecular dynamics (DCMD) algorithm by introducing a new type of “ghost” collision, which was in fact an implementation of Andersen’s thermostat.24 From then on, DCMD was effectively used to study the thermodynamics of homopolymers and clusters,25 the phase behavior of copolymers,26,27 and the phase behavior of polymer/small molecule mixtures.28 Due to the similarity of the discontinuous potential models used in DMD and in SAFT, it is straightforward to utilize DMD to simulate associating fluids. Liu and co-workers investigated one bond site and four bond sites models in their DMD simulation as a model for water, using the same models in the SAFT theory.29 They obtained good agreement between the thermodynamic averages from the DMD simulations and the SAFT predictions. They also studied the effects of the centerto-site bond length and the size and mass of the specific site on the dynamics of the system. The velocity autocorrelation and the self-diffusion coefficient were obtained, and it was shown that the presence of specific interactions decreases the selfdiffusion coefficient compared to that of the hard sphere system, especially at low densities. This study showed that DMD is an alternative simulation method to study the associating fluids and is capable of providing dynamic properties. Although the DMD method makes it possible to study the dynamics there are also some drawbacks of the method. One problem is directly related to the foundation of the DMD algorithm. In DMD, particles do not experience continuous forces but only elastic collisions. If we believe particles follow nonlinear Newtonian motion, then their trajectories in DMD are far from reality. This may not be a problem for the thermodynamic properties but will certainly affect the dynamics. Conservatively speaking, very likely the short-time dynamical properties will be different. And further careful investigation of this issue is required. Another problem is what ensembles can be implemented in the DMD algorithm. So far, for the canonical ensemble only Andersen’s constant temperature algorithm has been implemented, which is a timeirreversible thermostat due to its stochastic characteristics.30 Moreover the DMD algorithm is not yet available for, for example, the NPT and the grand-canonical ensembles. The aforementioned issues are not present in the conventional molecular dynamics method. Surprisingly, as far as we know, the conventional MD method has not been applied in the study of coarse-grained associating fluids. Although conventional MD simulations are abundantly applied at the atomistic level for hydrogen-bonding systems,
J. Phys. Chem. B, Vol. 111, No. 9, 2007 2275 the fully atomistic model is not suited for the direct comparison with the model used in the SAFT theory and many coarsegrained MC simulations of associating fluids. Therefore, in this work we develop a continuous coarse-grained potential model for the study of associating fluids that is applicable in conventional continuous MD or off-lattice MC simulations. Coarse-grained interaction potentials are important to accelerate the simulations of complex fluids.31-33 Coarse-grained models can be studied in simulations on larger time and/or length scales than fully atomistic models, and in a recent review many possible coarse-graining schemes are discussed.34 A key point in the coarse-grained simulations is the selection of appropriate functional forms for the coarse-grained interaction potentials and the choice of correct values for the coarse-grained model parameters. It is, for instance, important that the coarse-grained model correctly reproduces the properties of the real system at higher length and time scales, and improper choices may result in unwanted behavior.35 Another motivation is related to the development of the SAFT theory. The original SAFT and many subsequent adaptations and refinements take the hard sphere system as a reference. Recently, Blas and Vega have proposed a modified version of the SAFT theory (soft SAFT), based on the extension of Wertheim’s theory to Lennard-Jones chains by Johnson et al. The repulsive and dispersive free energy terms are combined by taking the Lennard Jones fluid as the reference system. It is anticipated that this modification will improve the accuracy of the theoretical prediction due to more realistic dispersion interactions. Obviously, to verify such a theory, a continuous potential model is needed in the simulation. Furthermore, comparing the vast amount of thermodynamics studies of associating fluids, dynamic studies are still rare, a direct consequence of the lack of an adaptive continuous MD simulation model. However, it is worth mentioning that continuous site-specific saturation potentials have been used in MD simulations in the study of the complex self-assembly of, for example, capsids as in Rapaport36 and Hagan.37 Last but not least, we must emphasize the importance of collective or cooperative motion in associating fluids. Multiple strong and directional specific interactions can create relatively stable networks. The relaxation and diffusion behavior is dependent on the dynamics of the destruction and reformation of specific bonding interactions. The reorganization of the network needs cooperative motion, as observed in liquid water.38,39 In the conventional MC method it is not easy to provide this kind of collective configurational changes, whereas the MD simulation technique is the method of choice when cooperative updates are needed. In the course of the study of associating fluids with the proposed coarse-grained interaction potential, we observed that the dynamical behavior of the underlying ideal particles (i.e., the particles with the specific interaction strength in the coarsegrained interaction potential switched off) is quite complicated, and therefore, the dynamic behavior of the ideal system should be investigated carefully as it serves as a reference system for the systems with the specific interaction strengths switched on. Therefore, in this work we report on a detailed study of the thermodynamic, structural, and dynamical properties of fluids consisting of particles interacting with the interaction potential developed in this work but with the specific interaction strength equal to zero. The paper is organized as follows. First, the potential model and simulation details are described. Then, the simulation results of the thermodynamical, structural, and dynamical properties are presented and discussed whereby the
2276 J. Phys. Chem. B, Vol. 111, No. 9, 2007
Li and Nies
effects of the different simulation parameters are addressed. Following this, some guidelines are provided that should be kept in mind when selecting appropriate parameters for a coarsegrained simulation. Finally, conclusions are made. 2. Potential Model and Theory 2.1. Potential Model. The specific interaction in associating fluids is short-ranged and highly oriented and in many cases also acts as a saturation interaction. The latter property implies that once the specific interaction is formed it is not available for another specific interaction with another molecule. A simple molecular model for saturation interactions commonly used in literature that finds its inspiration in the successful Wertheim theory for saturation interactions10-13 consists of a spherical particle “labeled” with an off-center site that accounts for the saturation interaction. The structural and thermodynamic properties of these simple models have been predicted using the Wertheim theory and are found to be in good agreement with molecular simulation results. In this paper we discuss molecular dynamics simulations of molecules that consist of a simple spherical center particle c with an (off-center) short-ranged saturation interaction site. In this work we call such a molecule a “Tmer”. A schematic representation of two Tmer molecules is depicted in Figure 1. Two molecules interact with a short-ranged specific interaction potential Φijss between sites si and sj as well as with a longerranged interaction potential Φijcc between the two centers. The bond vector lics ) ris - ric between the off-center site and the center particle defines the orientation of the molecule. The potential energy of a collection of N such molecules is assumed to be pairwise additive and is given by N
Utotal ≡ UvdW + Uspec + Ubond ) N
N
Φijcc(rijcc) + ∑ ∑ i)1 i
