J. Phys. Chem. 1995,99, 15518-15522
15518
Molecular Dynamics Simulation of the Association of Model Colloidal Particles in Two Dimensions U. Havemann" Universitaet Potsdam, Institut fuer Computerintegrierte Systeme, Karl Liebknecht-Strasse 24/ 25, Universitaetsgelaende, Geb. IO, 14476 Golm b. Potsdam, Germany
A. G. Grivtsov (deceased) and N. N. Merkulenko Institute of Physical Chemistry, Russian Academy of Science, 1I7321 Moscow, Russia Received: March 17, 1995; In Final Form: June 29, 1995@
We investigate a simplified model of a colloidal system-a microscopic 2D system containing two "colloidal" particles and a few hundred smaller solvent molecules-by the MD method. The purpose of the present work is to study the process of association of the colloidal particles at small distances. An anisotropy of the mobility of the solvent in the vicinity of the colloidal particles is found. The potential of mean force of the two dispersed particles is also calculated.
1. Introduction In spite of many efforts both of experimentalists and theoreticians, coagulation which occurs in colloidal dispersions due to thermodynamic instability remains a phenomenon whose mechanism is not fully understood. Since the pioneering work of von Smoluchowski,' many new approaches have been attempted. In the present study a very simple model for a colloidal system is investigated by computer simulation. This method has already been successfully used to study various problems of colloidal c h e m i ~ t r y . ~ . ' ~ . ~ ' , ~ ~ From the viewpoint of statistical mechanics the process of the colloid association can be described by kinetic (stochastic) Liouville- Smoluchowski-Fokker-Plank equation^.^ Because of the structure of these equations, a numerical solution requires a dynamical Monte Carlo method,6 where usually a constant coagulation rate7s8is supposed. Another method for the theoretical description of the association of colloids is based on the direct calculation of the dynamics of the colloid particles instead of considering the time dependence of the distribution functions of these particles. The particle dynamics is calculated using the Langevin equation, which follows from the Fokker-Planck equation? Solution of the Langevin equation for the motion of the colloidal particles is achieved by the Brownian dynamicst0." method. Apparently this method is most suitable for colloidal systems. For Brownian dynamics, however, it is necessary to know a kinetic parameter: either the viscosity or the diffusion coefficient. A dependence of this parameter on the distance between the colloidal particles is due to hydrodynamics and solvent ordering near colloids, which takes place when colloidal particles are coming close to each other. In this paper we derive the potential of mean force between two large soft disks immersed in a solvent of smaller disks as well as the diffusion behavior of the solvent molecules in the vicinity of the solutes. These quantities can be used as input for a subsequent Brownian dynamics simulation.
2. Formulation of the Model Since all statistical mechanical methods pointed out above represent a simplification of the description of dispersed systems, @
Abstract published in Advance ACS Abstracrs, October 1, 1995.
an attempt is made to calculate trajectories of all interacting particles and molecules by means of numerical integration of classical equations of motion, i.e. of Newtonian equations. Thus, a molecular dynamics approach was applied to the colloidal coagulation problem.I2 We used periodic boundary conditions to make the system pseudoinfinite. The calculation gives such information as coordinates, forces, and velocities at successive time steps. In our investigation these equations describe a fairly simple system: a two-dimensional microscopic system containing two colloidal particles (or rather because of their size nuclei of colloidal particles) and a few hundred solvent molecules (Figure 1). Of course, the 2D results have at most semiquantitative predictive power for the more interesting 3D case. We have chosen a two-dimensional model for the following reasons. Firstly, visualization of the solvent structure is easier, making the essential features of the association process more transparent. Secondly, the eventual extension of these calculations to three dimensions will require very large amounts of computing time. It is therefore advisable first to test our approach in two dimensions, at the same time narrowing the region of parameter space (solute-solute distance, solute-solvent mass, temperature, etc.) that will be used in a later, more costly, 3D The results obtained with colloidal and nucleation particles are qualitatively identical; therefore, most of the calculations were performed with rather small particles to reduce the computational time. The ratio of diameters between colloid and solvent was varied from 2.0 to 5.0; the ratio of masses, from 20 to 100; and the ratio of depths of interaction potential functions, from 5 to 50. The number of solvent molecules was varied from 422 to 610; square and rectangular boxes were used, depending on the colloid particle difference. All necessary information on the details of the method of our computer experiments is given in refs 13-15. The Stoermer-Verlet algorithm was used in our simulation^^^ (time step At = 0.01tS; see below). Initially, the solvent molecules were placed at the sites of a hexagonal lattice (see the inset of Figure 1). The velocity of the solvent molecules was increased randomly in small amounts until the desired temperature was reached, after approximately 1000 time steps (Figure 2). Thereafter for several thousand steps the system was equilibrated.
0022-365419512099-155 18$09.00/0 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 42, 1995 15519
MD Simulation of the Association of Model Colloidal Particles
-48
I
1.o
.'o
I
2 .o
I
I
I
3 .O
4 .O
5 .O R/ass
Figure 3. Solvent-solvent potential Us,(see eq l), solvent-colloid potential U,, colloid-colloid potential U,,(see eq 2; parameters of the potential: b = 1, s = 20.0, a,, = 3uSr,d = 0.45,A = 5.0.
Although the potentials used in this work are much simpler then those commonly applied in colloid c h e m i ~ t r y , ~ we -~~~,'~ consider them to be sufficiently representative for special colloidal dispersion systems like highly dispersed gold in methane or neopentane, for instance.*' In the simulation, parameters of the interaction potential of solvent molecules and their mass were set to unity: os,= = m, = 1. In the chosen units the standard measurement of time is also equal to t, = u,,(ms/~ss)~'2 = 1. 2.1. Trajectory Results. A general idea of the system as well as of the behavior of the solvent near the colloidal particles at different temperatures is given in Figures 1 and 4. The positions of the solvent molecules between the 5000th and the 25OOOth step are shown (using every 50th step). A high mobility of the solvent molecules as well as an anisotropy of this mobility in the vicinity of the colloidal particles is observed. This anisotropy results from the formation of an adsorption layer around the colloidal particle. Solvent ordering of a LennardJones liquid near planar surfaces was also found in Monte-Carlo studies.l 6
s
s
-
3 .O
VI
2.5
m
2.0
1.5 1 .o
0.5 0.0 Rss/u
Figure 2. Radial distribution function g(R,,) of the solvent for reduced temperatures T = 2.0 and T = 4.0.
The solvent-solvent interaction is described by the LennardJones pair potential:
USS(RS,) = 4 ~ s , [ ( ~ , ~ R s-s )(~s,/R,361 '2
(1)
R,, being the distance between the solvent particles. In some simulations, the interaction between the colloidal particles is described by a Lennard-Jones potential, in others by a potential of the following f ~ r m : ~ ~ , ~ ~
R,, being the distance and uccthe diameter of the colloidal particles; b, s, A (Hamaker), and d are constants. The solvent-colloid interaction is also described by a Lennard-Jones potential. Parameters are calculated using the Lorentz-Berthelot rule. Figure 3 shows the interaction potentials used in the simulations discussed in sections 2.1, 3, and 4.
-8.3
-6.3
-4.3
-2.3
-0.3
1.7
3.7
5.7
7.7
./U
Figure 4. Trajectories of motion of solvent molecules in the vicinity of the colloidal particles during 20 OOO time steps (two trajectories are marked by bold lines) and reduced temperature T = 4.0.The parameters of the colloid-colloid interaction are the same as in Figure 1.
Havemann et al.
15520 J. Phys. Chem., Vol. 99, No. 42, 1995 fn fn
colloidal particles R,,(O), close to 3us, and 4u,, the distance R,,(r) remains practically stationary during the simulation. This can be considered as an indication for aggregation either in close contact or solvent separated. For Rcc(0) w 6u,, the colloid particles approach each other until a more stable aggregate structure is reached.
8.0 -
r, \
u u
a
3. Potential of Mean Force (PMF)
1.0,.
.L-.
1
1-
-
. _1
3
2
1-.
4
-1 -
-1-
5 6 [ x l O 0001
-.
1
7
steps
Figure 5. Time dependence of the distance between colloidal particles I?,,(?) for different initial values R,,(O) at temperature T = 2.0. The colloid-colloid interaction is described by the Lennard-Jones potential with a,, = 2u,, and E,, = SE,,. Number of solvent molecules, 610; number of time steps, 70 OOO; mass of colloids, m, = 20m,.
The motivation for calculating the “potential of mean force” (PMF) was 2-fold. On the one hand our intention was to determine more quantitatively the preferred aggregation distances which should coincide with the minima of the PMF.21 On the other hand the PMF is of interest as input quantity for “Brownian dynamics” simulations, especially regarding the fact that fully atomic computer simulation methods are not applicable to low concentrated systems of colloidal particles. 3.1. Calculation of the PMF. The method is based on a forced slow approach (Ar = 10-5us,) of the two colloidal particles from large distances to contact by application of slowly varying During the corresponding simulation run, the mutual force acting between the two particles along their connectivity vector 71,2 is calculated:
-F1 and
AF(r12)= (F2- F1)712/r12
(3)
F2 are the total forces exerted on particles 1 and 2 by all other solvent moleculesin the systes. The total work done against the acting force AF(r21) AF(R,,) along the line of approach of the particles is identical with a free energy change and is called the potential of mean f o r ~ e : ~ ~ , ~ ~
-b 15
(4)
u
B
10
51
! I l
l
I1
1
I
I
I
I
I
Y\
I
I /
I I
I
0
(5) 3.2. Alternative Evaluation of the Radial Pair Distribution Function. To be able to estimate the quality, error, and numerical effort of the computations described in section 3.1, the potential of mean force was also calculated using the “umbrella sampling” technique (Beme et al.,19 Rossky et al.*O). For this, the colloidal-particle interaction is superimposed by a harmonic potential V(R,,) to constrain the distance between colloidal particles to a certain distance interval:
1
/
O‘
Wr(R,,) is related to the radial pair distribution function g(R,,) via the “reversible work theorem”:I8
I
I
I
I
310
3!5
4!0
4!5
‘
5!0
I
515
Rcc/o
Figure 6. Potential of mean force: (a) resulting force AF(R,), eq 3; (b) potential of mean forces Wr(Rcc),eq 4; and (c) radial distribution function g(R,,), eq 5 , as a function of the colloid-colloid distance R,,. The solid lines are functions fitted to the computed data.
From all computer experiments it turned out that, although the solvent molecules show high mobility, structured adsorption layers can be observed. This will lead to an increase of the effective diameter of the colloidal particles moving in the dispersion. Figure 5 shows the time evolution of the distance-starting from some selected initial conditions-between the colloidal particles R,,(t) for different initial values R,,(O); the following observations can be made: For initial distances between the
The complete function g(R,,) is then obtained by superposition of the pair distributions from several slightly overlapping umbrella windows, corrected for the influence of the additional umbrella potential. 3.3. Discussion. 1. In Figure 6a,b,c AF(Rcc),Wr(Rcc),and g(R,,) gained from our computer experiments according to section 3.1 are depicted. The used parameters are contained in Table 1. We find that the maxima of the pair correlation functions have equal spacings, being equivalent to the solvent molecule diameter uss.The location of the minima of the PMF corresponds well to the aggregation distances observed in Figure 5. 2. In Figure 7 the reverse procedure is described. Here, g(Rcc)is directly determined by use of the umbrella sampling
J. Phys. Chem., Vol. 99, No. 42, I995 15521
MD Simulation of the Association of Model Colloidal Particles
2.0 1.6
1.2 0.8
.~ I
;
-.-4,.
I .
,
c 4 d
4
5 2.8
3.2
3.6
4.0
4.4
5.2
4.8
Rcc/o
Figure 7. Umbrella sampling: (a) radial distribution function g(R,); (b) potential of mean forces WdRcc),eq 7, as a function of the colloidcolloid distance Rcc.The solid lines are functions fitted to the computed data.
I
I
TABLE 1: Simulation Parameters colloidal particle diameter colloidal particle mass potential depth density average temperature number of integration steps number of solvent molecules
acc
= 3%
m, = 500ms Ec
OOV
= 20€,,
p = 0.69 T = 4.0
n = 300.000 Nc = 382
technique, while Wr(Rcc) is calculated via the “reversible work” theorem:
(7) The computational effort is-depending on the chosen force constant-a factor of 3-5 higher compared to the “direct” evaluation of Wr(Rcc). Both methods lead to the same PMF and g(Rcc). Consequently, a significant consistency can be concluded.
4. Mean Square Displacement
To give an estimation for the self-diffusion coefficients of solvent molecules in one of the following environments (see the inset of Figure 8b), molecules between colloidal particles (type 1); molecules in the vicinty of colloidal particles (type 2); or bulk solvent molecules (type 3), the mean square displacement has been calculated and averaged explicitly for molecules being in such a state. Simulations over a wide range of temperatures with different box sizes have been performed (changing from a solid-like to a highly fluid state). At low temperatures a strong adsorption behavior was found. This corresponds well to the reduced mobility of molecules of type 1 and 2 (Figure 8a). In the highly fluid state, the mobility of type 2 molecules increases up to the bulk-type value (Figure 8c), although the strong adsorption of solvent molecules between colloidal particles remains. Figure 4 gives an explanation for this behavior. The depicted trajectories show that the once adsorped molecules move
I
I
I
100
200
300
I
I
400 so0 stepxl0
Figure 8. Mean square displacement of the solvent particles for temperatures T = 2.0 (a) and T = 4.0 (b). Curves 1,2, and 3 correspond to the solvent particle types described in the inset (b). The inset illustrates the position and classification of the solvent molecules (number of solvent molecules: 233). The regions of type 1 (13 solvent molecules), type 2 (35 solvent molecules), and type 3 (185 solvent molecules) are marked. The colloid-colloid interaction is described by the Lennard-Jones potential: acc= 5ass and cCc= 5ess.
predominantly in close proximity to the colloidal particles and thus remain in their current “state”.
5. Summary The formulation and testing of the described methods to gain the “potential of mean force” from computer experiments is important for very low concentrated systems of colloidal particles. The free energy needed to seperate solvent particles from contact distance to a certain distance Rcc(t)and thus the pair correlation function g(Rcc) can easily be accessed and calculated with good accuracy. In this way a comparison between theoretical work and experiment (single particle, volume light scattering) is made possible. The “potential of mean force” method reduces the “manyparticle” problem to a simplified level, containing only the most important variables. This leads conceptually to a modified and extended stochastic Brownian dynamics approach.
Acknowledgment. One of the authors (U.H.) gratefully acknowledges the opportunity of a stay at the N.K.Balabaev group (Russian Academy of Science) for common investigations and H. Sonntag (former Academy of Science of GDR) for fruitful discussions at the beginning of the work. The authors are also obliged to G. Peinel and C. Werge (University of Leipzig, Germany) for the valuable assistance during the
15522 J. Phys. Chem., Vol. 99, No. 42, 1995
Havemann et al.
development of the MD program and R. Zahradnik and P. Hobza (Czech Academy of Science) for discussions of the colloid interaction potential. We also thank A. Geiger, D. Paschek (University Dortmund, Germany), and F. Vesely (University Vienna, Austria) for valuable discussions. Financial support from the Deutsche Forschungsgemeinschaft is gratefully acknowledged.
References and Notes (1) Von Smoluchovski, M. Z. Phys.Chem. 1918, 92, 129. (2) Shchukin, E. D.; Yushchenko, V. S. J. Materials Sci. 1981, 16, 313. ( 3 ) Overbeek, J. Th. G. In Colloid Science; Kruyt, Ed.; Elsevier: Amsterdam, 1952. (4) Gedan, H.; Lichtenfeld, H.; Sonntag, H.; Krug, H. J. Colloid Surf: 1984, 11, 199. ( 5 ) Balescu. R. Eauilibrium and Noneauilibrium Statistical Mechanics: Wileyl New York, 1475. (6) Monte Carlo Method in Statistical Phvsics; Binder, Ed.: Suringer Verlag: New York, 1979. (7) Snook, I. K.: Van Megen, W.; Gaylor, K. J.; Watts, R. 0. Adv. Colloid Interface Sci. 1982, 17, 33. (8) Van Megen, W.; Snook, I. K. Adv. Colloid Interface Sci. 1984, 21, 119. (9) Isihara, A. Statistical Physics; Academic Press: New YorkLondon, 1971. (IO) Fixman, M. J. Chem. Phys. 1978, 69, 1527, 1538. .
L
(11) Ermak, D. L.; Bucholz, H. J. Comput. Phys. 1980, 35, 169. (12) Grivtsov, A. G.; Merkulenko, N. N.; Havemann, U.; Sonntag, H. Powerchnostnye Sily; abstracts of papers of the VI11 conference; Nauka: Moscow, 1985. (13) Fedoceew, D. W.; Chuzko, R. K.: Grivtsov, A. G. Geterogennaya Kristallizatsiya iz Gazovoy Fazy; Nauka: Moscow, 1978. (14) Vesely, F. Computerexperimente an Fluessigkeitsmodellen; Phys. Verlag: Weinheim, 1978. (15) Verlet, L. Phys. Rev. 1967, 159, 98. (16) Snook, I. K.; Van Meegen, W. J. Chem. Phys. 1981, 75, 4738. (17) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Chem. Phys. 1977, 23, 327. (18) Chandler, D. Introduction to Modem Statistical Mechanics; Oxford Univ. Press: New York, 1987. (19) Beme, B. J. J. Chem. Phys. 1979, 71, 2975. (20) Rossky, P. J. Chem. Phys. Lett. 1984, 105, 577. (21) Hansen, J. P. J. Phys.: Condens. Matter 1993, 5, B117-B126. (22) Pusey, P. In Liquids, Freezing and Glass Transition; Hansen, J. P., Levesque, D., Zinn-Justin, J., Eds.; North Holland: Amsterdam, 1991. (23) Wang, C. X.; Lin, H. J.; Shi, Y. Y.; Huang, F. H. Chem. Phys. Lett. 1991, 179, 475. (24) Pearlman, D. A. J. Chem. Phys. 1993, 98, 8946. (25) Zahradnik, R. Private communications. (26) Hamaker, H. C. Physica 1937, 4 , 1058. (27) Sonntag, H. Lehrbuch der Kolloidwissenschaft; VEB Deutscher Verlag der Wissenschaften: Berlin, 1977. (28) Vesely, F. Private communications. Jp9507612
