- Email: [email protected]
Monte Carlo Metropolis simulation of interacting anisotropic polarizable spins on a lattice
Computer Physics Communications 134 (2001) 47–57 www.elsevier.nl/locate/cpc
Monte Carlo Metropolis simulation of interacting anisotropic polarizable spins on a lattice G. Danese a , I. De Lotto a , A. De Marchi a , F. Leporati a,∗ , T. Bellini b , M. Buscaglia c , F. Mantegazza d a INFM – Dipartimento di Informatica, Università di Pavia, Via Ferrata 1, 27100 Pavia, Italy b INFM – Dipartimento di Fisica, Politecnico di Milano, 20133 Milano, Italy c INFM – Dipartimento di Elettronica, Universitá di Pavia, 27100 Pavia, Italy d INFM – Dipartimento di Fisica, Università di Milano “Bicocca”, 20126 Milano, Italy
Received 2 June 2000; accepted 11 July 2000
Abstract We present a Monte Carlo Metropolis computer simulation study of the equilibrium configurations of a spin lattice model of anisotropically polarizable particles. In the presence of an external field the spins polarize and couple via dipolar interactions. To simulate the system we used both the complete Hamiltonian, requiring heavy calculations of the local field, and a simplified version, obtained by a first order expansion of the local field. We find that in both cases, the dipolar interaction induces a new and unexpected orientational order characterized by a depressed value of the orientational order parameter hP2 i and by a corresponding enhanced value of hP4 i. 2001 Elsevier Science B.V. All rights reserved. PACS: 61.20.J; 61.30; 33.55; 32.10.9; 02.50.N Keywords: Computer modeling and simulation of liquid structure; Nematic liquid crystals; Polarizability; Kerr effect; Monte Carlo Metropolis algorithm; Acceleration; Lattice calculations
1. Introduction Steady progress in computation power and simulation techniques has enabled new investigations into the origins of peculiar collective behavior in systems of interacting entities. The idea of the work reported here was suggested by an interesting phenomenon which could conceal non-trivial collective effects deriving from interacting dipoles: the recurrently observed, but never fully understood, “anomalous birefringence” of dense aqueous dispersions of electrically charged rod-like colloids [1]. DNA fragments, Tobacco Mosaic Virus and various other kinds of natural and artificial rod-like submicronic particles, show an intriguing response to external electric fields. When in diluted aqueous suspensions, these particles align themselves along the field, while when they are highly concentrated, they tend to align, on average, perpendicularly to the field. This anomalous birefringence, which has not yet * Corresponding author.
E-mail address: [email protected] (F. Leporati). 0010-4655/01/$ – see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 0 ) 0 0 1 9 0 - 9
48
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
been explained, could be a consequence of the large and anisotropic electric polarizability of the electrically charged rod-like particles surrounded by the layer of dissociated ions. We thus started to investigate the effects of dipolar interaction in systems of free to rotate anisotropic polarizable particles in presence of an external field. More specifically, we studied particles carrying no permanent dipolar moment but which were easily polarizable along the direction of their long axis. Our interest was focussed on the orientational configuration of the system. Since no investigation of this kind has been so far presented in the literature, we decided first to explore the situation where the particles are constrained on a cubic lattice. Extensions of this problem to spin lattices of various geometric configurations and to hard sphero-cylinders without positional constraints are in progress. As a further simplification we considered nearest neighbour interactions only. This work forms part of a broader series of computer simulations devoted to understanding the specific role that dipolar interactions have on the configurations and phase behavior of self-assembled structures. Although dipolar interactions are not normally considered a key element for the structuring of molecular materials, in some cases they are needed to understand the phase diagram of some liquid materials [2]. Moreover, recent computer simulations of simplified particles carrying permanent dipoles, have shown that their particular spatial and orientational order structure depends on dipolar interactions [3]. We thus executed Monte Carlo Metropolis (MM) [4] simulations of systems ranging from a few hundreds to tens of thousands of particles, with field amplitude and particle density as variable parameters. Computationally, most of the effort was devoted to the problem of calculating the local field Ei,loc at the ith lattice site, which is necessary to evaluate the induced dipole on each particle. Since Ei,loc is the sum of the external field E0 and of the dipolar fields of the neighboring particles, the calculation involves intricate implicit equations, which we computationally processed through iterative refreshing of the whole set of the electric dipoles induced on the spins [5]. Since this is a necessary step towards calculating the total energy of the system, such dipole refreshing iterations had to be performed at every Monte Carlo move. The main results presented in this paper are that: (i) dipolar interactions decrease the overall quality of the field induced spin alignment quantified by hP2 (cos ϑi )i, the standard order parameter for the non-polar orientational order, where ϑi is the angle formed by E0 and si , the unit vector directed along the ith particle; (ii) associated with the decrease in hP2 (cos ϑi )i we find a drastic increase in hP4 (cos ϑi )i, indicating a quite anomalous angular distribution, where a large fraction of the particles aligns perpendicularly to the field. The brackets h·i indicate averaging over the particle population. The precise definition of P2 (·) and P4 (·) will be given later. In this paper we show that the same behavior is obtained even after introducing an approximation in the expression for the local field, which has the advantage of enabling a direct calculation of the total energy without need for dipole refresh cycles. Adoption of this scheme considerably reduces the processing times, thus making it possible to study much larger lattices than those accessible with the unreduced Hamiltonian.
2. The model We considered a system of N spins located on a cubic lattice and characterized by a finite polarizability α|| along their axis and by a negligible perpendicular polarizability, i.e. α⊥ = 0. An external field E0 in the z direction polarizes the system. The induced dipolar moment pi on the ith spin can be expressed through the tensorial polarizability αi , pi = αi · Ei,loc . When expressed in the frame of the local coordinates (xi0 , yi0 , zi0 ), zi0 being along the spin direction, ! 0 0 0 αi = 0 0 0 ≡ α|| Ai , 0 0 α||
(1)
(2)
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
where the dimensionless tensor Ai has been introduced for future use. The local field Ei,loc is X Ej →i , Ei,loc = E0 +
49
(3)
j 6=i
where Ej→ i is the field in i due to the dipole pj . Given these definitions, the Hamiltonian of the system is: 1X (αi Ei,loc ) · Ei,loc . 2 N
H =−
(4)
i=1
Within the limit of extreme dilution, the contributions from Ej→ i become negligible and we simply obtain pi = α|| E0 cos ϑi si .
(5)
Within this limit, the system energy is exclusively due to the external electric field, so that X 1 cos2 ϑi . H0 = − α|| E02 2 N
(6)
i=1
When the density of spins is not negligible, the effect of Ej→ i should be considered. Calculating Ej→ i from the dipolar field of pj , we obtain α|| sj · Ej,loc 3(sj · uij ) · uij − sj , (7) Ej →i = 3 4πεr where ε is the dielectric constant of the medium the particles are in, r is the distance between the ith and the j th spin, uij is a unit vector along the direction of the vector joining site i with site j . If we limit our attention to nearest neighbor effects, r corresponds to the simulated lattice unit and uij is parallel to the x, y or z axis. When Eqs. (1), (3) and (7) are combined, we obtain, for a fixed orientational configuration, the amplitude of pi depending on the amplitude of all the other dipoles (which in turn depend on pi ). Such a set of equations cannot always be analytically inverted to give an explicit equation of pi as a function of external field and spin orientations. As will be discussed in the next section, this problem was solved by performing iterative calculations of pi . By inspecting this set of equations it can easily be recognized that all the physical quantities appear in the simulation through two independent dimensionless quantities only. First, we define α|| , (8) ρ= 4πεr 3 ρ is a dimensionless parameter proportional to the spin density and to its polarizability. For simplicity, in the following we will refer to ρ as “density”. Second, note that the ratio of the Hamiltonian to the thermal energy can be re-expressed as N 1 X H =− (Ai Gi,loc ) · Gi,loc , kB T 2τ
(9)
i=1
where Gi,loc is the local field when E0 = 1 (in any units) and where τ=
kB T α|| E02
(10)
is a dimensionless parameter proportional to the temperature and inversely proportional to the squared amplitude of the electric field. For simplicity, in the following we will refer to τ as “temperature”. It is possible to introduce a simplified and approximate expression for the Hamiltonian, which, as shown in the following, produces the same basic results as the one described above. By using Eqs. (3), (7) and (8), and after
50
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
some algebraic manipulation, we can obtain a concise equation relating the 3N vector Eloc = {E1,loc , E2,loc , . . .} to the 3N vector E00 = {E0 , E0 , . . .} (1−K)Eloc = E00 ,
(11)
where K is a non-dimensional operator expressing the effect of the β Cartesian component of Ej,loc on the α Cartesian component of Ei,loc . Specifically X αβ β X β Kij Ej = ρ 3(sj · uij )uα sj,β − sj,α sj,β Ej . (12) (KEloc )αi = j,β
j,β
By inverting Eq. (11) we obtain Eloc = (1−K)−1 E00 ∼ = E00 + KE00 + · · · .
(13)
If this expansion is truncated at the first order in K, i.e. assuming that ρ 1, we obtain an explicit expression of the local field as a function of spin orientation and the external field, from which an explicit expression for the energy can also be obtained X (14) 3(si · E0 )(si · uij )(sj · uij )(sj · E0 ) − (si · E0 )(si · sj )(sj · E0 ) . Hs = H0 − 2α|| ρ (ij )
Here the sum is extended over the nearest neighbor pairs. Finally, for successive use, it is helpful to introduce two parameters, the second and the fourth Legendre polynomials, that quantify the degree and quality of the orientational order of the system: (15) P2 (cos ϑi ) = 12 3cos2 ϑi −1 , (16) P4 (cos ϑi ) = 18 35 cos4 ϑi − 30 cos2 ϑi + 3 . These two quantities express the different characteristics of the orientational distribution of the system since hP2 i, also known as the nematic order parameter, is large only when the spins are aligned along the z axis, while hP4 i also has high values when a considerable part of the spins are oriented perpendicularly to the z axis. In specific physical situations, hP4 i can be expressed in terms of hP2 i. In particular, when the spins are independent of each other, but still interacting with the external field as in Eq. (5), it can be shown that 5 200 35650 hP2 i3 − hP2 i4 + · · · . (17) hP4 i ≈ f hP2 i = hP2 i2 − 7 539 49049 Eq. (16) has a range of validity which goes beyond the independent particle regime, extending to those situations where the mean effect on each lattice location of the neighboring spins can be represented as a uniform field. For example, Eq. (16) is found to be quite accurate in the case of the Lebwhol–Lasher model commonly used to describe nematic order [6].
3. The Monte Carlo Metropolis algorithm and its implementation According to the MM algorithm, equilibrium in a N particle system is reached through a sequence of moves, carried out by randomly selecting a spin, changing its orientation by a random angular jump and by evaluating the corresponding change in the energy. Each move can be accepted or rejected depending on the variation of the energy 1U associated with it [4]. To accelerate the convergence of the simulation, the amplitude of the mean angular displacement is periodically adjusted so as to obtain a mean acceptance probability of 50%. In the MM simulation reported here, the relevant properties of the system (hP2 i, hP4 i and H ) are evaluated every 2N moves, i.e. every macrostep.
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
51
Random numbers (for spin selection and move displacement) were generated using the Linear Congruential Method [7]. Since particles are fixed on lattice sites, every MM move involves a change in spin orientation. Two angles, γ and φ, are thus randomly selected: the first represents the angle between the old and the new spin orientation while the second is used to select the new orientation within the cone defined by γ . The probability distribution of φ is flat in [0, 2π], while γ is obtained as the arccosine of a number randomly distributed in the range [−δ, δ], where δ determines the maximum amplitude of the move. Periodic boundary conditions were implemented in the simulation model, the cubic box containing the particle system being replicated throughout the space so that spins adjacent to a side of the box interact with spins belonging to its periodic image, i.e. with those adjacent to the opposite side [8]. Summation of the potential contribution due to particles belonging to all the replicated boxes may be carried out either by taking into account the single interactions with neighboring particles (short-range interactions) or by resorting to asymptotic formulas, like the Ewald or the reaction field method (long-range interactions) [8]. In our case, however, only first nearest neighbor interactions were considered, testing only a few second neighbor ones since they gave approximately the same results. A critical point was the assessment of the Ei,loc consistency, since each induced dipole generates a reaction field which, in turn, modifies the original dipole. It is therefore necessary to evaluate the whole {Ei,loc } set of the local field amplitudes on all the sites. To achieve this goal, after each MM move and before the application of the Metropolis criterion, we defined an iterative procedure to evaluate the dipole strength at each spin site while keeping the spin orientation fixed. At step n of such a “dipole equilibration” routine, {Ei,loc }n were calculated from values {Ei,loc }n−1 through Eqs. (1) and (6). When the set of values {Ei,loc }n was stationary, the energy was computed and the Metropolis criterion applied. Typically, after ten cycles the system energy was defined with precision in the order of 10−3 kB T , which we took as the “convergence” requirement for the dipole equilibration process. The “equilibration” routine greatly increased computational time, which grew more than proportionally to the spin number. The increase became evident when the simulation times for the complete Hamiltonian were compared with the simulation times for the approximated Hamiltonian (Eq. (13)).
Fig. 1. Evolution of hP2 i as a function of the macrostep number, for two different initial conditions: perfectly aligned spins (continuous line); random orientations (dashed line). Complete Hamiltonian. ρ = 0.09, τ = 0.14, N = 216.
52
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
To accelerate the computation, a neighbor list was employed and a look-up table for the trigonometric terms in the spin polar coordinates was also introduced. These expedients are well known in literature, already proving their efficacy in standard MM simulations [9]. MM simulations were performed on particle systems ranging from 216 up to 30,000 at different temperatures and densities. Two different initial conditions were used in this work: the system was allowed to evolve either from an initial condition of spins ordered along the field or from a randomized configuration. The final state was found to be independent of the starting condition, as shown in Fig. 1, where we show the evolution of the parameter hP2 i as a function of the number of macrosteps in the case of the complete Hamiltonian and for the choice of parameter τ = 1/7 and ρ = 1/11. 4. Results We performed several simulations for different values of τ and ρ and for both the Hamiltonians in Eqs. (4) and (14). The results of the simulations are described through the quantities hP2 i, hP4 i and hH i/N all evaluated when the system reached its equilibrium stationary state. In Fig. 2, the dependence of the order parameter hP2 i on ρ at fixed τ is shown for both the complete and approximate Hamiltonian. This result is quite intriguing and unexpected since it indicates that the dipolar interactions have the overall effect of decreasing the average orientational order compared with the independent spin (ρ ≈ 0) system. The high density value of hP2 i is found to be in the interval 0.25–0.35. In the same manner, in Fig. 3 we show τ dependence hP2 i at fixed ρ. Here the dependence is not monotonic, with hP2 i growing after increasing τ from quite a low value (about 0.25), reaching a maximum and then decreasing to zero for large τ . This behavior supports the notion that, at low τ , under the conditions of large ρ where the dipolar interactions are expected to have the largest effect, the orientational order is unexpectedly small. To better understand the origin of this effect we show, in Fig. 4, the orientation distribution of the spins at the same temperature and for two different values of density. While the low density distribution exhibits, as expected, two peaks close to ϑ = 0 and
Fig. 2. hP2 i plotted as a function of the non-dimensional density ρ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, τ = 0.14, N = 1000) and the approximate Hamiltonian (open dots, τ = 0.1, N = 8000).
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
53
Fig. 3. hP2 i plotted as a function of the non-dimensional temperature τ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, ρ = 0.1, N = 1000) and the approximate Hamiltonian (open dots, ρ = 0.14, N = 8000).
Fig. 4. Orientational order P (θ ) of the spins with the same temperature τ = 0.14 and density ρ = 0.066 (nematic order in the upper part) or ρ = 0.09 (anomalous order in the lower part). The order is expressed in arbitrary units (a.u.).
54
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
Fig. 5. hP4 i plotted as a function of the non-dimensional density ρ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, τ = 0.14, N = 1000) and the approximate Hamiltonian (open dots, τ = 0.10, N = 8000). In the figure P4 = f (hP2 i) are portrayed for both cases (dashed line: approx. Hamilt.; continuous line: complete Hamilt.).
180, the high density distribution clearly shows that, because of dipolar interactions, a large fraction (about half of the total) of the particles lies perpendicular to the electric field. This is probably the result of the fact that dipolar interaction, which would typically favor side–side anti-pairing, cannot in this case induce anti-pairs because the dipoles are induced. As a result, to minimize the energy cost of side–side pairs, the system “prefers” to somehow alternate dipoles oriented along the field with dipoles oriented perpendicularly to it. To quantify the onset of this peculiar orientational order, we show, in Figs. 5 and 6, the dependence of hP4 i on ρ and τ . As can be seen, at low τ and large ρ, where an anomalous small hP2 i was found, here we have an enhancement of hP4 i. This can be even better appreciated if we compare hP4 i with f (hP2 i) in Eq. (16). We find, and show in Figs. 5 and 6, that, while hP4 i ≈ f (hP2 i) at low ρ or high τ , when the anomalous orientational order appears, hP4 i f (hP2 i). This fact clearly indicates that the collective behavior of an assembly of anisotropically polarizable particles cannot be mapped into an effective local uniform field, enhancing or depressing the overall orientational order, but that more complex phenomena related to pair correlations should also be taken into account. Quite interesting is also the dependence of ρ and τ on the system energy, especially when compared to HN = −α|| E02 N/2, the value of H evaluated for perfectly aligned spins. In Figs. 7 and in 8 we plot the averaged value of H /N and of Hs /N as a function, respectively, of ρ and τ . In both figures there is a range of parameters where, despite the positive contribution to the energy of the system due to thermal disorder, hH i and hHs i become smaller than HN . This fact confirms that the orientational order described by the deviation of hP4 i with respect to f (hP2 i), corresponds to a stable configuration which is energetically more favourable than the nematic-like order. Although the results obtained by studying the full vs. the approximate problem exhibit a quantitative difference, it is quite evident that their qualitative behavior is exactly the same. Hence, the simplified Hamiltonian constitutes an easier tool for studying the statistical properties of this new orientational order, such as size dependence and phase transition related fluctuations. Indeed, as anticipated, the absence of the “dipole equilibration” process considerably reduces computing times. This reduction factor scales with the particle number. Table 1 shows the computing times
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
55
Fig. 6. hP4 i plotted as a function of the non-dimensional temperature τ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, ρ = 0.1, N = 1000) and the approximate Hamiltonian (open dots, ρ = 0.14, N = 8000). In the figure P4 = f (hP2 i) are portrayed for both cases (dashed line: approx. Hamilt.; continuous line: complete Hamilt.).
per macrostep obtained for systems of a different size by using the full and the approximated Hamiltonians. The acceleration is quite evident.
5. Conclusions This paper is devoted to a study, by Monte Carlo Metropolis computer simulations, of the collective behavior of a system of anisotropically polarizable particles under the action of an external electric field. We found new unexpected and previously unobserved behavior. Our results clearly show that, as an effect of the induced dipole – induced dipole interaction, the mean alignment of the system in the direction of the field decreases. This is an important finding because it may shed new light on the electro-optic response of the colloidal suspension of rod-like particles, which, at high particle density, was found to be unexpectedly small. In our simulations, the state of the spin system is defined by two independent and dimensionless quantities qualitatively representing its temperature and density. We performed simulations as a function of both parameters. We found that, at low temperature and large density, spin angular distribution shows a peak indicating that an important fraction of the spins aligns perpendicularly to the field. The onset of such an anomalous preferred configuration is captured well by comparing the values of the two moments hP2 i and hP4 i of the angular distribution. By studying the energy of the system we showed that this anomalous orientational configuration corresponds to a stable state of the system. Further investigation is required to properly characterize the statistical properties of this system and to investigate the phase transition which is probably associated with the onset of this anomalous state. In the present study, we used two different approaches: to calculate the energy of the system we used both the complete Hamiltonian (although limited to nearest neighbor interactions) and a simplified Hamiltonian obtained under the assumption of small density. From the standpoint of the simulation, the difference between the two models is large because the complete Hamiltonian uses an implicit definition of local fields which requires an
56
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
Fig. 7. Normalised mean energy as a function of the non-dimensional density ρ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, τ = 0.4, N = 1000) and the approximate Hamiltonian (open dots, τ = 0.1, N = 8000).
Fig. 8. Normalised mean energy as a function of non-dimensional temperature τ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, ρ = 0.1, N = 1000) and the approximate Hamiltonian (open dots, ρ = 0.14, N = 8000).
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
57
Table 1 Comparison between the computation time (seconds per macrostep) needed by the complete and approximated equations for interaction energy Number of spins
Complete Hamiltonian (Eq. (4))
Approx. Hamiltonian (Eq. (14))
1000
49.04
0.023
512
13.22
0.011
256
0.747
0.005
64
0.198
0.004
iterative procedure to compute the energy, while in the simplified model the Hamiltonian is an explicit function of the spin orientation. By comparing the results of the two models it can be seen that the general behavior of the system is preserved even when the simplified approach, which drastically diminishes the computing times, is adopted.
Acknowledgements We thank Professor A. Maritan for having proposed the simplified Hamiltonian. We acknowledge useful discussions with M. Glaser, P. Pasini and C. Zannoni.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
M.E. Cates, J. Polyelectrolyte Phys. 2 (1992) 1109–1119. K. Kiyohara, K.E. Gubbins, A.Z. Panagiotopoulos, Mol. Phys. 94 (1998) 803–808. R. Berardi, S. Orlandi, C. Zannoni, Chem. Phys. Lett. 261 (1996) 357–362. N. Metropolis, A. Rosenbluth, M. Rosenbluth et al., J. Chem. Phys. 21 (1953) 1087–1092. M. Hloucha, U.K. Deiters, Fluid Phase Equil. 149 (1998) 41–56. U. Fabbri, C. Zannoni, Mol. Phys. 58 (1986) 763–788. D.E. Knuth, The Art of Computer Programming, Addison-Wesley, Reading, MA, 1981. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, 1987. P. Ewald, Ann. Phys. 64 (1921) 253–287. G. Danese, I. De Lotto, D. Dotti, F. Leporati, Comput. Physics Commun. 108 (1998) 211–217.
Monte Carlo Metropolis simulation of interacting anisotropic polarizable spins on a lattice G. Danese a , I. De Lotto a , A. De Marchi a , F. Leporati a,∗ , T. Bellini b , M. Buscaglia c , F. Mantegazza d a INFM – Dipartimento di Informatica, Università di Pavia, Via Ferrata 1, 27100 Pavia, Italy b INFM – Dipartimento di Fisica, Politecnico di Milano, 20133 Milano, Italy c INFM – Dipartimento di Elettronica, Universitá di Pavia, 27100 Pavia, Italy d INFM – Dipartimento di Fisica, Università di Milano “Bicocca”, 20126 Milano, Italy
Received 2 June 2000; accepted 11 July 2000
Abstract We present a Monte Carlo Metropolis computer simulation study of the equilibrium configurations of a spin lattice model of anisotropically polarizable particles. In the presence of an external field the spins polarize and couple via dipolar interactions. To simulate the system we used both the complete Hamiltonian, requiring heavy calculations of the local field, and a simplified version, obtained by a first order expansion of the local field. We find that in both cases, the dipolar interaction induces a new and unexpected orientational order characterized by a depressed value of the orientational order parameter hP2 i and by a corresponding enhanced value of hP4 i. 2001 Elsevier Science B.V. All rights reserved. PACS: 61.20.J; 61.30; 33.55; 32.10.9; 02.50.N Keywords: Computer modeling and simulation of liquid structure; Nematic liquid crystals; Polarizability; Kerr effect; Monte Carlo Metropolis algorithm; Acceleration; Lattice calculations
1. Introduction Steady progress in computation power and simulation techniques has enabled new investigations into the origins of peculiar collective behavior in systems of interacting entities. The idea of the work reported here was suggested by an interesting phenomenon which could conceal non-trivial collective effects deriving from interacting dipoles: the recurrently observed, but never fully understood, “anomalous birefringence” of dense aqueous dispersions of electrically charged rod-like colloids [1]. DNA fragments, Tobacco Mosaic Virus and various other kinds of natural and artificial rod-like submicronic particles, show an intriguing response to external electric fields. When in diluted aqueous suspensions, these particles align themselves along the field, while when they are highly concentrated, they tend to align, on average, perpendicularly to the field. This anomalous birefringence, which has not yet * Corresponding author.
E-mail address: [email protected] (F. Leporati). 0010-4655/01/$ – see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 0 ) 0 0 1 9 0 - 9
48
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
been explained, could be a consequence of the large and anisotropic electric polarizability of the electrically charged rod-like particles surrounded by the layer of dissociated ions. We thus started to investigate the effects of dipolar interaction in systems of free to rotate anisotropic polarizable particles in presence of an external field. More specifically, we studied particles carrying no permanent dipolar moment but which were easily polarizable along the direction of their long axis. Our interest was focussed on the orientational configuration of the system. Since no investigation of this kind has been so far presented in the literature, we decided first to explore the situation where the particles are constrained on a cubic lattice. Extensions of this problem to spin lattices of various geometric configurations and to hard sphero-cylinders without positional constraints are in progress. As a further simplification we considered nearest neighbour interactions only. This work forms part of a broader series of computer simulations devoted to understanding the specific role that dipolar interactions have on the configurations and phase behavior of self-assembled structures. Although dipolar interactions are not normally considered a key element for the structuring of molecular materials, in some cases they are needed to understand the phase diagram of some liquid materials [2]. Moreover, recent computer simulations of simplified particles carrying permanent dipoles, have shown that their particular spatial and orientational order structure depends on dipolar interactions [3]. We thus executed Monte Carlo Metropolis (MM) [4] simulations of systems ranging from a few hundreds to tens of thousands of particles, with field amplitude and particle density as variable parameters. Computationally, most of the effort was devoted to the problem of calculating the local field Ei,loc at the ith lattice site, which is necessary to evaluate the induced dipole on each particle. Since Ei,loc is the sum of the external field E0 and of the dipolar fields of the neighboring particles, the calculation involves intricate implicit equations, which we computationally processed through iterative refreshing of the whole set of the electric dipoles induced on the spins [5]. Since this is a necessary step towards calculating the total energy of the system, such dipole refreshing iterations had to be performed at every Monte Carlo move. The main results presented in this paper are that: (i) dipolar interactions decrease the overall quality of the field induced spin alignment quantified by hP2 (cos ϑi )i, the standard order parameter for the non-polar orientational order, where ϑi is the angle formed by E0 and si , the unit vector directed along the ith particle; (ii) associated with the decrease in hP2 (cos ϑi )i we find a drastic increase in hP4 (cos ϑi )i, indicating a quite anomalous angular distribution, where a large fraction of the particles aligns perpendicularly to the field. The brackets h·i indicate averaging over the particle population. The precise definition of P2 (·) and P4 (·) will be given later. In this paper we show that the same behavior is obtained even after introducing an approximation in the expression for the local field, which has the advantage of enabling a direct calculation of the total energy without need for dipole refresh cycles. Adoption of this scheme considerably reduces the processing times, thus making it possible to study much larger lattices than those accessible with the unreduced Hamiltonian.
2. The model We considered a system of N spins located on a cubic lattice and characterized by a finite polarizability α|| along their axis and by a negligible perpendicular polarizability, i.e. α⊥ = 0. An external field E0 in the z direction polarizes the system. The induced dipolar moment pi on the ith spin can be expressed through the tensorial polarizability αi , pi = αi · Ei,loc . When expressed in the frame of the local coordinates (xi0 , yi0 , zi0 ), zi0 being along the spin direction, ! 0 0 0 αi = 0 0 0 ≡ α|| Ai , 0 0 α||
(1)
(2)
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
where the dimensionless tensor Ai has been introduced for future use. The local field Ei,loc is X Ej →i , Ei,loc = E0 +
49
(3)
j 6=i
where Ej→ i is the field in i due to the dipole pj . Given these definitions, the Hamiltonian of the system is: 1X (αi Ei,loc ) · Ei,loc . 2 N
H =−
(4)
i=1
Within the limit of extreme dilution, the contributions from Ej→ i become negligible and we simply obtain pi = α|| E0 cos ϑi si .
(5)
Within this limit, the system energy is exclusively due to the external electric field, so that X 1 cos2 ϑi . H0 = − α|| E02 2 N
(6)
i=1
When the density of spins is not negligible, the effect of Ej→ i should be considered. Calculating Ej→ i from the dipolar field of pj , we obtain α|| sj · Ej,loc 3(sj · uij ) · uij − sj , (7) Ej →i = 3 4πεr where ε is the dielectric constant of the medium the particles are in, r is the distance between the ith and the j th spin, uij is a unit vector along the direction of the vector joining site i with site j . If we limit our attention to nearest neighbor effects, r corresponds to the simulated lattice unit and uij is parallel to the x, y or z axis. When Eqs. (1), (3) and (7) are combined, we obtain, for a fixed orientational configuration, the amplitude of pi depending on the amplitude of all the other dipoles (which in turn depend on pi ). Such a set of equations cannot always be analytically inverted to give an explicit equation of pi as a function of external field and spin orientations. As will be discussed in the next section, this problem was solved by performing iterative calculations of pi . By inspecting this set of equations it can easily be recognized that all the physical quantities appear in the simulation through two independent dimensionless quantities only. First, we define α|| , (8) ρ= 4πεr 3 ρ is a dimensionless parameter proportional to the spin density and to its polarizability. For simplicity, in the following we will refer to ρ as “density”. Second, note that the ratio of the Hamiltonian to the thermal energy can be re-expressed as N 1 X H =− (Ai Gi,loc ) · Gi,loc , kB T 2τ
(9)
i=1
where Gi,loc is the local field when E0 = 1 (in any units) and where τ=
kB T α|| E02
(10)
is a dimensionless parameter proportional to the temperature and inversely proportional to the squared amplitude of the electric field. For simplicity, in the following we will refer to τ as “temperature”. It is possible to introduce a simplified and approximate expression for the Hamiltonian, which, as shown in the following, produces the same basic results as the one described above. By using Eqs. (3), (7) and (8), and after
50
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
some algebraic manipulation, we can obtain a concise equation relating the 3N vector Eloc = {E1,loc , E2,loc , . . .} to the 3N vector E00 = {E0 , E0 , . . .} (1−K)Eloc = E00 ,
(11)
where K is a non-dimensional operator expressing the effect of the β Cartesian component of Ej,loc on the α Cartesian component of Ei,loc . Specifically X αβ β X β Kij Ej = ρ 3(sj · uij )uα sj,β − sj,α sj,β Ej . (12) (KEloc )αi = j,β
j,β
By inverting Eq. (11) we obtain Eloc = (1−K)−1 E00 ∼ = E00 + KE00 + · · · .
(13)
If this expansion is truncated at the first order in K, i.e. assuming that ρ 1, we obtain an explicit expression of the local field as a function of spin orientation and the external field, from which an explicit expression for the energy can also be obtained X (14) 3(si · E0 )(si · uij )(sj · uij )(sj · E0 ) − (si · E0 )(si · sj )(sj · E0 ) . Hs = H0 − 2α|| ρ (ij )
Here the sum is extended over the nearest neighbor pairs. Finally, for successive use, it is helpful to introduce two parameters, the second and the fourth Legendre polynomials, that quantify the degree and quality of the orientational order of the system: (15) P2 (cos ϑi ) = 12 3cos2 ϑi −1 , (16) P4 (cos ϑi ) = 18 35 cos4 ϑi − 30 cos2 ϑi + 3 . These two quantities express the different characteristics of the orientational distribution of the system since hP2 i, also known as the nematic order parameter, is large only when the spins are aligned along the z axis, while hP4 i also has high values when a considerable part of the spins are oriented perpendicularly to the z axis. In specific physical situations, hP4 i can be expressed in terms of hP2 i. In particular, when the spins are independent of each other, but still interacting with the external field as in Eq. (5), it can be shown that 5 200 35650 hP2 i3 − hP2 i4 + · · · . (17) hP4 i ≈ f hP2 i = hP2 i2 − 7 539 49049 Eq. (16) has a range of validity which goes beyond the independent particle regime, extending to those situations where the mean effect on each lattice location of the neighboring spins can be represented as a uniform field. For example, Eq. (16) is found to be quite accurate in the case of the Lebwhol–Lasher model commonly used to describe nematic order [6].
3. The Monte Carlo Metropolis algorithm and its implementation According to the MM algorithm, equilibrium in a N particle system is reached through a sequence of moves, carried out by randomly selecting a spin, changing its orientation by a random angular jump and by evaluating the corresponding change in the energy. Each move can be accepted or rejected depending on the variation of the energy 1U associated with it [4]. To accelerate the convergence of the simulation, the amplitude of the mean angular displacement is periodically adjusted so as to obtain a mean acceptance probability of 50%. In the MM simulation reported here, the relevant properties of the system (hP2 i, hP4 i and H ) are evaluated every 2N moves, i.e. every macrostep.
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
51
Random numbers (for spin selection and move displacement) were generated using the Linear Congruential Method [7]. Since particles are fixed on lattice sites, every MM move involves a change in spin orientation. Two angles, γ and φ, are thus randomly selected: the first represents the angle between the old and the new spin orientation while the second is used to select the new orientation within the cone defined by γ . The probability distribution of φ is flat in [0, 2π], while γ is obtained as the arccosine of a number randomly distributed in the range [−δ, δ], where δ determines the maximum amplitude of the move. Periodic boundary conditions were implemented in the simulation model, the cubic box containing the particle system being replicated throughout the space so that spins adjacent to a side of the box interact with spins belonging to its periodic image, i.e. with those adjacent to the opposite side [8]. Summation of the potential contribution due to particles belonging to all the replicated boxes may be carried out either by taking into account the single interactions with neighboring particles (short-range interactions) or by resorting to asymptotic formulas, like the Ewald or the reaction field method (long-range interactions) [8]. In our case, however, only first nearest neighbor interactions were considered, testing only a few second neighbor ones since they gave approximately the same results. A critical point was the assessment of the Ei,loc consistency, since each induced dipole generates a reaction field which, in turn, modifies the original dipole. It is therefore necessary to evaluate the whole {Ei,loc } set of the local field amplitudes on all the sites. To achieve this goal, after each MM move and before the application of the Metropolis criterion, we defined an iterative procedure to evaluate the dipole strength at each spin site while keeping the spin orientation fixed. At step n of such a “dipole equilibration” routine, {Ei,loc }n were calculated from values {Ei,loc }n−1 through Eqs. (1) and (6). When the set of values {Ei,loc }n was stationary, the energy was computed and the Metropolis criterion applied. Typically, after ten cycles the system energy was defined with precision in the order of 10−3 kB T , which we took as the “convergence” requirement for the dipole equilibration process. The “equilibration” routine greatly increased computational time, which grew more than proportionally to the spin number. The increase became evident when the simulation times for the complete Hamiltonian were compared with the simulation times for the approximated Hamiltonian (Eq. (13)).
Fig. 1. Evolution of hP2 i as a function of the macrostep number, for two different initial conditions: perfectly aligned spins (continuous line); random orientations (dashed line). Complete Hamiltonian. ρ = 0.09, τ = 0.14, N = 216.
52
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
To accelerate the computation, a neighbor list was employed and a look-up table for the trigonometric terms in the spin polar coordinates was also introduced. These expedients are well known in literature, already proving their efficacy in standard MM simulations [9]. MM simulations were performed on particle systems ranging from 216 up to 30,000 at different temperatures and densities. Two different initial conditions were used in this work: the system was allowed to evolve either from an initial condition of spins ordered along the field or from a randomized configuration. The final state was found to be independent of the starting condition, as shown in Fig. 1, where we show the evolution of the parameter hP2 i as a function of the number of macrosteps in the case of the complete Hamiltonian and for the choice of parameter τ = 1/7 and ρ = 1/11. 4. Results We performed several simulations for different values of τ and ρ and for both the Hamiltonians in Eqs. (4) and (14). The results of the simulations are described through the quantities hP2 i, hP4 i and hH i/N all evaluated when the system reached its equilibrium stationary state. In Fig. 2, the dependence of the order parameter hP2 i on ρ at fixed τ is shown for both the complete and approximate Hamiltonian. This result is quite intriguing and unexpected since it indicates that the dipolar interactions have the overall effect of decreasing the average orientational order compared with the independent spin (ρ ≈ 0) system. The high density value of hP2 i is found to be in the interval 0.25–0.35. In the same manner, in Fig. 3 we show τ dependence hP2 i at fixed ρ. Here the dependence is not monotonic, with hP2 i growing after increasing τ from quite a low value (about 0.25), reaching a maximum and then decreasing to zero for large τ . This behavior supports the notion that, at low τ , under the conditions of large ρ where the dipolar interactions are expected to have the largest effect, the orientational order is unexpectedly small. To better understand the origin of this effect we show, in Fig. 4, the orientation distribution of the spins at the same temperature and for two different values of density. While the low density distribution exhibits, as expected, two peaks close to ϑ = 0 and
Fig. 2. hP2 i plotted as a function of the non-dimensional density ρ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, τ = 0.14, N = 1000) and the approximate Hamiltonian (open dots, τ = 0.1, N = 8000).
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
53
Fig. 3. hP2 i plotted as a function of the non-dimensional temperature τ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, ρ = 0.1, N = 1000) and the approximate Hamiltonian (open dots, ρ = 0.14, N = 8000).
Fig. 4. Orientational order P (θ ) of the spins with the same temperature τ = 0.14 and density ρ = 0.066 (nematic order in the upper part) or ρ = 0.09 (anomalous order in the lower part). The order is expressed in arbitrary units (a.u.).
54
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
Fig. 5. hP4 i plotted as a function of the non-dimensional density ρ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, τ = 0.14, N = 1000) and the approximate Hamiltonian (open dots, τ = 0.10, N = 8000). In the figure P4 = f (hP2 i) are portrayed for both cases (dashed line: approx. Hamilt.; continuous line: complete Hamilt.).
180, the high density distribution clearly shows that, because of dipolar interactions, a large fraction (about half of the total) of the particles lies perpendicular to the electric field. This is probably the result of the fact that dipolar interaction, which would typically favor side–side anti-pairing, cannot in this case induce anti-pairs because the dipoles are induced. As a result, to minimize the energy cost of side–side pairs, the system “prefers” to somehow alternate dipoles oriented along the field with dipoles oriented perpendicularly to it. To quantify the onset of this peculiar orientational order, we show, in Figs. 5 and 6, the dependence of hP4 i on ρ and τ . As can be seen, at low τ and large ρ, where an anomalous small hP2 i was found, here we have an enhancement of hP4 i. This can be even better appreciated if we compare hP4 i with f (hP2 i) in Eq. (16). We find, and show in Figs. 5 and 6, that, while hP4 i ≈ f (hP2 i) at low ρ or high τ , when the anomalous orientational order appears, hP4 i f (hP2 i). This fact clearly indicates that the collective behavior of an assembly of anisotropically polarizable particles cannot be mapped into an effective local uniform field, enhancing or depressing the overall orientational order, but that more complex phenomena related to pair correlations should also be taken into account. Quite interesting is also the dependence of ρ and τ on the system energy, especially when compared to HN = −α|| E02 N/2, the value of H evaluated for perfectly aligned spins. In Figs. 7 and in 8 we plot the averaged value of H /N and of Hs /N as a function, respectively, of ρ and τ . In both figures there is a range of parameters where, despite the positive contribution to the energy of the system due to thermal disorder, hH i and hHs i become smaller than HN . This fact confirms that the orientational order described by the deviation of hP4 i with respect to f (hP2 i), corresponds to a stable configuration which is energetically more favourable than the nematic-like order. Although the results obtained by studying the full vs. the approximate problem exhibit a quantitative difference, it is quite evident that their qualitative behavior is exactly the same. Hence, the simplified Hamiltonian constitutes an easier tool for studying the statistical properties of this new orientational order, such as size dependence and phase transition related fluctuations. Indeed, as anticipated, the absence of the “dipole equilibration” process considerably reduces computing times. This reduction factor scales with the particle number. Table 1 shows the computing times
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
55
Fig. 6. hP4 i plotted as a function of the non-dimensional temperature τ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, ρ = 0.1, N = 1000) and the approximate Hamiltonian (open dots, ρ = 0.14, N = 8000). In the figure P4 = f (hP2 i) are portrayed for both cases (dashed line: approx. Hamilt.; continuous line: complete Hamilt.).
per macrostep obtained for systems of a different size by using the full and the approximated Hamiltonians. The acceleration is quite evident.
5. Conclusions This paper is devoted to a study, by Monte Carlo Metropolis computer simulations, of the collective behavior of a system of anisotropically polarizable particles under the action of an external electric field. We found new unexpected and previously unobserved behavior. Our results clearly show that, as an effect of the induced dipole – induced dipole interaction, the mean alignment of the system in the direction of the field decreases. This is an important finding because it may shed new light on the electro-optic response of the colloidal suspension of rod-like particles, which, at high particle density, was found to be unexpectedly small. In our simulations, the state of the spin system is defined by two independent and dimensionless quantities qualitatively representing its temperature and density. We performed simulations as a function of both parameters. We found that, at low temperature and large density, spin angular distribution shows a peak indicating that an important fraction of the spins aligns perpendicularly to the field. The onset of such an anomalous preferred configuration is captured well by comparing the values of the two moments hP2 i and hP4 i of the angular distribution. By studying the energy of the system we showed that this anomalous orientational configuration corresponds to a stable state of the system. Further investigation is required to properly characterize the statistical properties of this system and to investigate the phase transition which is probably associated with the onset of this anomalous state. In the present study, we used two different approaches: to calculate the energy of the system we used both the complete Hamiltonian (although limited to nearest neighbor interactions) and a simplified Hamiltonian obtained under the assumption of small density. From the standpoint of the simulation, the difference between the two models is large because the complete Hamiltonian uses an implicit definition of local fields which requires an
56
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
Fig. 7. Normalised mean energy as a function of the non-dimensional density ρ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, τ = 0.4, N = 1000) and the approximate Hamiltonian (open dots, τ = 0.1, N = 8000).
Fig. 8. Normalised mean energy as a function of non-dimensional temperature τ as obtained by a MM simulation performed using the complete Hamiltonian (full dots, ρ = 0.1, N = 1000) and the approximate Hamiltonian (open dots, ρ = 0.14, N = 8000).
G. Danese et al. / Computer Physics Communications 134 (2001) 47–57
57
Table 1 Comparison between the computation time (seconds per macrostep) needed by the complete and approximated equations for interaction energy Number of spins
Complete Hamiltonian (Eq. (4))
Approx. Hamiltonian (Eq. (14))
1000
49.04
0.023
512
13.22
0.011
256
0.747
0.005
64
0.198
0.004
iterative procedure to compute the energy, while in the simplified model the Hamiltonian is an explicit function of the spin orientation. By comparing the results of the two models it can be seen that the general behavior of the system is preserved even when the simplified approach, which drastically diminishes the computing times, is adopted.
Acknowledgements We thank Professor A. Maritan for having proposed the simplified Hamiltonian. We acknowledge useful discussions with M. Glaser, P. Pasini and C. Zannoni.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
M.E. Cates, J. Polyelectrolyte Phys. 2 (1992) 1109–1119. K. Kiyohara, K.E. Gubbins, A.Z. Panagiotopoulos, Mol. Phys. 94 (1998) 803–808. R. Berardi, S. Orlandi, C. Zannoni, Chem. Phys. Lett. 261 (1996) 357–362. N. Metropolis, A. Rosenbluth, M. Rosenbluth et al., J. Chem. Phys. 21 (1953) 1087–1092. M. Hloucha, U.K. Deiters, Fluid Phase Equil. 149 (1998) 41–56. U. Fabbri, C. Zannoni, Mol. Phys. 58 (1986) 763–788. D.E. Knuth, The Art of Computer Programming, Addison-Wesley, Reading, MA, 1981. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, 1987. P. Ewald, Ann. Phys. 64 (1921) 253–287. G. Danese, I. De Lotto, D. Dotti, F. Leporati, Comput. Physics Commun. 108 (1998) 211–217.