- Email: [email protected]
Grand canonical auxiliary field Monte Carlo: a new technique for simulating open systems at high density
Computer Physics Communications 157 (2004) 201–206 www.elsevier.com/locate/cpc
Grand canonical auxiliary field Monte Carlo: a new technique for simulating open systems at high density S.A. Baeurle a,b a Institut für Physikalische und Theoretische Chemie, D-93053 Regensburg, Germany b Max-Planck-Institut für Festkörperforschung, D-70569 Stuttgart, Germany
Received 1 July 2003
Abstract The computation of open many-particle systems at high densities is a major challenge since many decades due to the inherent limitations of grand canonical simulation methods based on particle exchange algorithms. In this paper we report on the statistical convergence behavior in the high density regime of a recently developed alternative called the grand canonical auxiliary field Monte Carlo method. We show on a common soft matter model widely used in polymer simulation that it possesses a more appropriate statistical behavior in the dense regime than the currently employed grand canonical Monte Carlo methods relying on particle exchange algorithms. 2003 Elsevier B.V. All rights reserved. PACS: 05.10.-a; 05.20.-y; 03.50.-z Keywords: Grand canonical auxiliary field Monte Carlo; Field-theoretic simulation method; Dense open systems; Numerical sign problem; Classical statistical mechanics
1. Introduction To perform simulations within the grand canonical ensemble is known to be very useful for the study of open systems, where matter and heat exchange between the system and its surroundings does occur. In such investigations one typically desires to extract information about the composition of the system as a function of the chemical potential and other thermodynamic properties. For instance, in a hydrogen E-mail address: [email protected] (S.A. Baeurle). URL: http://www-dick.chemie.uni-regensburg.de/group/ Stephan.Baeurle.html (S.A. Baeurle). 0010-4655/$ – see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.comphy.2003.11.001
storage study of carbon nanotubes one would like to know the amount of hydrogen molecules adsorbed as a function of the pressure of the hydrogen reservoir with which the nanotubes are in contact [1]. It is clear that one could alternatively describe such systems in the standard microcanonical or canonical ensemble. However, this approach is very inefficient due to a great deal of computational effort wasted on the particle bath and the system-bath interface, as well as on lengthy equilibration times, particularly at high density [2]. For the simulation within the grand canonical ensemble many techniques have been devised [3]. The most widely used method was developed by Norman et al. [4] and relies on the Metropolis Monte Carlo
202
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
(MC) algorithm to perform the particle displacements. To simulate the particle exchange between the physical system and the particle bath, it incorporates a supplementary particle creation/destruction step into the algorithm. Based on a similar strategy several alternatives have been developed, such as, e.g., the method of Rowley et al. [5]. The grand canonical Monte Carlo (GCMC) techniques have particularly been used to study phenomena, such as adsorption [2], capillary condensation [6] and vapor–liquid equilibrium [7]. They work best if the acceptance of trial moves by which the particles are created or destructed is not too low. However, this can represent a severe limitation when they are applied to condensed phases [3]. Through various investigations one has established that the particle exchange (creation/destruction) probabilities become continuously smaller and therefore the particle exchange algorithm increasingly more impracticable with growing density in the physical system. Another difficulty with these techniques is that the insertion or deletion of a particle is a local nonequilibrium event, and it takes some time for the system to relax. As a consequence recently inserted particles are more likely to be removed, and particles are also more likely to be inserted back into the cavity left by a recently removed particle [8]. All these reasons limit the applicability of the GCMC techniques, in case of atomic fluids, to about twice the critical density. To overcome these difficulties special strategies have been conceived to extend the GCMC techniques to somewhat higher densities, such as, e.g., the cavity-biased method of Mezei et al. [9]. New developments essentially based on this approach have recently provided some improved sampling efficiency [10]. However, there is an obvious inherent limitation of the particle-based approaches in their extensibility to the high density regime caused by their underlying particle exchange algorithm. Other methods make use of extended sampling schemes, in which particles are gradually inserted into the physical system, such as, e.g., the grand canonical ensemble molecular dynamics method of Cagin and Pettitt [11] or the method of Attard [12]. However, these methods are unphysical in nature, since they do not sample the true grand canonical distribution function. As a consequence the convergence to the correct thermodynamic averages can never be guaranteed and these methods
are found to provide wrong results in several important cases [13]. In this paper we investigate the usefulness of the auxiliary field functional integral approach in the regime of high particle density. In particular we demonstrate that the recently introduced classical grand canonical auxiliary field Monte Carlo (GCAFMC) method [14–16], which relies on a fieldtheoretical description of the many-particle system, shows a more convenient statistical behavior for computing open systems at high densities than the conventional grand canonical simulation methods relying on the particle description. To this end, we compare the GC-AFMC approach to the standard method of Norman et al. by calculating a model system widely used in polymer simulation. We emphasize that the GCAFMC method does not make use of a particle exchange algorithm and is also adaptable to related techniques, such as, e.g., the Gibbs ensemble method of Panagiotopoulos [17] and the method of Widom [18], which are plagued by similar problems. The paper is structured as follows. In Section 2 we review the GC-AFMC methodology and recall the basic aspects of the GCMC technique of Norman et al. Afterwards, in Section 3 we introduce the model system used in our simulations and, then, present and discuss the results in Section 4. Finally, we end the paper with a brief summary and the conclusions.
2. Grand canonical simulation methodologies To begin, let us first recall the basic steps of the classical auxiliary field methodology within the grand canonical ensemble [15]. In the primary step we perform an exact transformation of the conventional particle representation of the grand canonical partition function ∞ eβµN Ξ (µ, V , T ) = dR N dP N h3N N! N=0 × exp −βH (R N , P N ) (1) into the basic field representation (BFR) Ξ (µ, V , T ) = Dσ exp −S[σ ]
(2)
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
with the effective action 1 ∗ −1 S[σ ] = σ (G)Φ (G)σ (G) 2β G χ − e−iσ (r ) dr , V
(3)
where σ (r ) represents function and the auxiliary field = 0) dσ R (G) dσ I (G) the integraDσ ∝ dσ (G G>0 tion measure depending on the set σ of field coef σ I (G), . . .}, restricted to half of the ficients {σ R (G), reciprocal lattice. The constant χ is a parameter which depends on the chemical potential µ through χ = V /λ3B exp[βµ + 12 βΦ(0)]. For the transformation we use the Hubbard–Stratonovich relation given by [15]
2 ∞ a 2 y 1 exp − x = exp − − ixy dy, (4) 2 2πa 2a −∞
where the constant a > 0. From our previous work [16] we know that the BFR of the grand canonical partition function does not permit the crude application of the standard multidimensional integration techniques, such as, e.g., the MC algorithm. The problem is related to the complex nature of the weight function exp[−S[σ ]], which causes a bad statistical convergence of the ensemble averages in the low temperature and/or large system size regime. Thus, to make the methodology amenable for computation, the contourdistortion technique based on Cauchy’s integral theorem has to be applied. In essence, the technique consists in reformulating the partition functional integral by performing a suitable contour-shift into the complex plane according to the prescription σ (r ) → σ (r ) + iψ(r ),
203
(5)
where ψ(r ) represents the shifting function. The goal of this procedure is to distort the original integration path in such way that one captures as many critical points as possible, providing a relevant contribution to the functional integral entity. In Fig. 1 we show a deformed path of the BFR, which crosses its homogeneous saddle point (SP) lying on the imaginary axis of the complex σ (r )-plane. This SP characterizes a mean field (MF) approximation to the liquid phase of the model system under consideration. From our previous investigations [16] we know that such a MF shift is not necessarily the most optimal shift with respect to
Fig. 1. Complex plane of the field σ (r ) at a fixed value of r. The original path of integration of the grand canonical partition function is along the real axis.
statistical convergence. In fact we could demonstrate that the method of Gaussian equivalent representation (GER) provides an improved shift, which causes a significant acceleration of the statistical convergence in the low temperature and/or large system regime and remains effective far beyond mean field approximation. For computational purposes the complex weight formulation of the grand canonical partition function has further to be recasted in a real weight formulation [15]. As a result one obtains a real non-positive definite weight function which cannot directly be used for sampling, since conventional numerical integration techniques, such as, e.g., MC, work only with real positive ones. The common approach in such a case is to employ a reweighting procedure [19], which consists in factorizing the original weight function into a real and positive part, called the reference weight function ρ ref [σ ], and a remainder, included in the estimator. The ensemble average of a property A is then calculated by evaluating the expression A =
A[σ ] sgn[ρ[σ ]]ref , sgn[ρ[σ ]]ref
(6)
where the brackets · · ·ref denote the averaging with respect to the absolute value of the real weight function ρ[σ ] ∝ cos(Im[S[σ ]]) exp[− Re[S[σ ]]] and A[σ ] the corresponding field estimator. A reliable indicator for the statistical convergence of GC-AFMC is the average in the numerator of (6), commonly known as the average sign. It can quite generally be said that, if this quantity amounts to one, ρ[σ ]
204
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
is positive definite and the statistical convergence is optimal, while approaching zero the convergence deteriorates. This phenomenon is commonly known as the numerical sign problem [16]. In order to assess the competitiveness of GCAFMC with regard to conventional grand canonical simulation methods, we make use of the GCMC technique of Norman et al., which is the most commonly employed GCMC technique in current research [3]. From the calculations we extract the particle creation/destruction probabilities Pcreate and Pdestruct , which are a measure for the statistical convergence of the method and therefore the analog to the average sign of GC-AFMC.
3. Model system
Fig. 2. Average density as a function of the excess chemical potential obtained from GER-GC-AFMC simulations of a system of Gauss-core particles in comparison to GCMC results. All statistical errors are smaller than symbol size.
For our investigations we chose a system of particles interacting through a purely repulsive Gaussian potential, i.e. Φ(r) = Φ(0) exp −(r/ l)2 . (7) The so-called Gauss-core (G) potential was recently found to be very useful in the description of a wide variety of soft matter systems [20], such as, e.g., nonintersecting polymers over the whole range of particle densities [21]. Moreover, it is well known from exact quantum ab initio calculations that a repulsive potential of exponential type generally reflects more realistically the repulsive part of the interaction in dense atomic systems than the standard Lennard-Jones potential [22,23]. Finally, this choice has also been motivated by the assumption that the convergence difficulties are mainly caused by the repulsive contribution to the interaction.
4. Results and discussion To study the convergence behavior of GC-AFMC with increasing density, we have performed several GER-GC-AFMC simulations of a system of Gausscore particles increasing the chemical potential successively. In these runs we chose a box size of VG∗ = 20.0 and a temperature of TG∗ = 0.1. Moreover, we ∗ = 60.0, which ensures set the field cutoff to Ecut/G a correct representation of the field function over the
Fig. 3. Average sign of GER-GC-AFMC in comparison to the particle exchange probabilities of GCMC as a function of the average density obtained for a system of Gauss-core particles. All statistical errors are smaller than symbol size.
whole density range. For completeness we specify that in our investigations we used the system of units characteristic for the model [24]. In Fig. 2 we have visualized the results obtained for the average density as a function of the excess chemical potential. For comparison, we show the ones yielded from GCMC simulations performed at the same temperature and chemical potentials, employing a box size of VG∗ = 512.0 [25]. From the plot, we conclude that the densities of both methods are in good agreement with each other for all chemical potentials considered. Next, we show in Fig. 3 the average sign as a function of the average
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
density resulting from the GER-GC-AFMC simulations mentioned previously in comparison to the particle creation/destruction probabilities obtained from GCMC. From these results we conclude that the average sign of GC-AFMC has a more favorable behavior than the probabilities of GCMC with growing particle density. The curve of the average sign possesses ∗ ≈ 0.25 a minimum of sgn[ρ[σ ]]ref ≈ 0.1 at ρG and tends to one in the low and high density limit. In contrast, we see that the GCMC probabilities continuously decrease to zero as the density of the system increases. Next, we show in Fig. 4 the pair repul∗ (r ∗ ) at the typical particle separation sion energy ΦG G ∗−1/3 ∗ [26] and its derivative with respect to the aG = ρG density, which characterize, respectively, the importance of the repulsion and the interaction strength between the particles. From these curves we deduce that approaching the low and high density limits the forces acting between the particles become smaller, while the pair repulsion energy grows continuously with increasing density. We can explain this phenomenon by the fact that in the low density regime the particles are far apart and do not strongly feel each other, whereas in the high density regime they are close together and lie on similar potential energy levels. Approaching the density of the minimum average sign, the repulsive forces between the particles increase, which leads to a stronger oscillation of ρ[σ ] and therefore results in an enhanced cancellation of the contributions to the average sign. A further indication confirming this analysis is that, at this density, the Gauss-core model also possesses a maximum in the melting temperature [24]. Thus, we deduce from the above discussion that the average sign reflects the interaction strength between the particles at the given external conditions. In contrast, we attribute the decrease of the GCMC probabilities to the fact that the energy difference between Nth- and (N + 1)th-particle configuration continuously increases with growing density, as we can deduce from the behavior of the pair repulsion energy given in Fig. 4. In the high density limit of the Gausscore model this energy difference goes over in a constant, causing that the probabilities converge to an asymptotic limit. Thus, we see that the problems are not related. We note that the Gauss-core model is also a useful example to demonstrate the importance of selecting the ideal contour of integration, because it possesses a
205
Fig. 4. Pair repulsion energy and interaction strength for the Gauss-core-, Yukawa- and Coulomb-potential as a function of the density. All quantities are expressed in the reduced units, which are natural for the respective model.
single dominant SP at this temperature over the whole density range. This SP is characterized by the spaceindependent MF solution [15] = 0), σMF (r ) = iψMF (G
(8)
= 0) represents the shift obtained by where ψMF (G solving the MF equation −
= 0) ψMF (G = χeψMF (G=0) . βΦ(G = 0)
(9)
Quite generally we can say that GC-AFMC is useful in the parameter range, where one can capture the major contributions of the relevant critical points of the partition function integral, like in the liquid state of the Gauss-core model considered here. Fortunately, it appears that a wide variety of soft matter models have similar mathematical properties such as, e.g., systems of polymers [20,21], dendrimers [20] as well as nonionic and charged colloids [24,26]. However, we emphasize that the usefulness of GC-AFMC is not limited to soft-core interactions. It can be extended to models with a hard-core at zero particle distance with the restriction that the potential must be invertible [27], such as, e.g., in case of the Coulomb potential v(r) = a/r or the Yukawa potential v(r) = (a/r) exp[−κr] where a, κ > 0. In Fig. 4 we have included the pair repulsion energy and the interaction strength of these two interaction models. We can easily see that, similarly to the Gauss-core potential, they are both characterized by an increase in the pair repulsion energy and a decrease
206
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
of the interaction strength in the high density range. In case of the Yukawa potential it is also worth mentioning that it possesses several SPs, which are all located on the imaginary axis [27]. Each of these SPs represents a MF approximation to a thermodynamic phase of the model. More specifically one can attribute a homogeneous SP to its fluid phase with much the same character as the SP of the Gauss-core model discussed previously [26]. In addition, the Yukawa model also possesses a number of inhomogeneous SPs, which represent MF approximations to various crystal phases and can be determined with a standard self-consistent field theory (SCFT) technique [27]. In three dimensions it possesses, e.g., SPs with body centered cubic (BCC) and face centered cubic (FCC) symmetries. In Fig. 1 we have schematically illustrated in the complex σ (r )-plane the field configurations corresponding to the three stable thermodynamic phases of the Yukawa model. If one orders these SPs according to the real parts of S[σ ], i.e. Re[S[σMF/1 ]] < Re[S[σMF/2 ]] < Re[S[σMF/3 ]], . . . , then one clearly recognizes that there exists a dominant SP for each thermodynamic phase, while the others are less important at the external conditions under consideration. Thus, we can safely predict that applying the GER approach for each phase would permit us to calculate the desired thermodynamics properties far beyond the accuracy of the corresponding MF approximations. 5. Summary and conclusions In conclusion, we have demonstrated in this paper that the grand canonical auxiliary field Monte Carlo method possesses the appropriate convergence behavior for investigating open soft matter systems in the high density regime. On the example of the Gauss-core model we have shown that the numerical sign problem of GC-AFMC disappears, while the particle exchange probabilities of GCMC decrease continuously with growing density. Finally, we have discussed the extensibility of GC-AFMC to hard-core type interactions and outlined the conditions under which one can expect an optimal convergence behavior. Acknowledgements We thank Prof. Michele Parrinello and Dr. Roman Martonak for offering stimulating discussions and en-
couragements. Moreover, we gratefully acknowledge the support of Prof. Bernhard Dick and Dr. Alkwin Slenczka.
References [1] R.H. Baughman, A.A. Zakhidov, W.A. de Heer, Science 297 (2002) 787; F.L. Darkrim, P. Malbrunot, G.P. Tartaglia, Int. J. Hydrog. Energy 27 (2002) 193; C. Liu, Y.Y. Fan, M. Liu, H.T. Cong, H.M. Cheng, M.S. Dresselhaus, Science 286 (1999) 1127. [2] D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, San Diego, 1996. [3] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1996. [4] G.E. Norman, V.S. Filinov, High Temp. (USSR) 7 (1969) 216. [5] L.A. Rowley, D. Nicholson, N.G. Parsonage, J. Comput. Phys. 17 (1975) 401. [6] S. Sokołowski, Mol. Phys. 75 (1992) 999. [7] D. Boda, J. Liszi, I. Szalai, Chem. Phys. Lett. 256 (1996) 474; D. Boda, K. Chan, I. Szalai, Mol. Phys. 92 (1997) 1067. [8] J.C. Shelley, G.N. Patey, J. Chem. Phys. 102 (1995) 7656. [9] M. Mezei, Mol. Phys. 40 (1980) 901. [10] R.M. Shroll, D.E. Smith, J. Chem. Phys. 110 (1999) 8295, and references therein. [11] T. Cagin, B.M. Pettitt, Mol. Sim. 6 (1991) 5; Mol. Phys. 72 (1991) 169; J. Ji, T. Cagin, B.M. Pettitt, J. Chem. Phys. 96 (1992) 1333. [12] P. Attard, J. Chem. Phys. 107 (1997) 3230. [13] S. Weerasinghe, B.M. Pettitt, Mol. Phys. 82 (1994) 897. [14] S.A. Baeurle, Ph.D. Thesis, University of Stuttgart, Stuttgart, 2000. [15] S.A. Baeurle, R. Martonak, M. Parrinello, J. Chem. Phys. 117 (2002) 3027. [16] S.A. Baeurle, Phys. Rev. Lett. 89 (2002) 080602. [17] A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813. [18] B. Widom, J. Chem. Phys. 39 (1963) 2808. [19] S.A. Baeurle, J. Comput. Phys. 184 (2003) 540. [20] C.N. Likos, N. Hoffmann, H. Löwen, A.A. Louis, J. Phys.: Condens. Matter 14 (2002) 7681. [21] A.A. Louis, P.G. Bolhuis, J.P. Hansen, E.J. Meijer, Phys. Rev. Lett. 85 (2000) 2522. [22] G. Sutmann, Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, in: NIC Series, 2002. [23] F.H. Stillinger, T.A. Weber, Phys. Rev. B 22 (1980) 3790. [24] F.H. Stillinger, D.K. Stillinger, Physica A 244 (1997) 358. [25] It is worth noting that we chose the volume large enough to permit the simulation over the whole density range with the same box size. [26] M. Dijkstra, R. van Roij, J. Phys.: Condens. Matter 10 (1998) 1219. [27] G.H. Fredrickson, V. Ganesan, F. Drolet, Macromolecules 35 (2002) 16.
Grand canonical auxiliary field Monte Carlo: a new technique for simulating open systems at high density S.A. Baeurle a,b a Institut für Physikalische und Theoretische Chemie, D-93053 Regensburg, Germany b Max-Planck-Institut für Festkörperforschung, D-70569 Stuttgart, Germany
Received 1 July 2003
Abstract The computation of open many-particle systems at high densities is a major challenge since many decades due to the inherent limitations of grand canonical simulation methods based on particle exchange algorithms. In this paper we report on the statistical convergence behavior in the high density regime of a recently developed alternative called the grand canonical auxiliary field Monte Carlo method. We show on a common soft matter model widely used in polymer simulation that it possesses a more appropriate statistical behavior in the dense regime than the currently employed grand canonical Monte Carlo methods relying on particle exchange algorithms. 2003 Elsevier B.V. All rights reserved. PACS: 05.10.-a; 05.20.-y; 03.50.-z Keywords: Grand canonical auxiliary field Monte Carlo; Field-theoretic simulation method; Dense open systems; Numerical sign problem; Classical statistical mechanics
1. Introduction To perform simulations within the grand canonical ensemble is known to be very useful for the study of open systems, where matter and heat exchange between the system and its surroundings does occur. In such investigations one typically desires to extract information about the composition of the system as a function of the chemical potential and other thermodynamic properties. For instance, in a hydrogen E-mail address: [email protected] (S.A. Baeurle). URL: http://www-dick.chemie.uni-regensburg.de/group/ Stephan.Baeurle.html (S.A. Baeurle). 0010-4655/$ – see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.comphy.2003.11.001
storage study of carbon nanotubes one would like to know the amount of hydrogen molecules adsorbed as a function of the pressure of the hydrogen reservoir with which the nanotubes are in contact [1]. It is clear that one could alternatively describe such systems in the standard microcanonical or canonical ensemble. However, this approach is very inefficient due to a great deal of computational effort wasted on the particle bath and the system-bath interface, as well as on lengthy equilibration times, particularly at high density [2]. For the simulation within the grand canonical ensemble many techniques have been devised [3]. The most widely used method was developed by Norman et al. [4] and relies on the Metropolis Monte Carlo
202
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
(MC) algorithm to perform the particle displacements. To simulate the particle exchange between the physical system and the particle bath, it incorporates a supplementary particle creation/destruction step into the algorithm. Based on a similar strategy several alternatives have been developed, such as, e.g., the method of Rowley et al. [5]. The grand canonical Monte Carlo (GCMC) techniques have particularly been used to study phenomena, such as adsorption [2], capillary condensation [6] and vapor–liquid equilibrium [7]. They work best if the acceptance of trial moves by which the particles are created or destructed is not too low. However, this can represent a severe limitation when they are applied to condensed phases [3]. Through various investigations one has established that the particle exchange (creation/destruction) probabilities become continuously smaller and therefore the particle exchange algorithm increasingly more impracticable with growing density in the physical system. Another difficulty with these techniques is that the insertion or deletion of a particle is a local nonequilibrium event, and it takes some time for the system to relax. As a consequence recently inserted particles are more likely to be removed, and particles are also more likely to be inserted back into the cavity left by a recently removed particle [8]. All these reasons limit the applicability of the GCMC techniques, in case of atomic fluids, to about twice the critical density. To overcome these difficulties special strategies have been conceived to extend the GCMC techniques to somewhat higher densities, such as, e.g., the cavity-biased method of Mezei et al. [9]. New developments essentially based on this approach have recently provided some improved sampling efficiency [10]. However, there is an obvious inherent limitation of the particle-based approaches in their extensibility to the high density regime caused by their underlying particle exchange algorithm. Other methods make use of extended sampling schemes, in which particles are gradually inserted into the physical system, such as, e.g., the grand canonical ensemble molecular dynamics method of Cagin and Pettitt [11] or the method of Attard [12]. However, these methods are unphysical in nature, since they do not sample the true grand canonical distribution function. As a consequence the convergence to the correct thermodynamic averages can never be guaranteed and these methods
are found to provide wrong results in several important cases [13]. In this paper we investigate the usefulness of the auxiliary field functional integral approach in the regime of high particle density. In particular we demonstrate that the recently introduced classical grand canonical auxiliary field Monte Carlo (GCAFMC) method [14–16], which relies on a fieldtheoretical description of the many-particle system, shows a more convenient statistical behavior for computing open systems at high densities than the conventional grand canonical simulation methods relying on the particle description. To this end, we compare the GC-AFMC approach to the standard method of Norman et al. by calculating a model system widely used in polymer simulation. We emphasize that the GCAFMC method does not make use of a particle exchange algorithm and is also adaptable to related techniques, such as, e.g., the Gibbs ensemble method of Panagiotopoulos [17] and the method of Widom [18], which are plagued by similar problems. The paper is structured as follows. In Section 2 we review the GC-AFMC methodology and recall the basic aspects of the GCMC technique of Norman et al. Afterwards, in Section 3 we introduce the model system used in our simulations and, then, present and discuss the results in Section 4. Finally, we end the paper with a brief summary and the conclusions.
2. Grand canonical simulation methodologies To begin, let us first recall the basic steps of the classical auxiliary field methodology within the grand canonical ensemble [15]. In the primary step we perform an exact transformation of the conventional particle representation of the grand canonical partition function ∞ eβµN Ξ (µ, V , T ) = dR N dP N h3N N! N=0 × exp −βH (R N , P N ) (1) into the basic field representation (BFR) Ξ (µ, V , T ) = Dσ exp −S[σ ]
(2)
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
with the effective action 1 ∗ −1 S[σ ] = σ (G)Φ (G)σ (G) 2β G χ − e−iσ (r ) dr , V
(3)
where σ (r ) represents function and the auxiliary field = 0) dσ R (G) dσ I (G) the integraDσ ∝ dσ (G G>0 tion measure depending on the set σ of field coef σ I (G), . . .}, restricted to half of the ficients {σ R (G), reciprocal lattice. The constant χ is a parameter which depends on the chemical potential µ through χ = V /λ3B exp[βµ + 12 βΦ(0)]. For the transformation we use the Hubbard–Stratonovich relation given by [15]
2 ∞ a 2 y 1 exp − x = exp − − ixy dy, (4) 2 2πa 2a −∞
where the constant a > 0. From our previous work [16] we know that the BFR of the grand canonical partition function does not permit the crude application of the standard multidimensional integration techniques, such as, e.g., the MC algorithm. The problem is related to the complex nature of the weight function exp[−S[σ ]], which causes a bad statistical convergence of the ensemble averages in the low temperature and/or large system size regime. Thus, to make the methodology amenable for computation, the contourdistortion technique based on Cauchy’s integral theorem has to be applied. In essence, the technique consists in reformulating the partition functional integral by performing a suitable contour-shift into the complex plane according to the prescription σ (r ) → σ (r ) + iψ(r ),
203
(5)
where ψ(r ) represents the shifting function. The goal of this procedure is to distort the original integration path in such way that one captures as many critical points as possible, providing a relevant contribution to the functional integral entity. In Fig. 1 we show a deformed path of the BFR, which crosses its homogeneous saddle point (SP) lying on the imaginary axis of the complex σ (r )-plane. This SP characterizes a mean field (MF) approximation to the liquid phase of the model system under consideration. From our previous investigations [16] we know that such a MF shift is not necessarily the most optimal shift with respect to
Fig. 1. Complex plane of the field σ (r ) at a fixed value of r. The original path of integration of the grand canonical partition function is along the real axis.
statistical convergence. In fact we could demonstrate that the method of Gaussian equivalent representation (GER) provides an improved shift, which causes a significant acceleration of the statistical convergence in the low temperature and/or large system regime and remains effective far beyond mean field approximation. For computational purposes the complex weight formulation of the grand canonical partition function has further to be recasted in a real weight formulation [15]. As a result one obtains a real non-positive definite weight function which cannot directly be used for sampling, since conventional numerical integration techniques, such as, e.g., MC, work only with real positive ones. The common approach in such a case is to employ a reweighting procedure [19], which consists in factorizing the original weight function into a real and positive part, called the reference weight function ρ ref [σ ], and a remainder, included in the estimator. The ensemble average of a property A is then calculated by evaluating the expression A =
A[σ ] sgn[ρ[σ ]]ref , sgn[ρ[σ ]]ref
(6)
where the brackets · · ·ref denote the averaging with respect to the absolute value of the real weight function ρ[σ ] ∝ cos(Im[S[σ ]]) exp[− Re[S[σ ]]] and A[σ ] the corresponding field estimator. A reliable indicator for the statistical convergence of GC-AFMC is the average in the numerator of (6), commonly known as the average sign. It can quite generally be said that, if this quantity amounts to one, ρ[σ ]
204
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
is positive definite and the statistical convergence is optimal, while approaching zero the convergence deteriorates. This phenomenon is commonly known as the numerical sign problem [16]. In order to assess the competitiveness of GCAFMC with regard to conventional grand canonical simulation methods, we make use of the GCMC technique of Norman et al., which is the most commonly employed GCMC technique in current research [3]. From the calculations we extract the particle creation/destruction probabilities Pcreate and Pdestruct , which are a measure for the statistical convergence of the method and therefore the analog to the average sign of GC-AFMC.
3. Model system
Fig. 2. Average density as a function of the excess chemical potential obtained from GER-GC-AFMC simulations of a system of Gauss-core particles in comparison to GCMC results. All statistical errors are smaller than symbol size.
For our investigations we chose a system of particles interacting through a purely repulsive Gaussian potential, i.e. Φ(r) = Φ(0) exp −(r/ l)2 . (7) The so-called Gauss-core (G) potential was recently found to be very useful in the description of a wide variety of soft matter systems [20], such as, e.g., nonintersecting polymers over the whole range of particle densities [21]. Moreover, it is well known from exact quantum ab initio calculations that a repulsive potential of exponential type generally reflects more realistically the repulsive part of the interaction in dense atomic systems than the standard Lennard-Jones potential [22,23]. Finally, this choice has also been motivated by the assumption that the convergence difficulties are mainly caused by the repulsive contribution to the interaction.
4. Results and discussion To study the convergence behavior of GC-AFMC with increasing density, we have performed several GER-GC-AFMC simulations of a system of Gausscore particles increasing the chemical potential successively. In these runs we chose a box size of VG∗ = 20.0 and a temperature of TG∗ = 0.1. Moreover, we ∗ = 60.0, which ensures set the field cutoff to Ecut/G a correct representation of the field function over the
Fig. 3. Average sign of GER-GC-AFMC in comparison to the particle exchange probabilities of GCMC as a function of the average density obtained for a system of Gauss-core particles. All statistical errors are smaller than symbol size.
whole density range. For completeness we specify that in our investigations we used the system of units characteristic for the model [24]. In Fig. 2 we have visualized the results obtained for the average density as a function of the excess chemical potential. For comparison, we show the ones yielded from GCMC simulations performed at the same temperature and chemical potentials, employing a box size of VG∗ = 512.0 [25]. From the plot, we conclude that the densities of both methods are in good agreement with each other for all chemical potentials considered. Next, we show in Fig. 3 the average sign as a function of the average
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
density resulting from the GER-GC-AFMC simulations mentioned previously in comparison to the particle creation/destruction probabilities obtained from GCMC. From these results we conclude that the average sign of GC-AFMC has a more favorable behavior than the probabilities of GCMC with growing particle density. The curve of the average sign possesses ∗ ≈ 0.25 a minimum of sgn[ρ[σ ]]ref ≈ 0.1 at ρG and tends to one in the low and high density limit. In contrast, we see that the GCMC probabilities continuously decrease to zero as the density of the system increases. Next, we show in Fig. 4 the pair repul∗ (r ∗ ) at the typical particle separation sion energy ΦG G ∗−1/3 ∗ [26] and its derivative with respect to the aG = ρG density, which characterize, respectively, the importance of the repulsion and the interaction strength between the particles. From these curves we deduce that approaching the low and high density limits the forces acting between the particles become smaller, while the pair repulsion energy grows continuously with increasing density. We can explain this phenomenon by the fact that in the low density regime the particles are far apart and do not strongly feel each other, whereas in the high density regime they are close together and lie on similar potential energy levels. Approaching the density of the minimum average sign, the repulsive forces between the particles increase, which leads to a stronger oscillation of ρ[σ ] and therefore results in an enhanced cancellation of the contributions to the average sign. A further indication confirming this analysis is that, at this density, the Gauss-core model also possesses a maximum in the melting temperature [24]. Thus, we deduce from the above discussion that the average sign reflects the interaction strength between the particles at the given external conditions. In contrast, we attribute the decrease of the GCMC probabilities to the fact that the energy difference between Nth- and (N + 1)th-particle configuration continuously increases with growing density, as we can deduce from the behavior of the pair repulsion energy given in Fig. 4. In the high density limit of the Gausscore model this energy difference goes over in a constant, causing that the probabilities converge to an asymptotic limit. Thus, we see that the problems are not related. We note that the Gauss-core model is also a useful example to demonstrate the importance of selecting the ideal contour of integration, because it possesses a
205
Fig. 4. Pair repulsion energy and interaction strength for the Gauss-core-, Yukawa- and Coulomb-potential as a function of the density. All quantities are expressed in the reduced units, which are natural for the respective model.
single dominant SP at this temperature over the whole density range. This SP is characterized by the spaceindependent MF solution [15] = 0), σMF (r ) = iψMF (G
(8)
= 0) represents the shift obtained by where ψMF (G solving the MF equation −
= 0) ψMF (G = χeψMF (G=0) . βΦ(G = 0)
(9)
Quite generally we can say that GC-AFMC is useful in the parameter range, where one can capture the major contributions of the relevant critical points of the partition function integral, like in the liquid state of the Gauss-core model considered here. Fortunately, it appears that a wide variety of soft matter models have similar mathematical properties such as, e.g., systems of polymers [20,21], dendrimers [20] as well as nonionic and charged colloids [24,26]. However, we emphasize that the usefulness of GC-AFMC is not limited to soft-core interactions. It can be extended to models with a hard-core at zero particle distance with the restriction that the potential must be invertible [27], such as, e.g., in case of the Coulomb potential v(r) = a/r or the Yukawa potential v(r) = (a/r) exp[−κr] where a, κ > 0. In Fig. 4 we have included the pair repulsion energy and the interaction strength of these two interaction models. We can easily see that, similarly to the Gauss-core potential, they are both characterized by an increase in the pair repulsion energy and a decrease
206
S.A. Baeurle / Computer Physics Communications 157 (2004) 201–206
of the interaction strength in the high density range. In case of the Yukawa potential it is also worth mentioning that it possesses several SPs, which are all located on the imaginary axis [27]. Each of these SPs represents a MF approximation to a thermodynamic phase of the model. More specifically one can attribute a homogeneous SP to its fluid phase with much the same character as the SP of the Gauss-core model discussed previously [26]. In addition, the Yukawa model also possesses a number of inhomogeneous SPs, which represent MF approximations to various crystal phases and can be determined with a standard self-consistent field theory (SCFT) technique [27]. In three dimensions it possesses, e.g., SPs with body centered cubic (BCC) and face centered cubic (FCC) symmetries. In Fig. 1 we have schematically illustrated in the complex σ (r )-plane the field configurations corresponding to the three stable thermodynamic phases of the Yukawa model. If one orders these SPs according to the real parts of S[σ ], i.e. Re[S[σMF/1 ]] < Re[S[σMF/2 ]] < Re[S[σMF/3 ]], . . . , then one clearly recognizes that there exists a dominant SP for each thermodynamic phase, while the others are less important at the external conditions under consideration. Thus, we can safely predict that applying the GER approach for each phase would permit us to calculate the desired thermodynamics properties far beyond the accuracy of the corresponding MF approximations. 5. Summary and conclusions In conclusion, we have demonstrated in this paper that the grand canonical auxiliary field Monte Carlo method possesses the appropriate convergence behavior for investigating open soft matter systems in the high density regime. On the example of the Gauss-core model we have shown that the numerical sign problem of GC-AFMC disappears, while the particle exchange probabilities of GCMC decrease continuously with growing density. Finally, we have discussed the extensibility of GC-AFMC to hard-core type interactions and outlined the conditions under which one can expect an optimal convergence behavior. Acknowledgements We thank Prof. Michele Parrinello and Dr. Roman Martonak for offering stimulating discussions and en-
couragements. Moreover, we gratefully acknowledge the support of Prof. Bernhard Dick and Dr. Alkwin Slenczka.
References [1] R.H. Baughman, A.A. Zakhidov, W.A. de Heer, Science 297 (2002) 787; F.L. Darkrim, P. Malbrunot, G.P. Tartaglia, Int. J. Hydrog. Energy 27 (2002) 193; C. Liu, Y.Y. Fan, M. Liu, H.T. Cong, H.M. Cheng, M.S. Dresselhaus, Science 286 (1999) 1127. [2] D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, San Diego, 1996. [3] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1996. [4] G.E. Norman, V.S. Filinov, High Temp. (USSR) 7 (1969) 216. [5] L.A. Rowley, D. Nicholson, N.G. Parsonage, J. Comput. Phys. 17 (1975) 401. [6] S. Sokołowski, Mol. Phys. 75 (1992) 999. [7] D. Boda, J. Liszi, I. Szalai, Chem. Phys. Lett. 256 (1996) 474; D. Boda, K. Chan, I. Szalai, Mol. Phys. 92 (1997) 1067. [8] J.C. Shelley, G.N. Patey, J. Chem. Phys. 102 (1995) 7656. [9] M. Mezei, Mol. Phys. 40 (1980) 901. [10] R.M. Shroll, D.E. Smith, J. Chem. Phys. 110 (1999) 8295, and references therein. [11] T. Cagin, B.M. Pettitt, Mol. Sim. 6 (1991) 5; Mol. Phys. 72 (1991) 169; J. Ji, T. Cagin, B.M. Pettitt, J. Chem. Phys. 96 (1992) 1333. [12] P. Attard, J. Chem. Phys. 107 (1997) 3230. [13] S. Weerasinghe, B.M. Pettitt, Mol. Phys. 82 (1994) 897. [14] S.A. Baeurle, Ph.D. Thesis, University of Stuttgart, Stuttgart, 2000. [15] S.A. Baeurle, R. Martonak, M. Parrinello, J. Chem. Phys. 117 (2002) 3027. [16] S.A. Baeurle, Phys. Rev. Lett. 89 (2002) 080602. [17] A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813. [18] B. Widom, J. Chem. Phys. 39 (1963) 2808. [19] S.A. Baeurle, J. Comput. Phys. 184 (2003) 540. [20] C.N. Likos, N. Hoffmann, H. Löwen, A.A. Louis, J. Phys.: Condens. Matter 14 (2002) 7681. [21] A.A. Louis, P.G. Bolhuis, J.P. Hansen, E.J. Meijer, Phys. Rev. Lett. 85 (2000) 2522. [22] G. Sutmann, Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, in: NIC Series, 2002. [23] F.H. Stillinger, T.A. Weber, Phys. Rev. B 22 (1980) 3790. [24] F.H. Stillinger, D.K. Stillinger, Physica A 244 (1997) 358. [25] It is worth noting that we chose the volume large enough to permit the simulation over the whole density range with the same box size. [26] M. Dijkstra, R. van Roij, J. Phys.: Condens. Matter 10 (1998) 1219. [27] G.H. Fredrickson, V. Ganesan, F. Drolet, Macromolecules 35 (2002) 16.