Layer-by-layer assembly of patchy particles as a route to non-trivial structures

TL;DR

This paper introduces a layer-by-layer assembly method for patchy particles using DNA functionalization, enabling the robust formation of complex periodic structures like the Double Diamond lattice through simulations.

Contribution

It presents a novel layer-by-layer protocol combined with specific patch interactions for high-quality self-assembly of complex structures, demonstrated via simulations.

Findings

Successful in silico assembly of Double Diamond lattice

Ability to convert lattice to Cubic Diamond by particle removal

Enhanced robustness and quality of self-assembled structures

Abstract

We propose a new strategy for robust high-quality self-assembly of non-trivial periodic structures out of patchy particles, and investigate it with Brownian Dynamics (BD) simulations. Its first element is the use of specific patch-patch and shell-shell interactions between the particles, that can be implemented through differential functionalization of patched and shell regions with specific DNA strands. The other key element of our approach is the use of layer-by-layer protocol that allows to avoid a formations of undesired random aggregates. As an example, we design and self-assemble in silico a version of a Double Diamond (DD) lattice in which four particle types are arranged into BCC crystal made of four FCC sub-lattices. The lattice can be further converted to Cubic Diamond (CD) by selective removal of the particles of certain types. Our results demonstrate that by combining the directionality, selectivity of interactions and the layer-by-layer protocol, a high-quality robust self-assembly can be achieved.

Figures

Equations

$$V_{PP}=-{ϵ}_{PP}{\left(\frac{2a}{r}\right)}^{6}g({\theta }_i)g({\theta }_j)$$

$$V_{SS}=-{ϵ}_{SS}{\left(\frac{2a}{r}\right)}^6(1-g({\theta }_i))(1-g({\theta }_j)),$$

$$V_{Rep}(r)=ϵe^{-\beta (\frac{r}{a}-2)}.$$

$${\varphi }_4=\frac{1}{1+(\frac{ϵ}{{ϵ}_0}{)}^{\frac{3}{2}}{e}^{-Zk(ϵ-{ϵ}_0)}},$$

$$\frac{(1-f{)}^2}{f}=e^{k({ϵ}_0-ϵ)},$$

Full text

I Introduction

In recent years, the field of nanoparticle (NP) and colloidal self-assembly had demonstrated an impressive progress. It was marked by an emergence of new approaches, such as the use of anisotropic particles and specific key-lock interactions. Among the most popular model systems with anisotropic potentials are so-called patchy particles, i.e. colloids with chemically distinct “patches” on their surfaces Zhang and Glotzer (2004); Wilber et al. (2007); Kern and Frenkel (2003); Smallenburg and Sciortino (2013); Pawar and Kretzschmar (2016); Manoharan, Elsesser, and Pine (2003). The key-lock interactions are typically introduced with the help of single-stranded DNA (ssDNA) attached to the particle. Such DNA-functionalized NPs and colloids are designed to selectively bind to each other, thanks to Watson-Crick hybridization of complementary ssDNA Nykypanchuk et al. (2008); Park et al. (2008); Martinez-Veracoechea et al. (2011). Even the simplest spherical particles with isotropic DNA-mediated interactions have been shown to form a remarkable range of crystalline lattices Macfarlane et al. (2011); Wang et al. (2015).

The major reason why the patchy particles are considered to be promising building blocks for self-assembly is that their local arrangement could be enforced by choosing a specific symmetry of the patch pattern. A classic example would be a diamond lattice, often considered to be the “Holy Grail” of colloidal self assembly thanks to its potential for fabrication of photonic bandgap materials Ho, Chan, and Soukoulis (1990); Yablonovitch, Gmitter, and Leung (1991). In a diamond, each particle is supposed to be bound to four other particles forming a tetrahedron around it. Diamond is notoriously hard to assemble, partially due to its low density: most self-assembly techniques favor compact arrangements, either to minimize interparticle potential, or to maximize system’s entropy. A particle with tetrahedral arrangement of patches is thought to be a natural building block that would favor local tetrahedral arrangement, consistent with a diamond lattice. However, both Monte Carlo (MC) and Molecular Dynamics (MD) simulations of this system reveal serious fundamental limitations of such an approach Romano, Sanz, and Sciortino (2010); Smallenburg and Sciortino (2013); Zhang et al. (2005); Vasilyev, Klumov, and Tkachenko (2013, 2015). First, there are two forms of the diamond lattices: cubic diamond (CD) and hexagonal diamond (HD), both having local tetrahedral arrangement, as shown in Figure 1. While CD is more symmetric and slightly more favorable thermodynamically, this competition would lead to a very slow self-assembly kinetics. Furthermore, CD lattice has to compete also against a liquid structure with local tetrahedral organization (so-called tetrahedral liquid) Smallenburg and Sciortino (2013). It is predicted that CD becomes a ground state of the patchy particle system only in the limit of sufficiently small patch size, i.e., very strong directionality of the bonds. As a result, the success of using patchy particles for self-assembly of a diamond or in fact any other non-trivial structure remains limited so far.

Among the encouraging recent developments is the experimental demonstration of self-assembled CD lattice in the system of DNA-functionalized NPs combined with tetrahedral “DNA cages” Liu et al. (2016). Those DNA cages are building blocks that are made with the help of DNA origami technique, and can be viewed as nanoscale patchy particles. However, it should be noted that this analogy is not straightforward: the “patches” at the vertices of the cage have anisotropic triangular shape which favors global CD arrangement. This was consistent with the earlier numerical results that highlighted an importance of restricting the rotational degree of freedom of the patch-patch interactions to achieve a robust self-assembly of a CD lattice Zhang et al. (2005); Romano and Sciortino (2012).

At present, there are several approaches to experimental implementation of colloidal patchy particles Manoharan, Elsesser, and Pine (2003); Pawar and Kretzschmar (2016). While achieving a non-circular patch shape remains challenging, the control over their overall arrangement and the particle symmetry is remarkably good. Furthermore, the patches can be decorated with DNA, thus enabling a selective DNA-based interactions between patches Wang et al. (2012). Ideally, one would like to design so-called “polychromatic” particles in which each individual patch is decorated with a specific ssDNA sequence Vasilyev, Klumov, and Tkachenko (2015). Our previous studies indicate that this would lead to a high level of control over the resulting self-assembled morphology, and in fact enable the design of a nearly arbitrary self-assembled structure Halverson and Tkachenko (2013). Unfortunately, these hypothetical polychromatic patchy particles are not yet available. Nevertheless, some level of differential DNA functionalization has been achieved. Namely, one can effectively “color” all the patches of a given particles with one type of DNA, and the rest of the particle surface (the “shell”) with another Zheng et al. (2016). This already opens a variety of new opportunities. Not only particles can be designed to have different symmetries (e.g. octahedral, tetrahedral, cubic etc.), but one can independently control the strength of patch-patch, patch-shell, and shell-shell attraction between any pair of them.

In this paper, we use Brownian Dynamics (BD) simulations to demonstrate the potential of such particles with differential DNA-functionalized patches and shells, as a platform for programmable self-assemblies of non-trivial periodic lattices. Another new element in our approach, compared to the previous numerical and experimental studies of the patchy particles, is the use of layer-by-layer growth protocol. This is a well known technique in material science and chemical engineering Richardson, Björnmalm, and Caruso (2015); Richardson et al. (2016); Xiao et al. (2016), and as we show below, it can be employed in the context of patchy particles to create periodic structures of exceptional quality.

The specific model system that we consider is motivated by the classic quest for a self-assembled diamond lattice. Similar to other proposals Zhang and Glotzer (2004); Wilber et al. (2007); Kern and Frenkel (2003); Smallenburg and Sciortino (2013); Pawar and Kretzschmar (2016); Manoharan, Elsesser, and Pine (2003), we use particles with tetrahedral patch arrangement, consistent with the diamond symmetry, but instead we target a BCC arrangement of four particle types, each forming an FCC sub-lattice. This structure can also be viewed as Double Diamond (DD) lattices, which enables its conversion to a regular diamond through partial particle deletion. Interestingly, self-assembly of DD lattice has been recently observed in conventional (non-patch) colloidal system with DNA-mediated interaction, but the mechanism behind its formation remains obscure Wang et al. (2017).

II Model and Methodology

The model structure that we target for self-assembly in this numerical study is a variation of DD lattice. The latter can be obtained by placing two types of particles: “A” and “B” into sites of a BCC lattice in such a way that each particle has four tetrahedrally arranged nearest neighbored of its type, and four more of the other type, as shown in Figure 2(a). This structure has a number of advantages over the regular CD lattice from the point of view of its self-assembly. First, it solves the problem of competition between cubic and hexagonal diamond: only the former allows an interpenetrating double lattice. Second, unlike the diamond itself, this structure is rather dense thus allowing extra interparticle contacts leading to an additional thermodynamic stability. In fact, this BCC arrangement is known to be a part of the phase diagram of tetrahedral patchy particles, but it requires substantial osmotic pressure to get stabilized Smallenburg and Sciortino (2013); Romano, Sanz, and Sciortino (2010).

Rather than relying on an external pressure, one can propose a scheme in which this structure is favored by interparticle interactions: e.g. patch-patch attraction between particles of the same types (A-A or B-B), and shell-shell attraction between those of the different types (A-B), as shown in Figure 2(b). By construction, a binary system of patchy particles with tetrahedral arrangement of patches and these interaction rules, would prefer DD lattice as that would maximize the number of favorable contacts. However, these particle will also have a tendency to random aggregation, which is a common obstacle to self-assembly of high quality crystalline structures, especially in colloidal and NP systems.

In this work, in order to avoid the undesired random aggregates, we propose a layer-by-layer protocol for the self-assembly process. Specifically, we introduce a model system of four particle types: A*⋆, B⋆, A, and B. The target structure, shown in Figure. 2(c), is essentially the same as a DD, except each of the diamond sub-lattices is now made of particles of two types. Together, the particles form a BCC lattice made of four FCC sub-lattices, each containing particles of only one type. In order to achieve the self-assembly of this complicated structure, a pairwise patch-patch interaction is introduced between particle types A⋆* and B, as well as B*⋆* and A, respectively, as illustrated in Figure  2(d). A pairwise shell-shell interaction is introduced between particle types A*⋆* and A, as well as B*⋆* and B. The advantage of this interaction scheme is that now the particles can be introduced in batches: A*⋆+B⋆* or A+B, and there is no attractive interactions between the particles of the same batch (there is a regular hard core repulsion between any two particles). This opens a possibility of a very clean layer-by-layer assembly, not affected by aggregation of the free particles in the suspension.

To study the system, we introduce the force field which is analogous to Kern-Frenkel Kern and Frenkel (2003) model of patchy interactions, but replaces step-like potentials with smooth functions making it suitable for Molecular Dynamics (MD) or BD simulations. A number of conceptually similar but formally distinct force fields have been previously proposed in Refs. Miller and Cacciuto (2009); Mahynski et al. (2016). The patch-patch (PP) and shell-shell (SS) pairwise interaction potentials are modeled as:

[TABLE]

[TABLE]

Here $\epsilon_{PP}$ and $\epsilon_{SS}$, are patch-patch and shell-shell interaction parameters, respectively. $g(\theta)=\left(\frac{1}{1+e^{\alpha(p_{0}-p)}}\right)$ is the angular potential and $\alpha$ is its sharpness parameter. $g(\theta)\approx 1$ when particles i and j are oriented so that the vector joining their centers ($\vec{r}_{ij}$) intersects an attractive patches on both particles; otherwise $g(\theta)\approx 0$. $p=cos\theta=\vec{n}_{i}\ .\ \hat{r}_{ij}$, where $\theta$ is an angle between vector $\hat{r}_{ij}=\vec{r}_{ij}/|\vec{r}_{ij}|$, and the unit vector ($\vec{n}_{i}$) that points to the center of a patch from the center of particle $i$, as shown in Figure 3a. $p_{0}=cos\theta_{0}$ is a constant that determines the patch size. The fundamental length scale in our problem is the particle radius, a = 1.00. The unit time is the self-diffusion time $\tau=\frac{a^{2}}{6D}$, where $D$ is the diffusion coefficient of the particle Halverson and Tkachenko (2013). All energies are given in units of $k_{B}T$.

In addition to attraction, there is a type-independent hard core repulsive potential:

[TABLE]

Here $\epsilon=(\epsilon_{PP}+\epsilon_{SS})/2$. In our simulations we keep the relation between $\beta$ and $\alpha$ as $\alpha=4\ \beta\ \frac{1}{\sqrt{(1-(p_{0}^{2}))}}$, to ensure that the maximum force generated by the attractive and repulsive potentials are comparable. Specifically, we use $\beta=15.00$ and $\alpha=124.04$. The particles undergo translational and rotational diffusion subjected to the forces and torques derived from the above interaction potentials, Eqs. 1-3.

III Layer-by-Layer assembly

We now introduce the layer-by-layer assembly protocol. This approach is extensively used in various sub-fields of material science and chemical engineering Richardson, Björnmalm, and Caruso (2015); Richardson et al. (2016); Xiao et al. (2016). A periodic square lattice consisting of particles A*⋆* and B*⋆*, acts as a template, as shown in Figure 4(a). Two patches of each particle are oriented upward, diagonally to the square unit lattice. Particles of the template layer are fixed during the simulation. The radius of the particle, $a=1.00$ and the lattice constant is $2.50a$ (the value is chosen so that a particle would fit well in a BCC unit cell). Next, we place $n=180$ free A- and B-type particles (number ratio 1:1) inside the cubic simulation box that already contains the template layer. The box side is $L=20.0$, and all the particles are subjected to wall-like boundary constraint that effectively limits the box height to $h=0.9L$. The volume fraction of the free particles is $\delta=\frac{4}{3}\pi a^{3}n\frac{1}{L^{2}h}$ = $0.10$. The interactions between the particles are given by the rules shown in Figure 2(d), and interaction potentials Eqs. 1-3 with the following parameters: $\theta_{0}=40^{\circ}$, $\epsilon=\epsilon_{PP}=\epsilon_{SS}=8.00$. The wall at position $z=h=0.9L$ is introduced in the form of a local downward force $f_{z}=-(\frac{a}{z-h})^{12}$.

At time $t=1000$, the system reaches a steady state with a new layer of A and B particles formed on top of the template, as shown in Figure 4(c). One can see that e.g. B particles (green color) binds with A*⋆* particles (blue color) through PP interaction and with B*⋆* (red color) through SS interaction, as dictated by our design. Figure 4(d), shows the the top view of A and B particles adsorbed on top of the template layer.

At the next step, the adsorbed particles are kept while all other A and B particles are discarded, and the new batch of particles of types A*⋆* and B*⋆* is introduced. To maintain the same volume fraction, the position of the upper boundary $h$ is shifted upward by the amount $1.3a$, which is determined as a vertical separation between the first and second layer. Similarly to the previous step, the BD simulation is run for time $t=1000$ sufficient for the particles to form the new layer. After that, all free particles are again discarded, and the new batch is introduced, this time of types A and B again. This procedure can be repeated indefinitely, each time resulting in a formation of the new layer. To speed up the simulation, we “freeze” positions and orientation of all the particles that are located deeper than two layers from the surface.

In order to check the quality of the obtained structure, we calculate several quantities that work as local order parameters. First, we can determine the filling factor $\phi$ in a particular layer, which is defined as the ratio of the number of particles in that layer to the number of sites that they can fill. Second, we can determine coordination number $Z$ of each particle by counting the contacts with its neighbors, but only those which are consistent with the target structure. This allows one to calculate anothe filling factor, $\phi_{8}$, by counting the particles in each layer that have the maximum number of such contacts, $Z=8$. Naturally, this parameter is only meaningful for the inner layers but not for the top one. The results are averaged over ten independent runs. As shown in Figure 5, the overall filling factor $\phi$ is as high as $0.99$ for all the layers, whereas parameter $\phi_{8}\geq 0.98$ is also quite close to 1.

By using the layer-by-layer protocol, we self-assembled the patchy particles into a ten-layer structure as shown in Figure 6(a) (the patches are removed for better view). As intended, the simulations produced a high quality lattice that can be interpreted as DD structure made of two interpenetrating CD sub-lattices. After removing one of them formed by the particles of types B*⋆* and A, we obtain the CD lattice shown in Figure 6(b).

IV Robustness of Layer-by-Layer Approach

In order to verify the robustness of our proposed strategy, we first investigate the adsorption of the particles on top of the template layer across a wide range of the model parameters. Namely, we choose different values of $\theta_{0}$: $10^{\circ}$, $25^{\circ}$, $40^{\circ}$, and $50^{\circ}$. We also vary the interaction strength of PP and SS interaction. To do this, we set their ratio $\gamma=\frac{\epsilon_{SS}}{\epsilon_{PP}}$ at values $0.0,0.5,1.0$, or $2.0$, and gradually vary the average interation strength, $\epsilon=\frac{\epsilon_{SS}+\epsilon_{PP}}{2}$. All other conditions are kept the same as in the original simulations.

Figure 7 shows the filling factor $\phi_{4}$ for the first adsorbed layer made of A and B particles as a function of average binding strength $\epsilon$, for various values of $\gamma$, and fixed patch angle $\theta_{0}=40^{\circ}$. Similarly to $\phi_{8}$, $\phi_{4}$ accounts for the particles that have the maximum number of the favorable contacts, which in this case is $Z=4$. As one can see from the Figure, almost all the data collapse onto a single curve, with an exception of the case of $\gamma=0$ (no shell-shell interactions) which is slightly shifted. The crossover from the empty to completely filled layer occurs around $\epsilon=4$, and the curves saturate for $\epsilon>7$. This result can be described theoretically with a simple two-state model, in which each site can be either empty or occupied, and the free energy of the “occupied” state is dominated by the strength of interaction of the particle with its Z neighbors, $Zk\epsilon$. Here $k\epsilon=k(\epsilon_{s}+\epsilon_{p})/2$ is the average of the depths of SS and PP potentials, with constant $k=0.520$ determined from Figure 3(b). We can successfully fit the simulation data with a regular Gibbs-Boltzmann formula:

[TABLE]

The best fit is achieved for $\epsilon_{0}=3.5$ and $Z=3.5$, which is close to the expected value of $Z=4$. Note that the prefactor $(\frac{\epsilon}{\epsilon_{0}})^{\frac{3}{2}}$ accounts for the fact that there are three translational degrees of freedom that are localized due to binding, and the rigidity of the bonds scale as $\epsilon$. Similarly, we investigated the first layer formation for various $\theta_{0}$ while keeping $\gamma=1.0$ constant, as shown in Figure 8. Once again, most of the data collapse to the same curve, with an exception of the smallest patch size, $\theta_{0}=10^{\circ}$.

We found that once the first layer is fully saturated (for $\epsilon>7$), the complete layer-by-layer procedure yields a near perfect lattice in a relatively wide window of parameter $\epsilon$, at least up to $\epsilon=10$. However, once the interaction strength is very large, the binding becomes nearly irreversible on the time scale of the simulations. This leads to the generation of errors. Figure 9(a-d) illustrate how the procedure fails outside the preferred parameter window. On the one hand, for $\epsilon=6$, each layer is not completely filled resulting in a large number of vacancies. On the other hand, in the limit of a very large interaction parameter, $\epsilon=20$, the structural BCC order is preserved. However, the structure becomes compositionally disordered, i.e. all the sites are occupied but not by the particles of the right type (see Figure 9(d)).

When proposing the layer-by-layer protocol and the 4-particle scheme, we argued that the simpler binary system has a generic tendency to form random aggregates during the self-assembly process. In order to demonstrate this explicitly, we simulate a system made of A and B particles with interaction rules shown in Figure 2(b). The template layer is also made of A and B particles. We add free particles by keeping volume fraction and other parameters the same as in the original simulations. As shown in Figure 10, the random aggregate indeed forms on top of the template layer. By varying the interaction strength $\epsilon$ one cannot achieve the regime in which the ordered structure is growing while the random aggregate does not. In principle, by choosing volume fraction substantially lower than in our simulations, one could expect that the ordered structure nucleated by the template layer, would require a lower value of $\epsilon$ to be formed, than the one needed for aggregation. However, that scenario appears to be unrealistic since the that volume fraction is too low and the kinetics is many orders of magnitude slower than in our scheme.

Finally, we investigate the robustness of the proposed self-assembly strategy by introducing various defects in the template layer, shown in Figure 11. We grow the layers using the above layer-by-layer protocol while keeping all the simulation conditions the same as before, and check the quality of the obtained layers by determining $\phi$ and $\phi_{8}$. As one can see in Figure 12(a), the average filling factor $\phi$ for the adsorbed layers reaches $0.98$ after just three layers are deposited, for all the defect arrangements. Similarly, $\phi_{8}$ is recovered, as shown in Figure 12(b). The obtained results show that such errors can be “healed” after depositing of just two to three layers.

V Conclusions

In this paper, we proposed and numerically investigated several strategies that allow to achieve a self-assembly of non-trivial high-quality structures. Specific model example was a BCC lattice made of four particle types, which can be viewed as a variation of DD lattices. These strategies include: (1) use of selective patch-patch and shell-shell interactions and (2) layer-by-layer assembly protocol. Thanks to a combination of the two approaches, at each stage the growing structure has an exposure to the particles that interact with its reactive surface, but not with each other. As a result, this type of self-assembly has a very low error rate. In fact, we have demonstrated that as the assembly progresses, the occasional errors get corrected and defects “healed”. The proposed scheme is certainly not limited to the specific morphology that we have explored, but this particular structure is of an additional interest since it can be converted to cubic diamond lattices upon selective deletion of half of the particles.

Not only does the layer-by-layer approach help to avoid formation of unwanted phases, such as a random aggregate, it provides a route to self-assembly of structures that do not have to be a true equilibrium state. For example, we successfully assembled DD lattice even in the regime when shell-shell interactions were completely switched off ($\gamma=0$). In that case, the attractive patch-patch interactions are acting only within pairs A-B⋆ and B-A⋆. Since the patch size is rather large, the expected equilibrium phase in both of these sub-systems would in fact be a tetrahedral liquid Smallenburg and Sciortino (2013) rather than the obtained DD-type lattice. One can imagine a variation of our approach that leads to programmable assembly of aperiodic “sandwich” structures made of even larger variety of particle types. In fact, it provides a hybrid strategy that combines the top-down fabrication with bottom-up assembly.

From the point of view of practical realization of this layer-by-layer assembly strategy, one has to solve two important experimental challenges. First, one has to develop an initiation procedure that creates the first layer for the future layer-by-layer growth. One can imagine a number of plausible approaches, e.g. use of DNA microarrays or 2D micro-patterning for positioning of those particles, or a combination of these techniques. Another technical challenge is the automation of the growth procedure: each layer is build by exposing the structure to a specific mixture of particles (either A+B or A⋆ + B⋆ in our case). Therefore, the environment has to be time-modulated. This can be achieved either by physically moving the substrate from one environment to another, or by using microfluidics. Development of such automatic protocols will also open new opportunities of manufacturing both periodic and aperiodic layered structures with new functionalities, e.g. for photonic applications.

ACKNOWLEDGMENT

The research was carried out at the Center for Functional Nanomaterials, which is a U.S. DOE Office of Science Facility, at Brookhaven National Laboratory under Contract No. DE-SC0012704.

APPENDIX: VERIFICATION AND CHARACTERIZATION OF THE MODEL

In order to verify that the force field derived from the potentials Eqs. 1-3 are working properly, we investigate the dimer formations in a binary systems of mutually attractive single-patch particles. First, we investigate the dimer formation when there is only patch-patch attraction between the particles of different types. Here, we choose $\theta_{0}=30{{}^{\circ}}$ ($p_{0}=\cos\theta_{0}=0.866$), $\epsilon=15$, an monomer volume fraction ($\delta=0.02$. We find that particles only form dimers when the $\theta\leq 30^{\circ}$. By analyzing the resulting histogram of $p_{i}=\cos\theta$, we conclude that the orientations of the particles is in the expected range for the PP interaction, $p_{i}\geq 0.866$ (see Figure 13 (black line)). As expected, we observe a very similar histogram, invert with respect to $p_{i}=0$ when only SS attraction is present, with $\theta_{0}=150{{}^{\circ}}$ (Figure 13 (red line)) .

Next, we study the dimer formation for different $\epsilon$ values where one type of interaction is present (either PP or SS). $\theta_{0}=30^{\circ}$ for PP case and $\theta_{0}=150^{\circ}$ for SS case, respectively. All other simulation conditions are same as before. Note that the SS system should behave exactly like the PP system since each particle have only one patch. Figure 14 shows the dimer formation/melting curves for different $\epsilon$ values. One can see that the melting curve for PP interaction and SS interaction overlap as expected.

A standard theoretical result for dimer-monomer coexistence can be used to fit the simulation data. At the steady state, the fraction of dimers ($f$) can be described as:

[TABLE]

where $k\epsilon_{0}=7.280$ is a constant and $k\epsilon$ is the depth of the potential. $k$ ($=0.520$) can be calculated from the Figure 3(b). A good agreement is found between the simulation data and Equation 5.