Continuum configurational bias Monte-Carlo studies of alkanes and polyethylene

Fluid Phase Equilibria, 83 (1993) 323-331 Elsevier Science Publishers B.V., Amsterdam

323

CONTINUUM CONFIGURATIONAL BIAS MONTE-CARLO STUDIES OF ALKANES AND POLYETHYLENE Juan J. de Pablo, Manuel Laso, Ulrich W. Suter Institut fiir Polymere, Eidgeniissische Technische Hochschule. ETH-Zentrum, CH-8092 Ziirich, Switzerland. Henry D. Cochran Chemical Technology Division, Oak Ridge National Laboratory. Oak Ridge, Tennessee 37831, USA.

ABSTRACT A novel class of Monte-Carlo simulation methods has recently been developed for the study of dense polymeric systems. Such methods allow the calculation of chemical potential and solubilities in macromolecular systems. In this work, results are presented for solubilities of long alkanes in polyethylene and for dilute solutions of long alkanes in near-critical solvents.

INTRODUCTION Until recently, simulations of dense polymeric systems have been hampered by the difficulties of sampling phase-space efficiently. On the one hand, moleculardynamics methods can only sample time scales that are generally short compared with most relaxation times in dense polymeric systems (Kremer and Grest, 1990). On the other hand, Monte-Carlo studies of polymers have generally been carried out on a lattice (Binder, 1987). Monte-Carlo simulations of dense polymers in a continuum have been scarce due to the difficulties associated with sampling configuration space adequately. We have recently developed a novel Monte-Carlo method that appears to overcome many of the limitations of previously available Monte-Carlo techniques for polymer simulations (de Pablo et al. 1992a). We have since applied the so-called Continuum-Configurational-Bias (or CCB) concept to a variety of realistic systems with promising results (de Pablo et al., 1992 a,b,c, Laso et al., 1992). In parallel work, Frenkel et al. (1992) have proposed alternative Monte-Carlo schemes that are similar to the CCB method and have applied them to investigate several aspects of

324

chains in solution. In this paper we discuss some of the systems that we’ve studied so far and we present the results of such studies.

CONTINUUM

CONFIGURATIONAL

BIAS (CCB) METHODS

In their simplest version, Monte-Carlo simulations of molecular systems are based on generating uniform random displacements of the molecules’ coordinates and accepting or rejecting these according to importance sampling criteria (Allen and Tildesley; 1991). For a Lennard-Jones fluid at constant density and temperature, for example, a particle is chosen at random and, after generating a random move of that particle, the new configuration is accepted or rejected using Metropolis sampling (Allen and Tildesley, 1987). ‘Ihe presence of additional degrees of freedom makes the simulation of polymer molecules more complicated. In addition to random displacements of the center-ofmass coordinates, the internal configuration of the molecules must also be changed. A small random change of a torsional angle, however, can give rise to large displacements of the molecules’ atoms, thereby generating steric overlaps and unacceptable trial configurations. To overcome this problem, trial configurations can be generated by introducing a bias that avoids overlaps and favours low-energy configurations, provided that such bias is subsequently removed from the simulation. In the CCB method, a chain molecule, chosen at random, is cut at a random point. The molecule is then regrown, segment by segment, to its original length. Every time a segment is appended to the growing chain, the torsional angle that dictates the position of the atom (or site) at the end of the segment is chosen in such a way as to favour low-energy configurations over high energy (or overlapping) configurations. Trial configurations are then accepted or rejected according to modified sampling criteria (de Pablo et al., 1992a) that take into account and exactly remove the bias introduced by our growing scheme. The details concerning this biased sampling of configuration space have appeared in several publications (de Pablo et al. 1992a,b,c, Laso et al., 1992, Frenkel et al., 1992) and need not be given here. Instead, in the following sections we describe some of the systems that we have studied using the CCB method.

HENRY’S

CONSTANTS

OF ALKANES

IN POLYETHYLENE

In this work, polyethylene is described according to a unified atom representation in which methylene groups are replaced by Lennard-Jones pseudoatoms. Interaction sites are connected by rigid bonds of 1.53 8, at rigid angles of 112’. Sites in different molecules or sites in the same molecule but separated by more than three bonds interact through a Lennard-Jones potential energy function with parameters o=3.94 8, and &=0.098 kcal/mol. In addition to Lennard-Jones

325

interactions, torsions around carbon-carbon bonds are hindered potential-energy function (Baschnagel, 1990).

by a torsional

Successive configurations for polyethylene melts are generated in an isothermalisobaric (Np7’) statistical-mechanical ensemble. The melts studied here consist of ten chains of 71 methylene groups each. These configurations are stored for subsequent studies, such as the simulation of Henry’s constants for intermediate molecular weight solutes in polyethylene. The CCB method can be used in conjunction with Widom’s test or ghost particle method to determine the chemical potential of a dense system of chain molecules or the chemical potential of an articulated solute in a dense system (Siepmann, 1990, de Pablo et al. 1992b, Frenkel et al., 1992). Since the chemical potential of a solute at infinite dilution is related to its Henry’s constant, the solubility can be estimated for low concentrations (or low pressures). Figure 1 shows weight fraction Henry’s constants for alkanes dissolved in polyethylene at 1 bar and at 513 K and 423 K. The results of our simulations compare favorably with experiment. Note that both the correct temperature and chain-length dependence are predicted by our calculations. Henry’s

law

020,’ o0

Dnh

*

Cal

01 Milloncy

k l+ilU4”ll,

Bh0YL nlolecLIIc.

(,11-/C)

MPT rnscm,,,r

r g

o-

+

t

2 c I

/ T = 423

4-

4 I<

T = 513 K

? 01

t

0

t

CCB Gnbbs-ensemble

2 2

1 Number

(1 or Cnrbon

0

10

30

40

sim. 4 50

hhrnrr P

(bar)

Fig. 1: Henry’s constants for alkanes in polyethylene. Fig. 2: Solubility of pentane in polyethylene as a function of pressure. SOLUBILITY OF ALKANES IN POLYETHYLENE While the solubility of a gas at low pressure in a polymer melt can be determined accurately according to Henry’s law, such a law is only valid for low concentrations of the solute in the polymer. We have therefore developed a CCB scheme to simulate phase equilibria for chain molecules in the Gibbs ensemble (Panagiotopoulos et al., 1988, Laso et al. 1992). Incorporation of a CCB sampling scheme into Gibbs-ensemble simulations (CCB-Gibbs) extends the usefulness of such calculations to the study of moderately long chain molecules. While such studies are

326

important in their own right, it is important to recognize that the CCB Gibbs method is not suitable for simulation of phase equilib~a for long polymer molecules. The real usefulness of the CCB Gibbs method, however, becomes apparent in simulations of the solubility of moderately Iong solute molecules (i.e. plasticizers) in dense polymer matrices (Laso et al., 1992). In such a method, the simulation is carried out in two boxes; one box contains the mixture of polymer and solute molecules; a second coexisting box contains an essentially pure gas or liquid of solute molecules. Such simulations can be carried out at different pressures, thereby generating a solubility curve for the solute in the polymer. Figure_ 2 shows the solubiIity of pentane in polye~ylene as a function of pressure. The slope of this curve at the origin corresponds to the inverse of the Henry’s constant for this system. Note that even at moderate pressures (10 bar), deviations form Henry’s law are significant; the CCB Gibbs methodology, however, provides a direct route towards estimating the correct solubility behavior of articulated solute molecules (e.g. pentane) in dense polymers. POLYMERS IN SUPERCRITICAL

SOLUTION

The potential usefulness of supercritical fluid extraction (SFE) processes has prompted a number of theoretical studies of the behavior of solutes in supercritical solutions (Debenedetti, 1987, Kumar et al., 1987, de Pablo, 1990, Wu et al., 1990). While such studies have significantly increased our understanding of near-critical mixtures, they have been restricted to spherical Lennard-Jones particles dissolved in near-critical Lennard-Jones solvents, However, it is important to recognize that SFE processes are generally aimed at separating or purifying mixtures of relatively large substances such as heavy oils (for enhanced oil recovery), polymers (for polymer fmctionation) or biomolecules (e.g. in supercritical fluid chromatogmphy~. It has been established that near-critical mixtures can exhibit IocaI sokent density enhancement (which we shall call clustering for short) or, alternatively, local solvent density depletion around a solute. Since we expect such formations to have an effect on the configurations that long, articulated solute molecules can assume, we have initiated a study aimed at elucidating the behavior of such long molecules in near-critical systems. The system that we discuss here consists of an alkane molecule of 24 segments at infinite dilution in near-critical methane and in near-critical cyclohexane (we use 2200 solvent particles). I At near-critical conditions, the chemical details of the molecules do not have significant qualitative consequences on the properties of the system; we have therefore chosen to represent both methane and cyclohexane by spherical Lennard-Jones particles.

327

Figure 3 shows the simulated solvent-solvent and the chain site-solvent radial distribution functions (rdf) for the system at T*=1.03 and p*=l.OO. Superscript * denotes a quantity that has been reduced with respect to the critical properties of the solvent. The reduced critical constants for the Lennard-Jones fluid that we use here are T,=T,k,/&=1.31 and p,=p,03=0.31 (Panagiotopoulos et al., 1988). To ensure a good resolution of solute-solvent. correlations, these data were collected over simulations of 100 x lo6 steps. Further, to improve the statistics for the conformation of the dissolved chain we use a sampling scheme that preferentially moves those molecules. located in a spherical region defined around the solute (Allen and Tildesley, 1987). The solid line shows the pure methane rdf calculated from the Omstein-amike equation with a Percus-Yevick closure (such an equation has been shown to yield remarkably accurate results for LJ spheres at conditions similar to those investigated here [de Pablo, 19901). The good agreement between theory and simulation indicates that the system of one chain in 2200 methane particles does correspond to infinite dilution. Further, simulations with 5000 solvent particles give results that are identical to those with 2200 solvent particles. The site-solvent rdf shown in Figure 3 shows that clusters of solvent molecules are formed around the chain. These clusters provide a good solvent environment which causes the chain to assume an extended configuration. Figure 4 depicts sitesolvent rdf’s, also at T*=1.03 but at several densities. The strongest tendency to form clusters is observed at densities slightly below the critical (p*=O.86), in good agreement with earlier observations (Debenedetti, 1987). At liquid-like densities, the clusters disappear. At low, vapor densities (p*=O.71), we observe the formation of entities composed of a chain surrounded by a few solvent particles; the correlation between the chain and the solvent particles decays relatively fast (after only 4 solvent diameters) and, as shown below, the partial molar volume of the chain is positive.

0

10

20 r 6,

30

0

10

20

30

r (4

Fig. 3: Simulated solvent-solvent and chain site-solvent rdf’s at T*=1.03 and p*=l.OO. Fig. 4: Chain site-solvent rdf’s at T*=1.03 and several densities. As solvent density and quality change, the chain appears to undergo configurational changes. Table I gives the mean square end-to-end distance
328 the mean square radius of gyration for the chains in methane at several densities. These data indicate that the best solvent for the chains (the largest
Figure.5 shows site-solvent rdf’s for the chain in metha2e and for fh,e chain in cyclohexane. The conditions for both systems are the same (T =1.03 and p =l.). We must point out, however, that for the systems studied here (2200 cyclohexane particles),’ the rdf for pure cyclohexane is not identical to that obtained when a chain is present. Effects of the size of the system simulated are still important, and in this regard we cannot say that the simulated C24-cyclohexane system is at infinite dilution. Furthermore, it must be remembered that the absolute (rather than reduced) temperature of the cyclohexane simulation is 570 K, compared with 115K for the methane simulations. It is nevertheless instructive to discuss briefly the behavior of this system. Note that in near-critical cyclohexane clustering is not as strong as in methane. Further, the chain in cyclohexane assumes the configuration that it would have in vacuum (at the same temperature). While methane provides a very good solvent for a SCFE based process, cyclohexane would provide a relatively poor solvent.

Fig. 5: Chain site-solvent rdf’s in methane and cyclohexane at T*=1,.03 and p*=l.OO. The effects of clustering in dilute mixtures can now be expressed in terms of thermodynamic, measurable quantities by means of a statistical mechanical theory. Here we use the theory of solutions of Kirkwood and Buff (1951), and we assume that the chain can be regarded as a collection of identical atoms that are fused together; the connectivity and structure of the chain are implicit in the site-solvent rdf, which is obtained directly from a molecular simulation. At infinite dilution the partial molar volume of the chain, Tsw , is determined according to

329

~5~” =

pnmiY,- = 1 + p (h,(O)

- h,,(O) 1

where &(O) denotes the Fourier transfomt of the net correlation function [h(r)=gfr)-11 between species 1 and 2 evaluated a zero wave vector. Subscript 1 indicates a site of the chain, subscript 2 corresponds to the solvent, and n, is the number of sites of the chain. For infinitely dilute mixtures, the isothermal compressibility (or) of the system is that of the pure solvent, and is given by

P

-% P

= 1 + p QO)

Together, the infinite-dilution partial molar volume and the isothermal compressibility can be used to define a cluster size r according to

Here r denotes the number of solvent molecules that are found around the chain in excess to those found for a unifom~ distribution of the solvent at the same density. Table II gives our estimates for compressibility, partial molar volume and cluster size for the systems studied here. The partial molar volume of the chain in near-critical methane is very large and negative. It assumes a minimum at p*=O.86; for this density, the size of a cluster is very large, in accordance with our earlier comments regarding the strength of critical effects at densities slightly below the critical. The so-called enhancement factor expresses the correction factor which must be applied to a simple (ideal-gas) expression to calculate the solubility of a heavy component in a gas. This factor is inversely proportional to the fugacity coefficient of the solute in the gas; since its logarithm can be written as an integral of the partial molar volume of the solute from zero to the pressure of the system, a large negative partial molar volume translates into an exceedingly small fugacity coefficient and a consequently high enhancement factor. With the methods discussed here, solute partial molar volumes can be predicted from knowledge of molecular interactions only; these methods could therefore provide the basis for a useful tool for selecting adequate solvents for SCFE processes.

CONCLUSIONS We have developed a novel Monte-Carlo technique that can sample efficiently the configuration space of large, articulated molecules. This method has allowed US to simulate, for the first time, chemical potentials and phase equilibria for

330

macromolecular systems. The method has useful applications to situations that are often encountered in engineering. In this work we have addressed the calculation of PVT properties for long alkanes, the solubility of moderately long chain molecules in polymers, and the selection of suitable solvents for supercritical fluid extraction. In these initial studies, however, we have focused our attention on systems with relatively simple inte~olecular interactions; their challenge and their interest resides in the variety of ~onfo~ations that the molecules can assume. It now remains to be seen if methods such as those described here will also be useful in describing systems of more compIex chains that interact through more complicated potential-energy functions.

ACKNOWLEDGMENTS Generous allocations of CPU time on the NEC SX-3 of the Centro Svizzero di Cafcolo Scientific0 and on the Cray Y-MP and IBM 3090 of the ETH-Ztirich are gratefully acknowledged. We are also grateful for partial support of this project from Grant No. 03M4028 of the Bundesministerium fur Forschung und TechnoIogie of the Federal Republic of Germany and to Bayer AG. HDC acknowledges support from the Division of Chemical Sciences of the U.S. Department of Energy under Contract NO. DE-AC05840R21400 with Martin Marietta Energy Systems, Inc.

Table I - Mean square end-to-end distance and radius of gyration for C,, dissolved in near-critical methane and cyciohexane

331 Table II - Isothermal compressibility, partial molar volume of the chains, and cluster size for infinitely dilute near-critical solutions of C24

a) accuracy: pKT/p ca. AO.1, pvS ca. +lOO., r ca. *loo.

REFERENCES Allen, M.P., Tildesley, D.J., “Computer Simulation of Liquids”, Clarendon Press, Oxford (1987). Baschnagel, J., Diploma&it, Johannes-Gutenberg-Universitat, Mainz, (1990). Binder, K. (ed.), “Applications of the Monte Carlo Method in Statistical Physics”, Springer, (1987). de Pablo, J.J., PhD Thesis, University of California at Berkeley (1990). de Pablo, J.J., Laso, M., Suter, U.W., J. Chem. Phys., 96, 2395 (1992a). de Pablo, J.J., Laso, M., Suter, U.W., J. Chem. Phys., 96, 6157 (1992b). de Pablo, J.J., Laso, M., Suter, U.W., Siepmann, J.I., submitted for publication (1992c). Laso, M., de Pablo, J.J., Suter, U.W., J. Chem. Phys., August 1992, in press. Debenedetti, P.G., Chem. Eng. Sci., 42, 2203 (1987). Dee, G.T., Ougizawa, T., Walsh, D.J., preprint (1991). Frenkel, D., Mooij, G.C.A.M., Smit, B., J.Phys.: Cond.Matt., 4, 3053, (1992). Kremer, K., Grest, G.S., J. Chem. Phys., 92, 5057 (1990). Kumar, S.K., Chabria, S.P., Reid, R.C., Suter, U.W., Macromolecules, 20, 2550 (1987). Kirkwood, J.G., Buff, F-P., J. Chem. Phys., 19, 774 (1951). Maloney, D.P., Prausnitz? J.P., AIChE Journal, 22, 74 (1976). Panagiotopoulos, A.Z.;Quirke, N., Stapleton, M., Tildesley, D.J., Mol.Phys., 63,527, (1988). Wu, R.-S., Lee, L.L., Cochran, H.D., Ind. Eng. Chem. Res. 29(6), 977 (1990).