Molecular dynamics study on solvent clustering in supercritical fluid solutions based on particle radial kinetic energy

[OlmlglllA ELSEVIER

Fluid Phase Equilibria 104 (1995) 317-329

MOLECULAR DYNAMICS STUDY ON SOLVENT CLUSTERING IN SUPERCRITICAL FLUID SOLUTIONS BASED ON PARTICLE RADIAL KINETIC ENERGY Chee Chin Liew, Hiroshi lnomata and Shozaburo Saito

Research Center of Supercritical Fluid Technology, Department of Molecular Chemistry and Engineering, Tohoku University, Sendai, JAPAN 980 Keywords: theory, computer simulation, solvation dynamics, supercritical fluids, Lennard-Jones. Received 16 May1994; accepted in final form 19 September 1994

ABSTRACT

The objective of this work was to study the dynamic structures of solute-solvent clusters in supercritical fluids with molecular dynamics (MD) simulation. A new definition tbr cluster formation, based on the relative radial kinetic energy, is proposed. The method defines a cluster member as a solvent molecule which has a relative solute-solvent radial kinetic energy less than or equal to their pair potential energy. MD Simulations of Lennard-Jones (L J) mixtures are shown for CO 2 + pyrene and CO 2 + model substances using one solute molecule and 863 solvent (CO2) molecules at T*=1.40, p*=0.35 via a Nose-Hoover thermostat. The average number of member CO 2 molecules in a sphere with diameter of 13.5Oco 2 surrounding the pyrene solute was found to be 58. By tagging cluster and non-cluster members, we observed that approximately half of the particles in the first cluster shell exchanged in about 4 ps. Our simulations show that the number of remaining initial cluster members decay very rapidly in the beginning and approach a constant exponential decay after around 1 ps in all systems. The behavior of the time dependent radial distribution function of the remaining initial members indicated that only those members which were initially located in the solvation shell remained as members after about 1 ps. We considered the number of cluster members on the constant decay curve as the number of solvation members. By using the exponential decay rate, the cluster half-life and the solvation number of the CO 2 + pyrene system was found to be around I ps and 16 solvent molecules, respectively. Other LJ mixtures were also studied and it was found that the LJ energy parameter strongly influences the decay rate while the LJ size parameter directly affects the solvation number. Molecular weight only had a small effect on solvation number.

0378-3812/95/$09.50 © 1995 - Elsevier Science B.V. All rights reserved SSDI 0378-3812 (94) 02657-2

318

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

INTRODUCTION Elucidation of the clustering phenomenon in supercritical fluid solutions is regarded as being the key to understanding of the dramatic solubility enhancement observed of chemical substances in supercritical fluids (Eckert et al., 1983; Johnston and Penninger, 1989; Bruno and Ely, 1991). Furthermore, since the cluster formation must affect transport properties such as diffusivities and mass transfer coefficients, determination of methods to simulate the clustering is also thought to be important for effectively using supercritical fluids as a reaction medium (Bright and McNally, 1992; Eckert and Knutson, 1992). Due to different experimental methods or theoretical viewpoints, dozens of different definitions of clusters in supercritical fluids have been proposed (Eckert et al., 1983; Kim and Johnston, 1992; Petsche and Debenedetti, 1989; Cochran and Lee, 1989; Wu et al., 1989). As a general concept of this paper, a cluster is referred as an aggregate of particles which (1) have neither solid structures nor chemical binding among each other, (2) are loosely bounded together by van der Waals or other attractive forces and (3) have a finite size. The clustering phenomenon consists of both long-ranged and short-ranged forces. These are discussed in our definition below. A statistical mechanic theory of cluster formation was proposed by Hill (1955, 1956), for description of imperfect gases. His theory provides a practical condition that can be used to calculate equilibrium constants for the cluster. Hill's theory was applied in a molecular dynamics (MD) study of clusters in supercritical fluids by Petsche and Debenedetti (1989). In that study, the first dynamic pictures of supercritical fluid clusters were presented for Lennard-Jones (L J) mixtures. However, the method had some difficulties in determining the time correlation function of the cluster members because of the necessity to compute the connectivity matrix of directly and indirectly bound particles in the entire system (Sevick et al., t988). In this work, we propose a new definition of a cluster which is simple and allows one to easily follow the time evolution of statistical functions. We use the definition to investigate the dynamic structure of a cluster formed around a solute in dilute supercritical fluids with a MD simulation technique. DEFINITION OF C L U S T E R

Hill's definition of a cluster judges whether a pair of particles is bound or not, by comparing their relative kinetic energy (EK12) with their pair potential energy (U12). The particles are considered to be a bound pair, only when their EKt2 is less than or equal to U~2. EK12 (Vx,y,) 5 -U12(r )

where

(1)

C. C. Liew et al. /Fluid Phase Equilibria 104 (1995) 317-329 m I *m 2

EKl2 (Vx,y,z)

=

o.5

m 1+ m 2

v= '"

319

(2)

From the viewpoint of computer simulation, Hill's definition of a cluster is defined without any empirical parameters except the intermolecular potential energy function. Furthermore, Hill's definition is the only one that considers configuration and velocities of the particles which can give more information about the dynamic structure of a cluster. For multi-particle system, the computation of the cluster connectivity matrix requires a computation of I/2N(N-I) loops. In practice, the definition of Hill has difficulties when applied to MD simulation at higher densities. Definition based on the relative radial kinetic energy

The key problem in cluster simulation is the determination of when a specified particle can be considered to be a cluster member as opposed to a solvent particle at the specified state point. If we choose the definition of Hill, the periodic boundary condition may cause another problem when a "supercluster" is formed over periodic image boxes. To circumvent all of these difficulties we use a definition based on the relative radial kinetic energy as explained below. We consider the relative motion of two particles as shown in Fig.l. Here, the relative velocity of the particles can be divided into two components, the relative radial motion (Vr) and circular (or angular) motion (Vc).

VCe:: . V12

/7 /

/

"Solute Fig. I. Relative velocity of two particles and its component. If we consider the pairwise interaction of LJ fluids, it is clear that only the radial relative motion is affected by their pair potential. Therefore, in a mixture, the solvents in the cluster around a solute can be defined as those whose radial component of the relative kinetic energy E K r is less than or equal to their pair potential energy (U12). EK12(V~) < -U12(r)

(3)

320

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

where

EK12(V~)

= 0.5 ml *m2 V,~2

(4)

ml+m 2

In comparison with Hill's definition, our definition is based on the radial component of the relative velocity, V,., instead of the total relative velocity V,-.,..: • Furthermore, we consider cluster members to be only those particles that are directly bound, which eliminates computation of the connectivity matrix. In terms of MD simulation, this cluster definition is simple, and we can determine g(r) or even the time dependent radial distribution function of cluster members via ordinary MD simulation for further studies of cluster dynamic structure. MD SIMULATION Constant N,V,T MD simulations were performed via a Nose-Hoover thermostat (Nose, 1984; Hoover, 1985). The state point was set at T*=l.40, p*=0.35, which is slightly above the critical point of LJ fluids. One solute molecule plus 863 CO 2 molecules with parameters listed in Table 1 were used to simulate an infinitely dilute supercritical fluid solution. We chose pyrene and three modified model molecules as solutes. These form an attractive mixture as defined by Petsche and Debenedetti (1991). The model mixture allowed study of mass, size and energy parameters on cluster formation. We also simulate a pure CO2 system, where one CO2 is selected to act as a solute, for reference. TABLE 1 Ratio of LJ parameters and mass for mixture studied in this work. Solute

o/oco 2

e/eco 2

m/moo2

CO 2 Pyrene Model A Model B Model C

1.00 1.88 1.88 1.00 1.00

1.00 2.94 1.00 2.94 1.00

1.00 4.59 1.00 1.00 4.59

o12=0.5(o1+o2) ; el:=(e 1% )0.5

The equations of motion were integrated via the velocity version of Verlet (Swope et al., 1982; Allen and Tildesley, 1987; Toxvaerd, 1991 ). Interactions were truncated at 3.00 for all particles. The time step was set at 2.5 fs. Over 10,000 steps were performed in the equilibrium runs, and 100,000 steps for the production runs. In order to analyze the dynamic cluster structure, the location and velocity of all particles in the system were saved in a data file at a fixed interval during the production run. For the counting of solvent members in the solute-solvent cluster, we limited the

321

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

distance to a sphere of diameter 13.5Oco 2 centered around the solute molecule. This corresponded to the lull length of the simulation box. Petsche and Debenedetti (1989) found that the solute-solvent correlation function decayed to the unity at 6o~ojvc,,,. The cutoff used in our work was (13.5/2)Oco 2 that is greater than that value. RESULTS

In Fig. 2, the instantaneous number of CO2 molecules as cluster members around a pyrene molecule, and its fluctuation around the average are shown. The dotted line in this figure shows the number of initial member molecules which remain as members around pyrene since time t=O.

8O

i

i

l

l

I

~

i

i

I

T

i

i

l

i

I

I

i

i

i

~a 60

E ~.

4t)

@

'~,,

9*=0.35

\,

"-,,Remaining

* 20 0

Z

J 0

i

i

i

I 1

i

i

T*=l.40

i

i

I

i

i

i

i

2

I 3

i

i

J

(4

Time [psi

Fig. 2. Number of CO 2 molecules in cluster surrounding a pyrene. As time proceeds, the number of CO 2 molecules remaining in the cluster decay rapidly as has been similarly shown by Petsche and Debenedetti (1989) for the Xe + Ne system. By comparing this decay to the number of instantaneous cluster members, we can see how rapidly the members in the cluster change from time to time.

Dynamic Structures For a better understanding of the dynamic structure of the cluster, we introduce some snapshots of the cluster in time series, The tbur snapshots in Fig. 3 (A) show the instantaneous cluster members around a pyrene molecule. The number of CO 2 molecules can be matched to the instantaneous values in Fig. 2.

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

322

As shown in Fig. 3(A), the cluster members are distributed roughly around the solute as the structure changes from time to time. Fig. 3(B) shows the time evolution of the initial cluster members which remain. As can be seen, the number of members decay very rapidly in the beginning, and decreases slowly after 1 ps. Also, members which were initially bound closely to the solute remain even after 2 ps. 'U.)'

~0~ ~ '

~- ,o.~ o b , °

t = 0.0 ps

0 0 0 o 0

0

I

-

'~'

~ o

© 0

0

5

0

0

-5

[

i

i

I

i

i

i

i

I

i

I

i

i

i

i

I

i

]

i

i

i

i

I

]

i

i

i

i

I

i

i

i

i

i -5

i

i

i

6 o~ t = 0.5 ps

' ~"

8 °

0 0

-5

-5

F

5

OO-

oo

,,~

, ~-

OdP

0

°°I

°

d

t = 1.0 ps

~ o

%

o

5

i

0

0

0 ,l

0

t = 2.0 ps

0

~ o

I

0 -5

O0

0 0

io i

-5

i

i

~

o, 0 X-axis

(A)

i

i

5 r

i

I

5

i X-axis

(B)

Fig. 3. Snapshots of (A) instantaneous and (B) remaining members in time series.

'

~

C. C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

323

One question to be answered was where do the initial members of the cluster go? They might still be located around the spot that they were in, but they may no longer be a member of the cluster. They may be just at a higher energetic state caused by a gain in KEj2 above U¢2 through collisions. In order to trace the cluster members from time t=O, we tagged the initial cluster and non-cluster members and plotted the time dependent pair correlation function, or the van Hove function (Van Hove, 1954), of them. The van Hove function of the initial cluster and non-cluster members are related to the total van Hove function by

G tota(r,t)

:

G i~n~a d ~ a ~ rncrnbers(r'l)

a iniaal n o n - c l u s t e r ,,e~e~( r't)

+

(5)

Basically, we fixed the position of the observed point at the initial position of the solute. As shown in Fig. 4, it takes about 2 ps for the first solvent to reach the position where the solute has been, and the cluster members flow into the position just slightly faster than the non-cluster members. - - T

q

oluU

"r

1

1

1

Cluster

1"

l

r

T

r

members

2

t=0.0ps ....... t = l . 0 p s ..... t=2.0ps t=3.0ps t=4.0ps

oJil

m

o

Olnl

I Y/,"/

=

Non-Cluster

.~-~

0.5 0

members

0

2 4 r*(=r/~co2)

Fig. 4. The van Hove function of initial cluster and non-cluster members.

6

324

C.C. Lieu, et aL / Fluid Phase Equilibria 104 (1995) 317-329

Another question to be answered was how do the cluster and non-cluster members exchange with time as they move through the solution? For this case, we fixed the position of the observed point to the solute and tagged the initial cluster and non-cluster members as before. We then plotted the time dependent radial distribution function, gt (r), for initial cluster and non-cluster members as shown in Fig. 5. The radial distribution of cluster members at time t=O, g,=Jr) which is known as g(r), shows that nearly 9 out of 10 solvents in the first shell of solute are cluster members, while only 2/5 of solvents in the second shell are members. After about 4 ps, the first peak of both cluster and non-cluster members become equal in height. This means approximately half of the cluster members at t=O which initially located in the first solvation shell exchange with the others in the first 4 ps.

r

T

T

[

"r

T

r

1

T

T

"f

¢

@

Cluster members t=0.0ps t=l.0ps t=2.0ps t=3.0ps t=4.0ps

2

@

,,~

1

emml

om~l

0

¢

Non-Cluster members 1

1

~0.5

0 opil

o

2

4

r*(=r/oco2)

Fig. 5. The gt(r) function of initial cluster and non-cluster members.

6

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329 I

2

~

I

I

325

I

t = 0.2ps

=l @

o~

@ o~

o~

--0 "== 2

8ps

~0 2

t = 2.0ps

E

o~

[,.

0

0

2

4 r*(=r/CYco2)

Fig. 6. Time dependent radial distribution function of cluster members .judged at t=O ( members that remain in the cluster since t=O ( . . . . . ).

) and

Fig. 6 shows the time evolution of gt (r) for the cluster. The solid lines show the gt(r) o f the initial members of the cluster, and the dotted lines show the gt (r) of the initial members that still remain as members at the snapshots. The number of members

326

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

in the second shell of the solute decay most rapidly in the beginning, whereas the members in the first shell remain in the cluster for a rather long time if compared with the other shells. After about 1 ps, only solvents members which are initially located in the first shell, or the nearest neighbor solvation shell, of the solute, remain in the cluster. Cluster Distribution

The frequency distribution of the instantaneous cluster size for various LJ systems studied in this work is shown in Fig.7. The instantaneous cluster which includes the short-ranged members and the long-ranged members, in a sphere with diameter of 13.5 Oco 2 , all showed a Gaussian distribution. As shown in the figure, the average number of members CO 2 molecules surrounding the pyrene solute was found to be approximately 58, which is nearly 3.5 times that of the pure CO2. Comparing the cluster size of model A molecule with that of model B molecule, we can see that the size parameter has a much greater contribution to the cluster size than the energy parameter. On the other hand, the mass of the solute had a small effect on the cluster size as shown by the result for model C.

I

I

i

I

0.1 r

•

'

--+-

,e

t

I

Pyrene

..... • ..... m o d e l

A

model

B

model

C

e

.......• ---

£~°

I

.........~>........ CO

2

0.05 oJ=l

m

0

20

40 Cluster

60

80

Size

Fig. 7. Distribution of members CO2 molecules in a cluster based on 5000 snapshots from 100000 MD steps for various models.

C. C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

327

Cluster Stabilit 3, A final question that we wanted to answer was, can we say anything about the cluster stability and its lifetime? As pointed out earlier, the decay of the number of remaining cluster members is rapid. When we plotted the data in semi-log coordinates as in Fig. 8, the number of remaining members initially had a large decay and approached a constant exponential decay after around 1 ps in all systems. 100

--+-- P y r e n e model A model B --o-- m o d e l C

r~

t_

-

+

~

-

. . . . .

+

+ + 10 . ~ L ~ + ~ - - - ~ + ~ _

-~-

C02

- --~0~-_ o""~.~

-~.~.~

-.~.~._o-

_~_

=

Z I

I

0

I

I

I

I

I

i

I

0.5

I

1

i

t

I

t

I

I

i

i

I

1.5

Time[ps] Fig. 8. Time evolution of solvent members remaining in the cluster since t=O. TABLE 2 Time constant, half-life of cluster and number of initial solvation members. Solute

~[ 1/ps]

z 1/2[ps]

Ns.(O)

CO 2 Pyrene Model A Model B Model C

i.26 0.71 1.43 0.89 1.07

0.55 0.98 0.48 0.78 0.65

3.8 16.2 9.8 5.4 3.6

Half-life : ~ ]/2 = ln2/-c As indicated in Fig. 6, the last half of the decay corresponds to the numbers of those solvents that interacted strongly with the solute and form the solvation shell. In other words, the number of cluster members on the decay curve should correspond to the

328

C. C. Liew et al. / F l u i d Phase Equilibria 104 (1995 ) 317-329

number of solvation members bound by short-ranged forces. In order to investigate cluster life time, lines in Fig.8 were fitted on the last half of the decay curve. We used the following function:

haNs(t) = lnNs(O)-x't

(6)

where Ns(t) is the number of solvation members in the cluster, Ns(O) is the initial number of solvation members, and "~ is the time constant of the decay. The life-time of the cluster can be determined from the time constant "~ (the slope) of the lines. The results are listed in table 2. The half-life of the cluster surrounding a pyrene is around 1 ps, that is nearly 2 times longer that of the pure system. It is clear, by comparing the half-life of all models, that the LJ energy parameter is directly related to the stability of the cluster. Furthermore, the LJ energy parameter has the greatest influence on the cluster lifetime. We can also determine the initial number of solvation members, Ns(O), from Fig. 8. On the average, about 16 CO2 molecules form a solvation shell around a pyrene molecule, whereas only about four CO 2 molecules form around a CO 2 molecule in the pure system. The LJ size parameter greatly influences the number of solvation members. Other calculations show that the LJ size parameter does not change the slope of the curve and therefore has a much smaller effect on cluster stability. CONCLUSIONS We have developed a new definition of a cluster in a supercritical fluid solution which is based on comparing the relative radial kinetic energy to the pair intermolecular potential function. The definition, which is especially convenient for MD simulation, eliminates calculation of the connectivity matrix and allows one to follow the time evolution of the solute-solvent cluster statistical functions easily. Cluster snapshots are shown and related to the statistical functions. By tagging cluster and non-cluster members, it is possible to follow their time evolution and see the rapid exchange between cluster and non-cluster members. In the pyrene + CO 2 system, approximately half of the particles exchange in the first shell after 4 ps. Our simulations show the time dependent radial distribution function of the initial cluster members and the initial members which remain as member of cluster. The result shows that, only those members which initially located in the solvation shell remain as members after about 1 ps. We also found that the number of remaining initial cluster members decay very rapidly in the beginning and approach a constant exponential decay after around 1 ps in all systems. We considered the number of cluster members on the exponential decay curve is corresponding to the number of solvation members. The half-life and the solvation number of clusters can be determined from the decay rate of the number of remain initial cluster members. By varying the LJ size and energy parameter and the molecular mass, we have found that the energy parameter strongly influences the decay rate, while the size parameter

C.C. Liew et al. / Fluid Phase Equilibria 104 (1995) 317-329

329

strongly influences the number of solvation members. More in-depth research is underway on studying cluster formation in other phases and to use the technique to study mass transfer and reaction kinetics.

ACKNOWLEDGEMENTS This work was supported by a Grant-in-Aid for Scientific Research on Priority Areas (Supercritical Fluid 224,1993, No. 04238102) from the Ministry of Education, Science and Culture of Japan.

REFERENCES Allen, M.P. and Tildesley, D.J., 1987. Computer Simulation of Liquid. Clarendon Press, Oxford University Press, Chapter 3, pp71 - 109. Bright, F.V. and McNally, ME.P,, 1992. Supercritical Fluid Technology: Theoretical and applied approaches to analytical chemistry. ACS Symposium Series 488. Bruno, T.J. and Ely, J.F., 1991. Supercritical Fluid Technology: Review in modern theory and application. CRC Press. Cochran, H.D. and Lee, L.L., 1989. Solvation structure in supercritical fluid mixtures based on molecular distribution functions. In: Johnston, K.P. and Penninger, J.M.L., Supercritical Fluid Science and Technology. ACS Symposium Series 406, Chapter 3, pp27-37. Eckert, C.A., Ziger, D.H., Johnston, K.P., and Ellison, T.K., 1983. The use of partial molar volume data to evaluate equations of state for supercritical fluid mixtures. Fluid Phase Equilibria, 14: 167-175. Eckert, C.A., and Knutson, B.L., 1993. Molecular charisma in supercritical fluids. Fluid Phase Equilibria, 83: 93-100. Hill, T.L., 1955. Molecular clusters in imperfect gases. J.Chem. Phys. 23: 617-622. Hill, T.L., 1956. Statistical Mechanics: Principle and Selected Applications. McGraw-Hill, New York, Chapter 5, pp 122-178. Hoover, W.G., 1985. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. 31: 1695-1697. Johnston, K.P. and Penninger, J.M.L., 1989. Supercritical Fluid Science and Technology. ACS Symposium Series 406. Kim, S., and Johnston, K.P., 1987. Clustering in supercritical fluid mixtures. AIChE J. 33: 16031611. Nose, S,, 1984. A unified formulation of the constant temperature molecular dynamics methods. J.Chem. Phys. 81:511-519. Petsche, I.B. and Debenedetti, P.G., 1989. Solute-solvent interactions in infinitely dilute Supercritical mixtures: A molecular dynamics investigation. J.Chem Phys. 91: 7075-7084. Petsche, I.B. and Debenedetti, P.G., 1991. Influence of solute-solvent asymmetry upon the behavior of dilute supercritical mixture. J.Chem. Phys. 95: 386-399. Sevick, E.M., Monson, P.A. and Ottino, J.M., 1988. Monte Carlo calculations of cluster statistics in continuum models of composite morphology. J.Chem. Phys. 88:1198-1206. Swope, W.C., Andersen, H.C., Berens, PH. and Wilson, K.R., 1982. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water cluster. J.Chem. Phys. 76: 637-649. Toxvaerd, S., 1991. Algorithms for canonical molecular dynamics simulations. Molec. Phys. 72: 159-168. Van Hove, L., 1954. Correlations in space and time and born approximation scattering in systems of interacting particles. Phys. Rev. 95: 249-262. Wu, R.S., Lee,L.L. and Cochran, H.D.,1990. Structure of dilute supercritical solution: clustering of solvent and solute molecules and the thermodynamic effects. Ind.Eng.Chem. Res. 29: 977-988.