An Improved Monte Carlo Scheme for Simulation of Synthesis of

Langmuir 2003, 19, 6317-6328

6317

An Improved Monte Carlo Scheme for Simulation of Synthesis of Nanoparticles in Reverse Micelles Reetu Singh, M. R. Durairaj, and Sanjeev Kumar* Department of Chemical Engineering, Indian Institute of Science, Bangalore 560 012, India Received January 20, 2003. In Final Form: April 14, 2003 An improved Monte Carlo technique is presented in this work to simulate nanoparticle formation through a micellar route. The technique builds on the simulation technique proposed by Bandyopadhyaya et al. (Langmuir 2000, 16, 7139) which is general and rigorous but at the same time very computation intensive, so much so that nanoparticle formation in low occupancy systems cannot be simulated in reasonable time. In view of this, several strategies, rationalized by simple mathematical analyses, are proposed to accelerate Monte Carlo simulations. These are elimination of infructuous events, removal of excess reactant postreaction, and use of smaller micelle population a large number of times. Infructuous events include collision of an empty micelle with another empty one or with another one containing only one molecule or only a solid particle. These strategies are incorporated in a new simulation technique which divides the entire micelle population in four classes and shifts micelles from one class to other as the simulation proceeds. The simulation results, throughly tested using chi-square and other tests, show that the predictions of the improved technique remain unchanged, but with more than an order of magnitude decrease in computational effort for some of the simulations reported in the literature. A post priori validation scheme for the correctness of the simulation results has been utilized to propose a new simulation strategy to arrive at converged simulation results with near minimum computational effort.

1. Introduction Formation of nanoparticles using reverse micelles requires two reactants premicellized in different micellar solutions. The reaction occurs when the two solutions are brought in contact by mixing. Micelles are dynamic entities which collide and coalesce with each other, through Brownian motion, and break up. During intermicellar coalescence, reactants come in contact with each other leading to reaction in the micellar core. A large number of inorganic particles have been synthesized this way. Some examples include the synthesis of silver nanoparticles through the reduction of silver salts by NaBH4,1-3 semiconductor nanoparticles such as CdS by reduction of cadmium perchlorate with sodium sulfide,4 PbS and CuS nanoparticles by reduction of the metal salt with Na2S,5 and ZnS nanoparticles by reaction between Na2S and ZnSO4.6 Nanocomposites of CdS and ZnS have been prepared by Sato et al.7 and Hirai et al.8 Petit et al.9 reported the use of functionalized surfactant silver sulfosuccinate, Ag(AOT), and reducing agent NaBH4 to prepare colloidal silver particles. Nanoparticles of Fe(OH)3 have been synthesized through the reduction of Fe(NO3)2 with NH4OH.10 Nanoparticles have also been prepared by solubilizing only one of the reactants in micelles initially * To whom correspondence should be addressed. E-mail: [email protected]. (1) Manna, A.; Imae, T.; Iida, M.; Hisamatsu, N. Langmuir 2001, 17, 6000. (2) Barnickel, P.; Wokaun, A.; Sager, W.; H.-F, E. J. Colloid Interface Sci. 1991, 148, 80. (3) Dirk, L.; Hyning, V.; Zukoski, C. F. Langmuir 1998, 14, 7034. (4) Lianos, P.; Thomas, J. K. Chem. Phys. Lett. 1986, 125, 299. (5) Lianos, P.; Thomas, J. K. J. Colloid Interface Sci. 1987, 117, 505. (6) Zhang, J.; Han, B.; Liu, J.; Zhang, X.; Liu, Z.; He, J. Chem Commun 2001, 24, 2724. (7) Sato, H.; Hirai, T.; Komasawa, I. Ind. Eng. Chem. Res. 1995, 34, 2493. (8) Hirai, T.; Sato, H.; Komasawa, I. Ind. Eng. Chem. Res. 1994, 33, 3262. (9) Petit, C.; Lixon, P.; Pileni, M.-P. J. Phys. Chem. 1993, 97, 12974. (10) Pillai, V. K. Ph.D. Thesis, University of Florida, Gainesville, 1995.

and introducing the other one subsequently. CaCO3 particles have been synthesized by adding CO2 gas to the reaction mixture phase consisting of orgranic phase, micelles, and lime.11,12 Similarly, Ca(OH)2, calcium phosphate, and calcium thiophosphate nanoparticles have been prepared by adding water slowly to the corresponding reaction mixture consisting of reverse micelles and the other reactants.13-15 Modeling of the latter type of systems with the assumption of Poisson distribution of the reactants and products in the micelles has been attempted through population balance equations.16,17 The simulation of reaction systems with two types of micelle populations present initially has been carried out using Monte Carlo simulations.18,19 In these simulations, the distribution of various species and the nanoparticles, all contained in the micelles, is allowed to evolve with time. Thus, these simulations require identification of the various events taking place in the system, the rate at which they occur, and the outcome of the events. Li and Park proposed a simulation method for the formation of Fe(OH)3 particles.10 In their simulation, they considered only fusion and redispersion of micelles to alter the state of the system. The state of a micelles is described by using a continuous variable. Any product formed is assumed to (11) Kandori, K.; Kon-No, K.; Kitahara, A. J. Colloid Interface Sci. 1987, 122, 78. (12) Roman, J.; P, H.; D, F.; C, B.; F, J.; Martin, J. J. Colloid Interface Sci. 1991, 144, 324. (13) Delfort, B.; Born, M.; Chive, A.; Barre, L. Colloid Interface Sci. 1997, 189, 151. (14) Chive, A.; Delfort, B.; Barre, L. J. Colloid Interface Sci. 1997, 186, 300. (15) Chive, A.; Delfort, B.; Born, M.; Barre, L. Langmuir 1998, 14, 5355. (16) Bandyopadyaya, R.; Kumar, R.; Gandhi, K. S.; Ramkrishna, D. Langmuir 1997, 13, 3610. (17) Bandyopadhyaya, R.; Kumar, R.; Gandhi, K. S. Langmuir 2001, 17, 1015. (18) Li, Y.; Park, C.-W. Langmuir 1999, 15, 952. (19) Bandyopadhyaya, R.; Kumar, R.; Gandhi, K. S. Langmuir 2000, 16, 7139.

10.1021/la034086w CCC: $25.00 © 2003 American Chemical Society Published on Web 06/26/2003

6318

Langmuir, Vol. 19, No. 15, 2003

Singh et al.

nucleate instantly, and the molecules present in a dimer are distributed equally among the daughter micelles. Bandyopadhyaya et al.19 have brought out the limitations in the simulation scheme of Li and Park.18 Some of these are: the inability of the scheme to predict the dynamics of the system as the rates of various processes occurring in the systems are not needed to carry out the simulation, identical predictions for systems with vastly different nucleation characteristics, and tendency to predict very wide particle size distributions because nucleation is permitted for any number (in fact fractional as well) of product molecules in the micelles. Bandyopadhyaya et al.19 have implemented a new technique which is general and overcomes many of the limitations in the technique of Li and Park. Instead of continuous variables used by Li and Park, it uses discrete variables to characterize micelles so that the populations of species in micelle can be represented realistically and correctly. The nucleation is assumed to proceed at a finite rate and only in those micelles which at the least have a critical number of product molecules in them. If a micelle containing a critical number of product molecules in it fuses with another micelle without a particle in it and then redisperses, the opportunity to nucleate may be lost because of the redistribution of product molecules in daughter micelles. If this micelle fuses with a micelle containing a particle, the particle grows instantly and opportunity to nucleate is again lost. Thus, formation of new nuclei, growth of the existing particles, and redistribution of the product molecules in micelles compete with each other. The competition among them determines the particle size distribution. Equal division of solubilizate molecules among daughter micelles, employed by Li and Park, is replaced by more general binomial redistribution. Since fusion-fission of micelles and nucleation in them occur simultaneously, the Interval of Quiescence technique of Shah et al.20 has been used to simulate the dynamics of the stochastic processes. Bandyopadhyaya et al.19 have presented a detailed comparison of simulation results obtained using their technique and that of Li and Park for the experimental system of Pillai.10 They have compared the results to determine the effect of removing one limitation at a time from the technique of Li and Park. The comparisons show that incorporation of all the new features in the scheme of Bandyopadhyaya et al.19 has an effect on simulation results. A comparison of the two simulation results with the experimental data further shows that the results obtained for the technique of Bandyopadhyaya et al.19 are in better agreement. Bandyopadhyaya et al.19 have also presented simulation results for the experimental system of Lianos and Thomas4 in which more parameters have been varied. Most of the effects studied in this experimental study could not have been captured with the simulation technique of Li and Park because of its many limitations. The comparisons provided by Bandyopadhyaya et al.19 for the experimental data of Lianos and Thomas4 show that these observations are also explained well by their simulation technique. Both Li and Park and Bandyopadhyaya et al. could obtain converged simulation results with the number of micelles as large as 100000. The technique of Bandyopadhyaya et al. makes rigorous and correct predictions for the dynamics of the process and the particle size distribution, but the use of a large number of micelles renders it extremely computation intensive. The simula-

tion results obtained for the CdS system are barely converged (shown later). The simulation of formation of CdS particles having a mean aggregation number of the order of 20 with a mean occupancy of the limiting reactant in the micelles at 0.3846 required 100000 micelles with one simulation to predict the size distribution for this system. A single simulation itself required a computation time of 4.0045 × 104 s (about 12 h) on a Compaq workstation running at 667 MHz. Pileni et al.21 have reported synthesis of silver particles using similar reverse micellar systems. The mean occupancy of the limiting reactant here is much smaller at 0.0675 and the mean aggregation numbers are of the order of 5000. Clearly, the minimum number of micelles needed to capture the whole size distribution should be well in excess of one million. Since the computation time increases as N2 with increase in number of micelles used in a simulation (shown later), a single simulation of silver particle formation using the simulation scheme of Bandyopadhyaya et al.19 can require more than 1200 h on the same machine. The objective of the present work is to develop an improved Monte Carlo technique which can be used to carry out rigorous simulations of nanoparticle formation in a realistic time frame. It is expected that development of such a technique will be useful in other areas as well.

(20) Shah, B. H.; Ramakrishna, D.; Borwanker, J. D. AIChE J. 1977, 23, 897.

(21) Pileni, M. P.; Lisiecki, I.; Motte, L.; Cizeron, J.; Moumen, N.; P, L. Prog. Colloid Polym. Sci. 1993, 93, 1.

2. Improved Monte Carlo Simulation Technique Monte Carlo simulation techniques have been used in the past in a wide variety of situations, involving a whole spectrum of computational effort. At one end of the spectrum lie problems in which every event occurring in the system affects the fate of the system at later times and needs to be incorporated in the simulation. An example of such problems is transfer of solute from dispersed phase drops to external phase with simultaneous breakup and coalesce of drops. Breakup and coalescence of drops change interfacial area for mass transport and also the state of mixing in the drops which affects internal mass transfer coefficient. Each event therefore needs to be accounted for in the simulation. At the other end of the spectrum lie simulations of those systems in which several processes occur simultaneously and the progress of some processes can be obtained without monitoring the progress of all the others. A trivial example to consider is simulation of particle breakup in turbulent liquid-liquid dispersions. Here, breakup of particles and coalescence and breakup of drops occur simultaneously. If the objective is to study particle breakup and if the latter is unaffected by the presence of liquid-liquid interface, the coalescence and breakup of drops need not be addressed at all. A more realistic example of the latter class of problems is simulation of breakup of large size particles in an agitator and their simultaneous solubilization by micelles, say, for nanoparticle formation. Both particle breakup and their dissolution need to be considered in a simulation. Since the particle breakup is dominant only for large size particles and dissolution is dominant only for small size particles, the two are nearly uncoupled parts of the problem. Particle breakup can therefore be simulated in isolation from the other process using Monte Carlo techniques. In comparison, a simulation scheme which incorporates both the processes occurring in the system at all times would be extremely computation intensive.

Simulation of Nanoparticle Synthesis

Langmuir, Vol. 19, No. 15, 2003 6319

Between the two extreme ends of the spectrumsevery event being important for the evolution of the system and the decoupling of the simultaneously occurring processess lies a class of problems in which certain types of simultaneously occurring events belonging to the same process of interest become decoupled. In general, the events that become decoupled during one stage of the process can be replaced by another set of decoupled events at a later time after the process has evolved to a new stage. Formation of nanoparticles by mixing of two micellar solutions, each one containing a reactant, has been used extensively in the recent past to produce a large variety of nanoparticles, and is a very relevant example of this class of systems. As indicated earlier, these processes have been simulated in the past by incorporating all the events in the simulation schemes, rendering the simulations extremely computation intensive and difficult to use. 2.1. Effect of Infructuous Events on System Dynamics. Monte Carlo simulation of a system can be visualized as stochastic movement of a collection of phase points (as if a cloud) in hyperdimensional space, each phase point representing a micelle present in the system. A micelle in this space is identified by all its attributes: the number of reactant molecules of A and B, the number of molecules of product C in it, and the number of product molecules present in it as a solid particle. If the system can reach a steady state, the shape of the cloud of phase points becomes invariant with time. If not, the shape continues to evolve with time. The state of the system at any instant depends on the state of the system at a previous time (Markov processes). When the system experiences some events, which do not change the distribution of phase points in the cloud, these are termed infructuous events. For the present case, the infructuous events are all those events in which the daughter micelles come out to be the same as the parent micelles. This is necessarily the case when one empty micelle coalesces with another empty micelle or those containing a single molecule or one particle in them. The removal of such infructuous events from the simulation procedure is not expected to influence the course of evolution of a Markov process. The particle size distribution obtained at the end of the process must, therefore, remain unchanged. Intuitively, it however appears that the removal of infructuous events can alter (accelerate) the dynamic evolution of the system. One can therefore expect such a simulation strategy to make the same predictions for particle size distribution as those obtained with the regular Monte Carlo simulation, but incorrectly capture the process dynamics, allowing a much too shorter process time with reduced computational effort. We show below the effect of elimination of infructuous events on prediction of process dynamics. Let us consider a process in which both the desired and infructuous events take place simultaneously, with frequencies fd and fin, respectively. The total frequency at which events occur is given by ft ) fd + fin. The interval of quiescence or the time between two successive events is a random quantity and follows Poisson distribution with mean being 1/ft. The ith random interval of quiescence can thus be obtained as20

tQi )

-ln(1 - u˜ i) ft

(1)

where u˜ i is a random number. The time required for process completion is obtained by summing up all the quiescent time intervals encountered during the process. Process completion time when only desirable events are considered

to occur in the system, and if there are Nd of them, is given by Nd

Td )

∑ i)1

-ln(1 - xi)

(2)

fd

Process completion time when both Nd desirable and Nin infructuous events occur is given by Nd+Nin

Tt )

∑ i)0

-ln(1 - xi)

(3)

ft

If Nd and Nin are very large, which is generally the case, then Nd

-

ln(1 - u˜ i) ) Nd ∫0 ∑ i)1

1

ln(1 - u) du ) Nd

(4)

Therefore, from eqs 2 and 4

Td ) Nd/fd

(5)

Similarly from eqs 3 and 4

Tt )

(Nd+Nin)f∞

lim

Nd+Nin

-

∑ i)0

-ln(1 - u˜ i) ft

)

Nd + Nin ft

(6)

The ratio of the time required for process completion for the two cases is

ftNd Td ) Tt fd(Nd + Nin)

(7)

fd Nd ) Nd + Nin ft

(8)

For Nd, Nin f ∞

Hence from eqs 7 and 8

Td/Tt ) 1

(9)

It is clear that for the interval of quiescence method, the process time remains independent of the number of events considered to take place in the system. The elimination of infructuous events from the system therefore should leave the process dynamics unaffected. An explanation for this behavior is provided by eq 1 which shows that when certain events are eliminated from the system, the frequency at which the system evolves is reduced. This in turn increases the interval of quiescence between two successive events. The events become more widely spaced in time and can now cover the same real time with a smaller number of quiescent time intervals. In a nutshell, if certain events are eliminated from the system, interval of quiescence technique allows the system to evolve with unaltered dynamics through a smaller number of widely spaced quiescent time steps. 2.2. Effect of Ignoring Decoupled Events. After the limiting reactant is exhausted from the system, the events occurring in the systemsfusion and redispersion of micelles and formation of fresh nucleisremain the same as before, but the outcome of the former becomes different. No new product molecules can form, henceforth. Prior to the exhaustion of the limiting reactant, characterization

6320


Singh et al.

of micelles required specification of one of the reactants and the presence of product molecules or nucleated particles (because two reactants cannot coexist in a micelle and a nucleated particle also cannot coexist with unnucleated product molecules). After the limiting reactant is exhausted from the system, the state of a micelle for subsequent evolution of the process can be specified merely by specifying whether it has product molecules in it or a particle. This allows a segment of micelle population to be treated as empty micelles, engaged in a large number of infructuous events, which need not be considered. The above simplification and the corresponding acceleration of the simulation can be achieved, provided it can be established that the presence of excess reactant molecules in micelles does not affect (i) the process of critical or more number of product molecules coming together to form a nuclei, (ii) growth of a nuclei, and (iii) the probability of formation of micelles containing a critical number of product molecules or more in them. The first two, because of the limited information available on nucleation and growth behavior in confined systems, are assumed to be independent of the presence of other molecules in the micelle, as is the case with nucleation and growth behavior in bulk medium at low concentrations. The third will hold if it can be shown that the redistribution of product molecules in daughter micelles is unaffected by the presence of other species in the dimer formed after fusion of the parent micelles. In an abstract sense, after the reaction reaches completion, the remaining processes still continue to move phase points in various directions of the hyperspace. The movement of phase points along the direction representing occupancy of the limiting reactant ceases because it is zero for all the micelles. If the movement of phase points along the direction representing occupancy of the excess reactant ceases to influence the movement of phase points in other directions, for example, the direction representing the occupancy of product molecules, the exchange of excess reactant then decouples from the rest of the events and can be dropped. If the likelihood of a molecule of any species in a dimer to go in one of the daughter micelles is independent of the presence of other molecules in the dimer, the redistribution of product molecules should remain unaffected by the presence of other molecules. This can be easily shown to be so. Consider a dimer formed by coalescence of two micelles containing b molecules of B and c molecules of C in it. When the dimer breaks into two daughter micelles, the total contents of the dimer are distributed binomially among the two micelles. The probability of getting a total of r molecules in a micelle after redistribution is given by

(21)

p[b + c:r] ) 2 (b+c)Cr

(b+c)

(10)

An additional factor of 2 appears because r molecules can be picked by a daughter micelle when the first micelle has either r or (b + c) - r molecules in it. The micelle containing r molecules could consist of any possible combination arising out of B and C molecules present in the dimer. The probability that they are made up of i number of B molecules and j number of C molecules, when r molecules are withdrawn from a collection of b + c molecules without distinguishing, is b

p[(i,j)|(b + c:r)] )

Ci × cCj (b+c)

Cr

(11)

The probability of getting one micelle containing i molecules of B and j molecules of C when a dimer consisting of b molecules of B and c molecules of C breaks is given by

p(b:i, c:j) ) p(b + c:r) × p[(i, j)|(b + c:r)] )

(21)

(b+c)

2 bCicCj

(12)

The probability that the daughter micelle could contain any number of C molecules but only i molecules of b is given by

()

j)c

p(b:i) ) 2

1

p(b:i, c:j) ) 2 bCi ∑ 2 j)0

b

(13)

If the species B and C are distributed among the daughter micelles independently, the probability of getting i molecules out of a pool of b molecules of B is

p(b:i) ) 2 bCi

(21)

b

(14)

Equations 13 and 14 are identical, indicating that the distribution of product molecules in the daughter micelles is independent of the presence of other species in the dimer. As the method of interval of quiescence ensures that the dynamics of evolution remains unaffected if infructuous events are dropped, and the above analysis shows that the probability of forming micelles with a given number of product molecules in it remains the same in absence of excess reactant, the simulation can be accelerated by eliminating the excess reactant postreaction. The contribution of this step to saving in computation time can vary from marginal to substantial, depending on the excess ratio used in experiments. For nearly stoichiometrically balanced reactant concentrations, the saving is negligible, whereas for systems with large excess ratio (of the order of 10), the elimination of excess reactant can decrease the computation time by more than 2 orders of magnitude. 2.3. Repeated Simulations with a Relatively Small Number of Micelles. The second-order process of fusion of micelles is central to the simulation effort. It controls both the growth of the existing particles and the buildup of supersaturation in some micelles for the nucleation to occur. The computational effort therefore increases with number of micelles N used in the simulation as CN2 (shown later). C is an undetermined constant. If the same simulation is carried out with N/M micelles M times (equal to effectively N micelles), the computational effort decreases with M as CN2/M. If the results of all M simulations carried out with N/M micelles are combined, results similar to those obtained through a single simulation with N micelles can be expected at M times smaller computational effort, provided N/M micelles can faithfully capture the process. If a single simulation, carried out with a system of 100000 micelles19 (12 h), is instead carried out with 10000 micelles 10 times and then averaged, a result of similar quality can be obtained with a 10 times reduced computational effort. These numbers suggest that one can complete the simulation in less than 5 min with a system of 500 micelles, and average taken over 200 simulations. Clearly, the system size cannot be decreased indefinitely. There must exist a certain minimum number of micelles that need to be taken to get accurate simulations results. An obvious method used in the literature to ensure accuracy is to carry out simulation with a larger number of micelles and check for independence from the system


Langmuir, Vol. 19, No. 15, 2003 6321

For coalescence between micelles belonging to the same class, the coalescence frequency is given by

βqNi(Ni - 1)N ˜ mic 2N

fii )

(18)

The factor 2 appears in the denominator to avoid double counting. The nucleation frequency is calculated the same way as by Bandyopadhyaya et al.,19 i.e., by summing up the nucleation rate for each micelle but with a slight difference. Since the micelles containing less than the critical number of product molecules or a solid particle cannot undergo nucleation, the summation is taken only over the relevant classes of micelles. Since

Figure 1. Meaningful coalescence events.

p3

size. A more effective and an internally consistency method may, however, estimate this number on the basis of the number of micelles required to capture some feature of the size distribution, say the size of the largest particle. Although the size distribution is not known, such a measure can be used to post priori verify the correctness of the simulation result. We will see in a later section if such a post priori check can be developed. 3. Implementation The simulation technique developed here builds on the Monte Carlo simulation technique of Bandyopadhyaya et al.19 A simulation starts with equal volumes of two microemulsions containing reactants A and B brought in contact at time t ) 0. The reactant species are Poisson distributed with their mean occupancies over the corresponding populations. Thus, reactant A is Poisson distributed over one-half of the population and B over the other half. The probability of finding a micelle with i molecules of a species having a mean occupancy λ for a Poisson distributed population is

λie-λ i!

Pi(λ) )

(15)

The interval of quiescence technique is used20 for simulation in conjunction with the measures considered in the previous sections to accelerate simulation. The entire population of N micelles is divided into four classes as shown in Figure 1. This division is introduced to facilitate incorporation of only meaningful coalescence events in the simulation. The infructuous events are dropped from the coalescence events. The total coalescence frequency of the events that change the state of the system is then calculated as 6

fc )

f(i) ∑ i)1

(16)

Frequencies f(1) through f(6) represent f12, f13, f23, f34, f22, and f33, respectively. Here, fij represents coalescence frequency between micelles belonging to the ith class and the jth class and is given by19

˜ mic βqNiNjN fij ) N

(17)

where q is the Brownian collision frequency, β is the coalescence efficiency, N is the total population size, and N ˜ mic is the total volume number density of the micelles.

∑

fn )

kn(li)

(19)

i)p1+1

where kn(li) is the nucleation rate of the ith micelle containing li product molecules in it. We have taken the lower limit at (p1+1) because a micelle containing a solid particle cannot undergo nucleation. The total frequency with which the system evolves is ft ) fc + fn. The probability that at any instant the next event to occur is coalescence is fc/ft. The probability for it to be a nucleation event is fn/ft. Here fn and fc both change with time. The ith quiescent time interval is calculated by

tQi ) -

ln(1 - u˜ 1) ft

(20)

where u1 is a random number obtained from a uniform random number generator U(0,1]. The event occurring at the end of a quiescent time interval is determined by generating another random number u˜ 2 in (0,1]. If u˜ 2 < fc/ft, then the event is coalescence, otherwise it is taken to be nucleation. If the event is nucleation, the nucleating micelle is chosen by generating a random number u˜ 3 in the interval. The index j satisfying the inequality j-1

∑ i)p1+1

kn(li) fn

j

< u˜ 3 e

∑ i)p1+1

kn(li) fn

(21)

is taken as the index of the micelle that nucleates and all the liquid product in it is converted to a solid particle. The scheme of Bandyopadhyaya et al.19 allowed collisions among all micelles irrespective of their state, and therefore all possible collision pairs had equal probability of coalescence. In the present scheme, only meaningful coalescence events are allowed and therefore coalescence between two micelles belonging to certain classes occur with a coalescence probability proportional to the frequency of coalescence between these classes. The probability of coalescence between classes i and j is given by

pij ) fij/fc

(22)

The same holds for coalescence among all other classes depicted in Figure 1. Thus, if coalescence is the event at hand, the classes to which the coalescing micelles belong are determined by using a random number u˜ 4 in association with the probability of coalescence among various

6322


Singh et al.

Figure 2. Flowchart for removing micelles from existing classes and adding them to suitable class after an event.

Figure 3. Flowchart for the improved Monte Carlo scheme.

classes. To do so, a cumulative probability variable F(i), defined as

1 F(i) ) fc

i

p(j) ∑ j)1

is generated using the coalescence frequencies of various classes identified in Figure 1. Here probabilities p(1) through p(6) represent probabilities p12, p13, p23, p34, p22, and p33, respectively. The index i satisfying the following inequality decides the type of the coalescence event chosen.

F(i - 1) < u˜ 4 e F(i)

(23)

The micelles from the chosen classes are picked randomly by scaling random numbers (u˜ 5 and u˜ 6) with the population size of the corresponding classes, while paying attention to the lower bound on the micelle index as shown in Figure 1. After each event the system properties are updated. The redistribution laws during fusion-redispersion are

the same as those given by Bandyopadhyaya et al.19 After an event, the state of the micelle(s), undergoing the event, is changed. The class to which the micelle(s) belongs after the event is determined. If this class happens to be the same as that of the micelle undergoing the event, the properties of the micelle are updated. If this class happens to be different from the class to which the micelle belonged before the event, the micelle has to be shifted to the new class. The micelle is then removed from the class to which it belonged prior to the event and placed in an appropriate class depending upon its new state as shown in Figure 2. The flowchart for the simulation scheme is shown in Figure 3. Before moving further, we would like to highlight an interesting point. An alternative strategy in which the first micelle is chosen randomly from any class and the second micelle is chosen randomly from only those classes with which meaningful fusion events are possible does not yield correct results. Furthermore, the simulation results change with increase in system size. 4. Results and Discussion The standard results needed to benchmark the results obtained using the present technique were generated using


Langmuir, Vol. 19, No. 15, 2003 6323

Figure 4. Particle size distribution obtained from the scheme of Bandyopadhyaya et al.19 λA ) 0.3846, λB ) 3.846, kc ) 2, N ) 105, and Nsim ) 10. Table 1. Parameters Used for Simulation parameter

value

parameter

value

b q h T Nmic

0.01 1.097 × 10-17m3 s-1 0.001 kg m-1 s-1 298 K 1.566 × 1023 m-3

B A Vm ksp Vmic R

600 278.42 s-1 5.24 × 10-29m3 3.6 × 10-23 mol2 m-6 1.8387 × 10-24 m3 32

Table 2. Standard Results Generated Using the Technique of Bandyopadhyaya et al.19 (λA ) 0.3846, λB ) 3.846, kc ) 2, N ) 105, Nsim ) 10) mean aggregation number (MAN) time for reaction completion (×104 s) process completion time (×102 s) standard deviation of MAN COV of particle size distribution time for 95% conversion to solid (×102 s)

13.84 8.09 4.01 0.0365 0.416 1.448

the technique of Bandyopadhyaya et al.19 A collection of 100000 micelles was simulated 10 times to obtain the predictions. The parameters used for simulation are given in Table 1, and the results obtained are presented in Table 2. Unless and otherwise specified, all our results are compared with these results. The results have been compared visually, as well as quantitatively using the chi-square (χ2) test of association22 to determine if the particle size distributions obtained using two techniques are statistically the same. The chi-square test is carried out on frequencies. The differences between the proportions or the fraction of particles containing 1, 2, 3, ... molecules are also plotted as a function of particle size to see if the differences between the two size distributions show any pattern consistently. The particle size distribution (corresponding to 100% conversion) obtained from the scheme of Bandyopadhyaya et al.19 is shown in Figure 4. The process dynamics is obtained by keeping track of the nucleated micelle, liquid, and nucleated product fraction as a function of time. 4.1. Repeated Simulations with a System of Smaller Size. One of the techniques proposed to accelerate the MC simulation is to consider a small population instead of a large one and take an average over a large number of simulations. A population of one thousand micelles was simulated two thousand times. The results, presented in Table 3 along with the standard results, show (22) Frank, H.; Althoen, S. C. Statistics Concepts and applications; Cambridge University Press: Cambridge, 1994.

Figure 5. Particle size distibutions: small number of simulations with large population vs large number of simulations with small population. Table 3. Comparison of Results: Small Number of Simulations with Large Population vs Large Number of Simulations with Small Population. λA ) 0.3846, λB ) 3.846, kc ) 2 process variables mean aggregation number (MAN) time for reaction completion (×104 s) process completion time (×102 s) standard deviation of MAN COV of particle size distribution time for 95% conversion to solid (×102 s)

N ) 1000 N ) 105 Nsim ) 2000 Nsim ) 10 13.82 4.75 2.66 0.119 0.434 1.456

13.84 8.09 4.01 0.037 0.416 1.448

that all the process variables match except for the reaction and process completion time. This is because the reaction completion time and the process completion time correspond to the times at which the last micelle, containing the limiting reactant molecule and the product molecule, respectively, are annihilated. In a larger system, the concentration of these micelles prior to their annihilation is smaller than that for a smaller system: lower the concentration that should be reached before their annihilation occurs, larger should be the time. We have also compared the time for 95% conversion of product molecules to solid. As expected, a fixed conversion is reached at the same time for both the populations. The particle size distribution is compared in Figure 5. As mentioned before, the quantitative tests, chi-square test of association and the differences between proportions, were carried out on the two PSDs. Although the size distributions look identical, the chi-square test shows that the hypothesis that the two size distributions are statistically similar cannot be accepted at a 5% level of significance. The estimated value of χ2 for ν ) 41 (number of bins) is 84.5, 2 (for 41 which is more than the standard value of χ41,0.05 degrees of freedom with 5% level of significance) at 56.94. The test on difference in proportions (Figure 6), however, shows that the differences in the two distributions are very small and distributed randomly. The above results are shown for 100% conversion. Simulation results were also obtained for 95% conversion to solid particles. The results show that prediction with small system size match fairly well the standard results at the same conversion. The χ2 test carried out on transient particle size distributions, shown in Figure 7, yields χ2 ) 49.45 for ν ) 35, which is smaller than 49.79, the standard 2 . The hypothesis that the two populations value of χ35,0.05 are statistically similar can be accepted at a 5% level of significance. The test on differences in proportions for the

6324


Singh et al.

Figure 6. Differences in proportions for PSDs obtained after 100% conversion: small number of simulations with large population vs large number of simulations with small population. λA ) 0.3846, λB ) 3.846, and kc ) 2.

Figure 7. Particle size distributions after 95% conversion: small number of simulations with large population vs large number of simulations with small population. λA ) 0.3846, λB ) 3.846, and kc ) 2. Table 4. Dependence of CPU Time on System Size and Number of Simulations: λA) 0.3846, λB ) 3.846, kc ) 2 varying population size

varying number of simulations

CPU time (s) N 800 1600 3200 6400 12800

Nsim ) 200 Nsim ) 500 9.93 35.77 138.96 636.96 2736.71

26.65 90.25 356.76 1585.55 7072.38

CPU time (s) Nsim

N ) 2000

N ) 5000

200 400 800 1600 3200 6400 12800

54.62 110.4 219.26 439.08 878.67 1754.5 3506.34

379.69 774.28 1552.93 3053.87 6125.41 12328.05 24515.34

two transient PSDs also gives satisfactory results, similar and still smaller in magnitude than those shown in Figure 6. It appears therefore that chi-square test for 100% conversion has amplified small differences between two large populations. We will address this issue in some detail toward the end of the paper. These comparisons show that a small population simulated a large number of times gives results similar to those obtained by simulating a large population a small number of times. Table 4 shows values of CPU time obtained for various population sizes (keeping the number of simulations constant) and different numbers of simula-

tions (keeping the number of micelles constant). The computations were carried out on a Compaq workstation with a processing speed of 667 MHz. The processing speed has no bearing on the estimation of the exponents describing the dependence of CPU time on the two variablesssystem size and the number of simulations. The absolute values of the CPU time will be different on machines running at different processor speeds. The data presented in Table 4 clearly show that CPU time increases linearly with the number of simulations and as square of the number of micelles (system size). Thus, increasing the number of simulations by a factor of 2000 and decreasing the micelle population by a factor of 100 still decreases the computation time on the whole. The issue of optimum combination of the two variables is discussed in detail later. 4.2. Species Dependent vs Species Independent Redistribution. Given that the probability of a molecule in a dimer to choose either of daughter micelles remains 0.5, irrespective of how many other molecules of what type are present, it was shown in a previous section (section 2.2) that the redistribution of species remains the same irrespective of how they are redistributed: first distribute the total number of molecules and then determine the constitutents of daughter micelles or first indentifying the constitutents of the dimer and then distribute the constitutents in the daughter drops one at a time. The simulations were carried using the parameters listed in Table 1 to confirm that the final particle size distribution also remains the same. The reaction is considered to be instantaneous so that only the excess reactant and the product are redistributed among the daughter micelles. The results obtained from the two redistribution methods are compared in Table 5. The results corresponding to the two redistribution methods are in agreement with each other including the reaction and process completion times (as the number of micelles in the two simulations is the same), mean aggregation number (MAN), and the coefficient of variance (COV). The two size distributions were tested quantitatively. The calculations show that for the two distributions χ2 ) 32.12 for ν 2 at ) 39, which is less than the standard value of χ39,0.5 54.56, indicating that the two populations can be stated to be statistically the same at a 5% level of significance. The same conclusion is reached by looking at the difference in proportions. Any differences between the two size distributions can therefore be attributed to the stochastic nature of the process. While all the variables pertaining to the converged results are identical, the standard deviation of the MAN itself, which gives an indication of simulation to simulation variation, is more for the species dependent redistribution. At this stage, the cause for this behavior is not known. It nonetheless indicates that the mimimum number of simulations needed to reached converged results for this case must be larger. 4.3. Effect of Elimination of Infructuous Events. Table 6 shows simulation results for two techniques, one incorporating all the events and the other incorporating only those that change the state of the system. The various infructuous events listed in section 3 were eliminated from the latter, and the simulations were carried out using the technique presented in this work employing the concept of classes. The excess reactant was not removed in the postreaction phase for both the simulations. The particle size distributions were tested quantitatively. The chisquare test gives χ2 ) 47.8 for ν ) 39, which is less than 2 at 54.56, indicating that two the standard value of χ39,0.05


Langmuir, Vol. 19, No. 15, 2003 6325

Table 5. Effect of Species Dependent and Species Independent Redistribution on Various Process Related Parameters: λA ) 0.3846, λB ) 3.846, kc ) 2, N ) 1000, Nsim ) 2000’ process variables

species dependent redistribution

species independent redistribution

mean aggregation number (MAN) reaction completion time (×104 s) process completion time (×102 s) standard deviation of MAN COV of particle size distribution time for 95% conversion to solid (×102 s)

13.90 4.73 2.68 0.12 0.235 1.47

13.79 4.49 2.66 0.084 0.237 1.46

Table 6. Comparison of Results: Effect of Removal of Infructuous Events from the Simulation, Excess Reactant Is Retained in the Postreaction Phase Also. λA ) 0.3846, λB ) 3.846, kc ) 2, N ) 1000, Nsim ) 2000 infructuous events process variables

included

eliminated

mean aggregation number (MAN) time for reaction completion (×104 s) process completion time (×102 s) standard deviation of MAN COV of particle size distribution time for 95% conversion to solid (×102 s)

13.82 4.75 2.66 0.119 0.434 1.450

13.93 4.78 2.71 0.261 0.434 1.448

Table 7. Comparison of Results: Effect of Elimination of Excess Reactant Postreaction (infructuous events have been eliminated). λA ) 0.3846, λB ) 3.846, kc ) 2, N ) 1000, Nsim ) 2000 excess reactant process variables

eliminated

retained

mean aggregation number (MAN) time for reaction completion (×104 s) process completion time (×102 s) standard deviation of MAN COV of particle size distribution time for 95% conversion to solid (×102 s) CPU time (s)

13.77 4.73 2.68 0.079 0.437 1.470 150.69

13.93 4.78 2.71 0.261 0.434 1.448 3.77924 × 103

populations can be taken to be statistically similar at a 5% level of significance. The results show that the simulation results and the dynamics of the system are not affected by the absence of events that do not change the state of the system. As observed before, simulation to simulation variation (indicated by standard deviation of MAN) is more for the simulation scheme in which infructuous events are eliminated. 4.4. Effect of Removal of Excess Reactant, Postreaction. After the reaction is completed, the redistribution of molecules of the excess reactant become decoupled from the processes affecting particle formation. Elimination of excess reactant after the limiting reactant is completely consumed renders micelles containing only the excess reactant empty and reduces the number of meaningful events in the system. This accelerates the simulation process as the quiescent steps become larger. A comparison of the results obtained with and without the removal of the excess reactant in the simulation, presented in Table 7, shows that the two simulation techniques produce identical results. Please notice that with elmination of excess reactant, simulation to simulation variation decreases significantly as indicated by a decrease in standard deviation of MAN from 0.261 to 0.079. The PSDs and the differences between the proportions for particle sizes were also analyzed. The value of χ2 for ν ) 40 is obtained to be 55.35, which is smaller than the 2 at 55.76. The differences in standard value of χ40,0.5 proportions were distributed nearly randomly. The two distributions are thus statistically identical at a 5% level of significance. These results show that simulation can be accelerated without any loss of accuracy. The saving in computation time with this technique increases signifi-

Figure 8. Differences between proportions for PSDs: comparison of results from the present technique (N ) 1000, Nsim ) 2000) with the standard results (N ) 105, Nsim ) 10). λA) 0.3846, λB ) 3.846, and kc ) 2.

cantly with an increase in the ratio of excess reactant to limiting reactant. 4.5. Accelerated Simulation vs Full Simulation. We present here simulation results for the new simulation technique in which all the three strategies presented so far to accelerate simulations are deployed and compare the results with those obtained using the technique of Bandyopadhyaya et al.19 Table 8 shows that for a given population size, the results obtained with the two techniques are identical. Different values obtained for reaction and process completion time for different population sizes are discussed earlier. The χ2 values for all the simulations presented in Table 8 are recorded in Table 9 for every pair of PSDs (one from the present technique and the other from the technique 2 of Bandyopadhyaya et al.19), along with the value of χν,0.05 . For size distributions in the pair to be statistically similar at a 5% level of significance, the former value should not exceed the latter. The table shows that except for one case, corresponding to χ2 ) 85.28 for ν ) 37, the above claim holds for every pair. The sequence of the differences between proportions for this case is shown in Figure 8. The figure shows that the differences are small and are nearly randomly distributed. A similar situation also arose on an earlier occasion (Figures 5 and 6). In both the cases, the size distributions look very similar, indeed. To probe the large values of χ2 that arose in some cases, a single case (N ) 1000, Nsim ) 2000) was run three times with different seeds for random number generation. The χ2 values for the three supposed to be identical size distributions were 37.9, 36.7, and 26.2 in comparison with 2 χ35,0.05 ) 49.79. Thus, although the distributions are accepted to be identical with 5% level of significance, if the level of significance were to be increased to 25%, the hypothesis that three distributions are identical will have to be rejected. This possibly is due to the large emphasis put on small deviations in large frequencies encountered

6326


Singh et al.

Table 8. Comparison of Results from the Present Technique with the Standard Results. λA ) 0.3846, λB ) 3.846, kc ) 2 N ) 1000 Nsim ) 2000

N ) 105 process variables

standard results Nsim ) 10

present technique Nsim ) 20

standard results

present technique

mean aggregation number (MAN) reaction completion time (× 104 s) process completion time (× 102 s) standard deviation of MAN COV of particle size distribution time for 95% conversion of solid (×102 s) CPU time (s)

13.83 8.09 4.01 0.037 0.416 1.448 4.4454 × 104 (1 simulation)

14.05 8.46 4.37 0.034 0.413 1.452 1.126 × 103 (1 simulation)

13.81 4.75 2.66 0.12 0.434 1.456 5.5562 × 103

13.77 4.73 2.68 0.0788 0.437 1.472 150.69

Table 9. χ2 Values Obtained by Comparing Various Particle Size Distributions: λA ) 0.3846, λB ) 3.846, kc ) 2 standard results results from present technique

Nsim ) 1

Nsim ) 10

Nsim ) 2000

N ) 1000 Nsim ) 2000 N ) 105 Nsim ) 1 N ) 105 Nsim ) 20

χ2 ) 30.25 2 χ39,0.05 ) 54.56 χ2 ) 37.53 2 ) 47.39 χ33,0.05 χ2 ) 29.01 2 χ34,0.05 ) 48.59

χ2 ) 85.28 2 χ37,0.05 ) 52.18 χ2 ) 42.8 2 χ34,0.05 ) 48.59 χ2 ) 42.8 2 χ35,0.05 ) 49.79

χ2 ) 50.77 2 χ43,0.05 ) 59.30

N ) 105

N ) 1000

χ2 ) 45.68 2 χ41,0.05 ) 56.94

Table 10. Effect of Various Techniques on CPU Time: λA ) 0.3846, λB ) 3.846, kc ) 2 CPU time (s) Bandyopadhyaya et al.19 scheme

Nmic ) 1000 Nsim ) 2000 Nmic ) 105 Nsim ) 1

infructuous events eliminated

all events considered excess reactant retained

excess reactant retained

5.5562 × 103

3.77924 × 103

4.00455 × 104

excess reactant eliminated 150.69 1.126 × 103

in these simulations. Acceptance or rejection of a hypothesis then should not be based on a χ2 test alone with a fixed level of significance, and some discretion can be exercised by the statistician with input from additional tests and considerations. We thus conclude that the present technique gives results similar to those obtained with full scale simulation for all cases. Although not shown here, the two simulation techniques have also been compared for temporal evolution23 of various quantities. The dynamics predicted using the present technique is in perfect agreement with the standard results. 4.6. Effect of Various Strategies on Acceleration of Simulations. The effect of various techniques suggested to accelerate the simulations was studied through their effect on CPU time, and the same are presented in Table 10. All the computations were carried out on a 667 MHz processing Compaq workstation. Although the absolute values of CPU time indeed depend on processor speed, the relative improvement seen here is expected to remain the same with any processor without parallelization. A comparison of computation times shows that Monte Carlo simulations carried out with a reduced population size, averaged over a large number of simulations, elimination of infructuous events, and the removal of excess reactant from the system in the postreaction phase (23) Singh, R. Modeling and simulation of batch reactors for the synthesis of nanoparticles using micellar route. Master’s thesis, Chemical Engineering, Indian Institute of Science, Bangalore, 2002.

Table 11. Comparison of the Present Technique with that of Bandyopadhyaya et al.19 for Parameters Listed in Table 1 λA ) 2, λB ) 10, kc ) 5 Nmic ) 10000, Nsim ) 50 process variables

Bandyopadhyaya et al.19

present code

mean aggregation number (MAN) reaction completion time (×104 s) process completion time (×102 s) COV of particle size distribution CPU time (s)

133.9 6.67 7.87 0.085 3.45 × 104

132.5 6.47 7.66 0.083 1.295 × 103

reduce the computation time by about 40 times without any loss of accuracy in capturing the system dynamics and the final size distribution. A major chunk of time saving however came from the removal of the excess reactant alone after exhaustion of the limiting reactant and then considering the coalescence of empty micelles as infructuous events. In the case studied by us, the initial ratio of excess reactant to limiting reactant was 10. For systems where this ratio is higher, a still larger saving in time can be expected. An example of such a system is formation of silver particles,21 mentioned earlier with initial occupancy of the limiting reactant 0.0675 and that of Ag+ salt as 8.1, making the excess reactant ratio 120! We have also carried out computations for another system. This system has higher initial mean occupancy of limiting reactant, less ratio of excess to limiting reactant and higher kc as compared to the previous system. The results from the present technique, recorded in Table 11, are the same as those from the technique of Bandyopadhyaya et al.19 The table also indicates that the saving in time is enormous even for systems with high initial mean occupancy of limiting reactant and smaller excess reactant ratio. Please notice that in both the simulations, the same number of micelles have been used and results are presented after averaging over 50 simulations. If a single simulation were to yield the same results, as was done by Li and Park and Bandyopadhyaya et al.,19 the number of micelles needed will be much larger and the saving in computational time will be much larger than that reported here. 5. Minimum Number of Micelles Needed for Simulation The issue of minimum number of micelles, Nmin mic , needed for simulations in order to obtain correct converged results is addressed here. The Nmin mic is expected to depend on the system parameters such as the mean occupancies of the reactants, critical nucleation number, nucleation rate, number density of micelles, and coalescence efficiency. It may not be possible to obtain Nmin mic as a function of all these parameters. While the number of parameters affecting the minimum number of micelles needed to get converged results is very large, we propose


Langmuir, Vol. 19, No. 15, 2003 6327

Table 12. Results for Varying Micelle Populations with System Parameters Given in Table 1. λA ) 0.3846, λB ) 3.846, kc ) 2 Nmic

[AN]dmax

MAN

COV

dmax Nmic/Nmic

% error in MAN

100 200 500 700 1000 2000 105

31 43 43 46 44 42 39

11.6 12.7 13.54 13.74 13.76 13.88 14.05

0.27 0.27 0.25 0.25 0.24 0.24 0.24

0.64 0.89 2.2 2.93 4.37 9.15 493

17.4 9.6 3.6 2.2 2.06 1.21 standard

Table 13. Results for Varying Micelle Populations with System Parameters Given in Table 1. λA ) 1, λB ) 5, kc ) 3 Nmic

[AN]dmax

MAN

COV

dmax Nmic/Nmic

% error in MAN

600 700 800 900 1000 10000

116 98 97 100 86 79

34.56 34.76 34.83 34.85 35.53 35.60

0.20 0.19 0.19 0.19 0.19 0.18

2.58 3.57 4.12 4.5 5.81 63.3

2.9 2.4 2.2 2.1 1.3 standard

Table 14. Results for Varying Micelle Populations with System Parameters Given in Table 15. λA ) 144.0, λB ) 144.0, kc ) 5 Nmic

[AN]dmax

MAN

COV

dmax Nmic/Nmic

% error in MAN

50 120 150 1000 2000

2391 2499 2341 2377 2480

712 752 759.5 769.5 771.01

0.27 0.27 0.27 0.27 0.27

1.51 3.46 4.6 30.3 58.06

7.6 2.46 1.50 0.19 standard

that all these variables affect the minimum number of micelles needed through a single variablesthe size of the largest particle formed in the simulation. The Nmin mic max should scale with Ndmic , the number of micelles needed to form the largest size particle present in the system. Thus dmax Nmin mic ) CNmic

(24)

where C is an unknown constant. The number of micelles needed to form the largest size particle (of size dmax and aggregation number [AN]dmax) would be max ) 2[AN]dmax/λA Ndmic

where λA is the mean occupancy of the limiting reactant. The factor of 2 is incorporated to account for the fact that at t ) 0, half the micelles are loaded with A and the remaining half with B. This should be changed if micellar solutions are added in different proportions. The parameters used to study the effect of number of micelles (Nmic) on predictions are listed in Table 1, and the results are shown in Table 12. The last two columns of this table max decreases, the percent show that as the ratio Nmic/Ndmic error in prediction of mean aggregation number increases. The table shows that for a value of C ) 5, the error is less max become smaller than than 2%. The values of Nmic/Ndmic unity because of the randomness associated with the initial amount of limiting reactant assigned for each simulation. A significant deviation from the mean for even one simulation changes dmax to a very large value. To test this hypothesis, simulations have been carried out for two very different systems with MAN as large as 35.6 and 771 compared with 14.05 for the previous case. The simulation results are shown in Tables 13 and 14. max > 5, the The results clearly show that for Nmic/Ndmic deviation from converged results is within 2%. Nmic and

Figure 9. Semilog plot for % error in MAN vs Nmicmin/Nmicdmax: / denotes data in Table 12, + denotes data in Table 13, and O denotes data in Table 14. Table 15. Parameters Used for Studying Systems with High Initial Mean Occupancy of Limiting Reactant parameter

value

parameter

value

β q η T Nmic Vmic

0.5 1.097 × 10-11 cm3 s-1 0.01 P 308 K 4.3 × 1016 cm-3 1.045 × 10-19 L

B A Vm ksp R

670 220.0 s-1 4.0 × 10-23 cm3 2.0 × 10-31 mol2 L-2 5

MAN for these cases are widely different. To test the possibility that error in the simulation results for insufmax ficient number of micelles is a function of ratio Nmic/Ndmic , as hypothesized earlier, the error reported in the last column of Tables 12 to 14 has been plotted against the max in Figure 9, using corresponding values of Nmic/Ndmic three different symbols for three cases. The figure shows that the data fall on a single curve for three very different systems. This indirectly supports the hypothesis that the dependence of Nmin mic on a large max number of parameters can be expressed through Ndmic alone. It is important to point out that [AN] dmax does not converge as rapidly with an increase in Nmic as MAN. The convergence for [AN] dmax is in fact so slow that even for the largest number of micelles used for each case presented in Tables 12-14, simulation with a different seed can result in more variation in [AN]dmax than that represented by the last two entries in the corresponding table. A New Approach To Obtain Converged Results. In view of the above findings, a better strategy to obtain converged simulation results is to start with a small max . If the results number of micelles, Nmic, and predict Ndmic do not satisfy the post priori check discussed above, simulate with a population 10 times larger. Repeat until max satisfies the post priori check the new value of Ndmic discussed above. 6. Conclusion The Monte Carlo simulation technique suggested by Bandyopadhyaya et al.19 to simulate formation of nanoparticle through a micellar route is general but computation intensive. The computational cost can become overwhelming to such an extent that for low occupancy systems, it poses a major handicap in simulating such systems. In view of this, several strategies, rationalized by simple mathematical analyses, are proposed to accelerate Monte Carlo simulations. These are elimination of infructuous events, removal of excess reactant postreaction, and use of smaller micelle population a large

6328


number of times. These strategies are incorporated in a new simulation technique which divides the entire micelle population into four classes and shifts micelles from one class to other as the simulation proceeds. The simulation results, throughly tested using chi-square and other tests, show that the predictions of the improved technique are the same as those obtained using the technique of Bandyopadhyaya et al.19 All the three strategies suggested accelerate the simulation process. Several test cases considered in the work show more than an order of magnitude saving in computational time. A reduction in number of micelles to minimum so as to require minimum computational time for simulation does

Singh et al.

not seem possible a priori. A new post priori validation scheme for the correctness of the simulation results has been utilized to arrive at converged results with near minimum computational effort.

Acknowledgment. The authors gratefully acknowledge many useful discussions with Professor K. S. Gandhi of Department of Chemical Engineering, IISc, Bangalore. These discussions have enriched this work significantly. LA034086W

An Improved Monte Carlo Scheme for Simulation of Synthesis of

Recommend Documents