J. Phys. Chem. B 2001, 105, 11275-11282
11275
Improving the Efficiency of the Aggregation-Volume-Bias Monte Carlo Algorithm Bin Chen† and J. Ilja Siepmann* Department of Chemistry and Department of Chemical Engineering and Materials Science, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455-0431 ReceiVed: June 12, 2001; In Final Form: September 14, 2001
The aggregation-volume-bias Monte Carlo (AVBMC) algorithm is reanalyzed, and on the basis of this analysis, two extensions of the AVBMC algorithm with improved sampling efficiency for super-strongly associating fluids are presented. The new versions of the AVBMC algorithm are based on the principle of super-detailed balance and retain the simplicity, generality, and robustness of the original AVBMC algorithm. The performances of the various versions of the AVBMC algorithm are compared via applications to the simple ideal-association model of van Roij and to the superheated vapor phase of hydrogen fluoride.
1. Introduction Associating fluids are an important class of systems wellknown for their nonideal fluid-phase behavior due to the formation of aggregates even in low-density gas phases and lowconcentration solutions.1 The “physical” nature of this type of association can be traced to particular intermolecular attractions with deep (a few tens of kJ/mol or thousands of Kelvin in thermal units) but narrow wells usually caused by the formation of intermolecular hydrogen bonds. Water, alcohols, organic acids, and other molecules which are capable of hydrogen bonding are typical examples of associating fluids. Molecular simulations of associating fluids are a challenging task because the presence of “bonded” configurations (or aggregates) leads to bottlenecks or traps in phase space. These configurations have a very large Boltzmann weight associated with the large favorable energies of aggregate formation. However, they represent only a small fraction of the total phase space compared to the large amount of nonbonded configurations. Thus, an efficient sampling algorithm for associating fluids has to be able to find these bonded configurations but not get trapped by them. Unfortunately, conventional Monte Carlo simulations employing the standard Metropolis algorithm2 might either miss forming bonded configurations (because of the extremely low probability to locate them using a symmetric underlying transition matrix) or get trapped once a bonded configuration is found (because of the large energy penalty and correspondingly low acceptance probability of moves that lead to aggregate dissociation). Over the last several years, various biased Monte Carlo algorithms3-9 have been proposed to overcome the shortcomings of the conventional Metropolis scheme. One general feature of these algorithms is to preferentially (more frequently) attempt moves that lead to the formation of bonded configurations, which in turn enhances the acceptance probability for moves that lead to the destruction of bonded configurations after removing the bias. These algorithms differ to a large extent in complexity and generality. The complexities and restrictions of * To whom correspondence should be addressed. E-mail: siepmann@ chem.umn.edu. † Present address: Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6323.
these biased algorithms have been recently discussed in detail by Kofke and co-workers7,8 and by us.9 It has been demonstrated that the aggregation-volume-bias Monte Carlo (AVBMC) algorithm9 and the unbonding and bonding (UB) algorithm proposed by Wierzchowski and Kofke8 are much simpler and more general than other biased schemes. Wierzchowski and Kofke have shown that the UB algorithm samples the phase space for super-strongly associating fluids (where the majority of molecules is bound to at least one other molecule) more efficiently than the original AVBMC algorithm.8 As suggested by its name, the UB algorithm makes use of two trial moves that explicitly lead to either unbonding or bonding of molecules, thereby sampling aggregate formation and destruction. The UB algorithm considers all possible ways that a biased move could arrive at a specific outcome to calculate the acceptance probability, for which it needs to determine the number of bound molecules for the trial position.8 Furthermore, a list of the bound molecules (associated with at least one other molecule) and the total number of bound molecules need to be tracked throughout an UB simulation. In contrast, the original AVBMC algorithm requires less information to determine the acceptance of a trial move and relies on the concept of “super-detailed balance” (introduced first for configurational-bias Monte Carlo moves) where only a subset of ways to move from one configuration to another is considered and paired with a corresponding counterpart for the reverse move. In the next section, the aggregation-volume-bias Monte Carlo algorithm will be reanalyzed. This analysis leads to two novel versions of the AVBMC algorithm. In section 3, the molecular models (force fields) and simulation details of this study are briefly described. The new versions of the AVBMC algorithm are compared to the original one for two different scenarios: the simple ideal-association model of van Roij10 and the superheated vapor phase of hydrogen fluoride. The results of these simulations are presented and discussed in section 4, and section 5 provides concluding remarks. 2. The AVBMC Algorithm For simplicity, the different versions of the AVBMC algorithms are presented using the canonical ensemble formalism, but the extension of them to other ensembles is straight-
10.1021/jp012209k CCC: $20.00 © 2001 American Chemical Society Published on Web 10/16/2001
11276 J. Phys. Chem. B, Vol. 105, No. 45, 2001
Chen and Siepmann
forward. The classical partition function of a system in the canonical ensemble can be written as follows:11
Q)C
[ ] E(r)
∫ exp - kBT
dr
(1)
where C is the standard constant containing the integration over momenta and the conversion factor required to relate the purely classical partition function to its quantum-mechanical equivalent. The remainder is the classical configuration integral that depends only on the set of nuclear coordinates, r. The symbols E, kB, and T stand for the potential energy, Boltzmann’s constant, and the temperature, respectively. The probability distribution entailed by this partition function is
[ ]∫ [ ]
E(r) / F(r) ) exp kBT
E(r) exp dr kBT
{
acc(Aout f Bin) ) min 1,
(2)
For completeness, we include a description of the conventional Metropolis scheme.2 Starting from a given configuration of the system, state A, the conventional Metropolis scheme proceeds as follows: (i) Attempt any type of displacement of nuclear coordinates (e.g., a translation or rotation of the nuclear coordinates of a randomly selected molecule with uniformly selected displacements) to the trial configuration, state B, (ii) calculate the potential energy difference, ∆E ) EB - EA, and (iii) accept this move with a probability of
acc(A f B) ) min[1, exp(- ∆E/kBT)]
9 for a description of how to efficiently select trial sites that are in the in or out regions); (iv) calculate the potential energy difference, ∆E ) EB - EA; and (v) accept this move with the following set of acceptance probabilities: (a) for the two cases in which the swap move does not involve molecule i entering or leaving the bonded region of particle j, that is, the cases out f out or in f in, the standard Metropolis acceptance rule (see eq 3) is used; (b) if the move results in the formation of a bonded configuration of i and j (irrespective of whether i was bonded to another molecule k before the move), that is the out f in case, the following acceptance rule is used:
(3)
It is straightforward to demonstrate that the Metropolis scheme samples the probability distribution in eq 2 by imposing the condition of microscopic reversibility:12
FA R(A f B) acc(A f B) ) FB R(B f A) acc(B f A) (4) where R(A f B) ) R(B f A) which are the components of the underlying transition matrix which is symmetric for the conventional Metropolis scheme. The symmetry of the matrix results in a low efficiency of locating the few bonded configurations in the vast sea of phase space characterized by nonbonded configurations and the standard acceptance rule makes the destruction of bonded configurations with their favorable energies a very slow process. To allow efficient hopping between the bonded and nonbonded configurations, the AVBMC algorithm9 introduced a biased intrabox swap move. First, we describe the AVBMC algorithm in its original form.9 Starting from a given nuclear configuration, state A, an AVBMC intrabox swap move proceeds as follows: (i) randomly selected a molecule i to be swapped; (ii) randomly selected a second molecule j (j * i) that acts as the target for the swap move; (iii) with a probability of Pbias, molecule i is only allowed to swap into the bonded region of molecule j, called the Bin state (e.g., for hydrogen bonding fluids, the bonded regions can be conveniently defined as the maximum hydrogen-bond distance between the heavy atoms which can be estimated a priori from energetic considerations or be determined from radial distribution functions calculated over short test simulations), whereas with a probability of 1 - Pbias, the molecule i is swapped into the nonbonded region of molecule j, called the Bout state (see ref
}
(1 - Pbias)Vin exp(-∆E/kBT) PbiasVout (5)
(c) if the move results in the destruction of a bonded configuration of i and j (irrespective of whether the swap move places i in the bonded region of another molecule k), that is the in f out case, the following acceptance rule is used:
{
acc(Ain f Bout) ) min 1,
}
PbiasVout exp(-∆E/kBT) (1 - Pbias)Vin
(6)
It is evident that the underlying transition matrix of the Markov chain for the AVBMC algorithm is asymmetric, and to ensure the detail balance condition, the acceptance rules written above have taken into account the different transition probabilities for trials to place a molecule in the in region (proportional to Pbias/ Vin) and trials to place a molecule in the out region (proportional to (1 - Pbias)/Vout). Thus, it is clear that the AVBMC algorithm improves the rate of bond formations and destructions on the following two aspects: (i) the transition probability for moves that create bonded configurations has been greatly increased since the intra-box swap move is directly targeted at either the in or out region and, correspondingly, (ii) the acceptance probability for moves that destroy bonded configurations has also been greatly enhanced (compare eq 6 to 3). Unlike the conventional Metropolis algorithm, in the AVBMC algorithm, the space surrounding a molecule is explicitly divided into in (bonded) and out (nonbonded) regions. This allows for other molecules to directly hop between these two regions (using the intrabox swap moves) and thereby to achieve more efficient sampling of both bonded and nonbonded configurations. It is immediately seen from this point that the AVBMC algorithm evokes a general idea which surpasses the conventional Metropolis scheme for simulation of associating fluids, that is to (i) acknowledge the presence of different microphase regions for such heterogeneous systems formed by the associating fluids and (ii) design moves that allow efficient hopping of particles from one region to the others. The challenge to the conventional Metropolis scheme posed by the heterogeneous systems is that these different microphase regions (e.g., bonded and nonbonded phase regions) may differ to a great extent on both energetic and entropic factors, whereas the transition probability in the Metropolis scheme is solely energetically based (the Boltzmann weight of the energy difference between the trial and present configurations). This challenge is met by the AVBMC algorithm in which both Boltzmann weight (energetic factor) and volume (entropic factor) are taken into account in the transition probability. In addition, we notice that during an intrabox swap move between the two subregions (in and out) of a randomly
Improving the Efficiency of the AVBMC Algorithm
J. Phys. Chem. B, Vol. 105, No. 45, 2001 11277
chosen molecule (say molecule j) the canonical partition function could be rewritten as follows (molecule j is put at the origin): N-1
Q)CV
∑ ∫V
Nin)0
in
drNin
∫V
out
{
drNout exp -
(a) if it is an out f in move, then the following acceptance rule is used:
acc(Aout f Bin) )
1
[Uin + Uout +
kBT
interactions between the two regions]
}
(7)
where Vin and Vout are the volumes of the bonded and nonbonded region, respectively, which sum to V. Nin and Nout are the number of molecules in the bonded and nonbonded region, respectively, which sum to N - 1 (j is excluded). Uin (or Uout) is the interaction energy among the molecules inside the in (or out) region. From this partition function, one may also derive the acceptance rules for the various intrabox swap moves used by the AVBMC algorithm. Notice the similarity between this equation and the partition function of the canonical version of the Gibbs ensemble.13-15 It is evident that there is an intrinsic connection between this scheme and the Gibbs ensemble Monte Carlo method. Both schemes involve explicit division of different phase regions and utilize the same type of moves to allow particles to cross the boundary of different phase regions (interbox swap in the Gibbs ensemble method versus intrabox swap in AVBMC). However, there are some differences between these two methods. In the Gibbs ensemble, the interactions between the different phase regions are actually zero because they are physically separated. This enforced separation of the different phase regions and the fact that each phase region is embedded into its own replicas with periodic boundary conditions also imply that the Gibbs ensemble is a suitable approach only when phase separation is on a macroscopic length scale, not for the systems which contain coexisting microphase regions that are only microscopic in length scale (e.g., aggregates or micelles). Whereas, AVBMC can be applied to both types of systems. As shown above, in the original AVBMC scheme (from now on called the AVBMC1 algorithm), there are four types of intrabox swap moves: out f in, in f out, in f in, and out f out. The first two types are supposed to contribute mainly to the formation and destruction of bonded configurations, whereas the latter two types may not help much in sampling compared to the conventional random swap moves. In addition, the majority of these intrabox swap moves belong to the out f in and out f out cases simply because of the very low probability of finding a particle i (randomly) which happens to be in the bonded region of the target particle j (see the intrabox move sequence above). These shortcomings greatly lower the efficiency of the AVBMC1 algorithm, which has also been pointed out by Wierzchowski and Kofke.8 To overcome these shortcomings, one can eliminate all of the in f in and out f out moves by rewriting the intrabox swap move sequence as follows: (i) Randomly select a particle j that acts as the target for the swap move; (ii) with a probability of Pbias, the move type is selected as out f in, whereas with a probability of 1 - Pbias, the move type is in f out; (iii) if it is selected as an out f in move, then a particle is chosen randomly in the out region of particle j and placed randomly inside the in region of particle j, or otherwise, if it is an in f out move, then a particle is chosen randomly in the in region of particle j and placed in the out region of particle j; (iv) calculate the potential energy difference between the initial state A and the trial state B, ∆E ) EB - EA; and (v) accept this move with the following set of acceptance probabilities:
[
min 1,
]
(1 - Pbias)VinNout exp(-∆E/kBT) PbiasVout(Nin + 1)
(8)
(b) if it is an in f out move, then the following acceptance rule is used:
[
acc(Ain f Bout) ) min 1,
]
PbiasVoutNin exp(-∆E/kBT) (1-Pbias)Vin(Nout + 1)
(9)
It should be noted that these acceptance rules are similar to the one used for the interbox swap moves in the Gibbs ensemble13-15 and the proof of the detailed balance condition is straightforward. Thus, this new version of the AVBMC algorithm (called AVBMC2) restricts all of the intrabox swap moves to either in f out moves or out f in moves, which leads to a more frequent sampling of bond formations and destructions. Because there is no limitation on the number of regions into which one can divide the whole space, the AVBMC algorithm can be extended in a similar spirit as the multiphase Gibbs ensemble.16 The AVBMC3 algorithm introduced below employs one intrabox swap move sequence (among many others) where a total of three regions are considered: (i) Randomly select a particle j and another particle k with j * k and their in regions nonoverlapping (otherwise, randomly select particle k again until the condition is satisfied), (ii) with a probability of Pbias, randomly select one particle either from the in region of particle k or from the out region of particle j (with equal probabilities, where the out region of particle j contains the in region of particle k, too) and move it to the in region of particle j, whereas with a probability of 1 - Pbias select one particle randomly inside the in region of particle j and move it either to the in region of particle k or to the out region of particle j (with equal probabilities), (iii) calculate the potential energy difference between the initial state A and the trial state B, ∆E ) EB EA, and (iv) accept this move with the following set of acceptance probabilities: (a) if it is an out f in move or an in f out move, the acceptance rules used in the AVBMC2 algorithm apply (see eqs 8 and 9) with Vjout ) V - Vjin - Vkin. (b) if it is an in f in move, e.g., move one particle from the in region of particle j to the in region of particle k, the following acceptance rule is used:
[
acc(Ain f Bin) ) min 1,
]
PbiasVkinNjin exp(-∆E/kBT) (1 - Pbias)Vjin(Nkin + 1)
(10)
(c) if it is a particle swap from the in region of particle k to the in region of particle j, the indices j and k as well as the biasing probabilities Pbias and 1 - Pbias in eq 10 need to be switched. Similar to the AVBMC2 algorithm, the probabilities of attempting in f out and out f in moves are enhanced in the AVBMC3 algorithm compared to the AVBMC1 algorithm. However, the main advantage of the AVBMC3 algorithm is that it allows for more frequent sampling of cluster hopping through the in f in moves between aggregates centered at particles j and k compared to both the AVBMC1 and AVBMC2 algorithms. It should be clarified that for the AVBMC2 and AVBMC3 intrabox swap moves, when the moving particle is selected from
11278 J. Phys. Chem. B, Vol. 105, No. 45, 2001
Chen and Siepmann
the in region of one particle (say molecule j), it is randomly selected from particle j’s neighbors inside the in region excluding particle j itself. If particle j does not have any neighbors, then this move is immediately rejected. Thus, compared to the AVBMC1 algorithm, in both AVBMC2 and AVBMC3 simulations, one needs to keep track of the neighbors of each molecule. Fortunately, this adds only a very small amount of computational cost and programming effort to the original AVBMC1 code because all of the pair distances required to determine the neighbors have been computed in the summation of the total energy at the beginning of the simulation and can be updated in the energy calculation during each trial attempt. 3. Simulation Models and Details A. The Ideal-Association Model (IAM). The IAM was initially introduced by van Roij10 and later used by Kofke and co-workers for testing their monomer-addition-subtraction algorithm (MASA)7 and UB method.8 An IAM molecule consists of two sites A and B, which are separated by a fixed distance. The only interactions defined in the IAM model are a squarewell type of attraction between an A site of one molecule and a B site of another molecule as follows
uij )
{
- 0
if rAi,Bj < Reff if rAi,Bj g Reff
(11)
An analytical expression of the aggregate distribution in the canonical ensemble has been provided by van Roij for this model under the constraint that branched or ring structures are disallowed:8,10
Ni )
i N 1 N 1 - (H - 1) ) Yi A 2A A
[
]
(12)
where Ni and N are the number of molecules in i-mers and the total number of molecules, respectively. H ) x4A+1, and A ) FVeff exp(β) where F is the number density defined by the ratio of the total number of molecules N to the total volume of the system V, Veff is the bonding volume, 4/3πR3eff, and β is the inverse temperature. Following from this aggregate distribution, the energy per molecule is described by
u)-
N
(i - 1)Yi ∑ N i A
(13)
Accordingly, the constant-volume heat capacity is given by
CV NkB
)
β22 A
[ (
∑i (i - 1) iYi-1
H-1 2A
-
1
) ]
H
- Yi
(14)
To avoid branching or ring structures in the simulations, we followed the approach by Wierzchowski and Kofke;8 that is, those attempts that lead to formation of branched and ring structures are immediately rejected. It should be noted that this treatment differed slightly from van Roij’s original approach, in which no additional energy is gained when extra B or A sites move into a nonvacant A or B site. It has been demonstrated by Wierzchowski and Kofke8 that this difference causes only a minor deviation for the energy and even a smaller effect on the heat capacity from van Roij’s original solution at the relatively low densities of interest. It should be noted that both sums in eqs 13 and 14, in principle, run from 1 to infinity (all possible
clusters), but in the simulations, the total number of molecules N is the upper bound (the largest possible cluster). Following Wierzchowski and Kofke,8 the simulated system consisted of 256 IAM molecules and the intramolecular separation RAB was set to 4Reff. Two sets of simulations were carried out. In the first set of simulations, the energy convergence was monitored using the three different versions of the AVBMC algorithm or the conventional Metropolis scheme starting from two different initial configurations: (i) the molecules are spaced far apart with initial energy of zero and (ii) the molecules are arranged as a linear n-mer with a total energy of (N - 1). The parameters β and A were set to 8 and 0.154, respectively. The second set of simulations were carried out using only the AVBMC algorithm at FVeff ) 3.43 × 10-5 and β ) 9.6, 10.4, 11.2, 12.0, 12.8, 13.6, 14.4, 15.2, 16.0, and 16.8. Each simulation included an equilibration of 200 000 Monte Carlo cycles followed by a production of 500 000 Monte Carlo cycles, and the statistical uncertainties were estimated by dividing the simulations into 10 blocks. The type of Monte Carlo move was selected at random according to the following probabilities: 0.25 (translation), 0.25 (rotation), and 0.5 (swap). For translations and rotations, the maximum displacements were adjusted to yield acceptance probabilities of 50%. It should be mentioned that the IAM model allows for a straightforward definition of the bonded region, the two spheres centered at sites A and B with a radius of Reff. This is the exact in region used by the AVBMC simulations. Considering site A (or B) of a molecule can only be associated with site B (or A) of another molecule, the two associating regions of the target molecule j (see the description above) were randomly chosen with equal probabilities as the in region, and thereafter the regular AVBMC intrabox swap move sequence proceeds. It should be noted that for an in f out move, only the neighbor associated with the selected in region (thus at most one neighbor) is chosen to move out, and for an in f in move in the AVBMC3 simulations once the in region for molecule j is selected, the in region for molecule k is also set. B. Superheated Hydrogen Fluoride Vapor. The OPLS (optimized potentials for liquid simulation) hydrogen fluoride model17 was used in this work. It consists of three Coulombic interaction sites, one on the fluorine nucleus and one on the hydrogen nucleus with equal partial charges of +0.725 e and one extra bond site with a compensating negative charge of -1.45 e. The hydrogen and fluorine nuclei are kept at a fixed separation of 0.917 Å, and the bond site is placed 0.166 Å away from the fluorine toward the hydrogen. There is an additional Lennard-Jones interaction at the fluorine site with /kB ) 75.75 K and σ ) 2.984 Å. A molecule-based potential truncation at fluorine-fluorine separations of 8 Å was used for both LennardJones and Coulombic interactions and long-range corrections were applied only for the LJ part of the potential. The system consisted of 108 hydrogen fluoride molecules. The initial configuration was constructed using a layered crystal structure. Isobaric-isothermal Monte Carlo simulations18 were performed over a range of pressures and temperatures using 5 × 104 and 106 Monte Carlo cycles for the equilibration and production periods at each state point (one Monte Carlo cycle consists of N ) 108 Monte Carlo moves). 4. Results and Discussion A. The IAM. As pointed out by Kofke and co-workers,7,8 the IAM is an ideal model representation for associating systems and therefore has been used often for testing the efficiency of various biased algorithms developed for simulations of associ-
Improving the Efficiency of the AVBMC Algorithm
Figure 1. Evolution of the energy per particle versus the CPU time for 256 IAM molecule at β ) 8 and A ) 0.154 starting from a noninteracting lattice arrangement (top) or from a fully interacting IAM chain (bottom). A single node of a Pentium III at 866 MHz computer was used. Results are shown for the simulations using the conventional Metropolis scheme (squares), the AVBMC1 (diamonds), the AVBMC2 (dotted lines), and the AVBMC3 algorithm (dashed lines). The analytical solution by van Roij10 is depicted as solid line.
ated fluids. There are two important parameters in this model, β and A, which can be used to adjust the strength of association for the system and, correspondingly, the level of difficulty for a conventional unbiased algorithm to obtain a converged answer for this system.7,8 First, we carried out the energy convergence test using the AVBMC algorithms (all three versions presented in section 2) or the conventional Metropolis algorithm on 256 IAM molecules with β ) 8 and A ) 0.154. This state condition was also used by Wierzchowski and Kofke in an energy convergence study of several biased algorithms, including the original version of the AVBMC algorithm and the UB algorithm.8 Following from their procedures, the simulations either started from an IAM system with zero interaction (the molecules were arranged on a lattice and thus far apart from each other) or started from a fully interacting IAM system (all the molecules participate in a chainlike aggregate and with the exception of the two ends, each molecule interacts with two others). For comparison with the simulations by Wierzchowski and Kofke,8 the energy convergence was plotted against the CPU time in Figure 1 (each point shows the average energy over a block size of 100 Monte Carlo cycles, and the overall simulation length is 10 000 Monte Carlo cycles). First, all three AVBMC algorithms allowed for convergence to the analytical solution, whereas the Metropolis scheme did not converge for the simulation started from the fully interacting configuration. For the simulations using the AVBMC1 algorithm, the energy convergence took about 30 CPU seconds for the noninteracting starting configuration and 50 CPU seconds for the fully interacting starting configuration. Wierzchowski and Kofke
J. Phys. Chem. B, Vol. 105, No. 45, 2001 11279 reported equilibration times of about 75 and 150 CPU seconds for their version of AVBMC1 and on a different computer. In addition, even after an equilibrium configuration was reached, the relatively large fluctuations of the block-averaged energies shown by the AVBMC1 simulation demonstrate that 100 Monte Carlo cycles are insufficient to yield a precise answer. In contrast, both AVBMC2 and AVBMC3 simulations required less than 4 CPU seconds to converge to an equilibrium distribution. These timings show a gain of about 1 order of magnitude in efficiency over the AVBMC1 algorithm and are also shorter than the timings reported by Wierzchowski and Kofke8 for the UB algorithm (these were about 40 and 85 CPU seconds for the noninteracting and fully interacting starting configurations, respectively). However, a direct comparison is difficult because of the different choices of the bonding region, see ref 8. It should be mentioned that after the short equilibration, a precise average energy can be obtained using a length of about 100 Monte Carlo cycles for both AVBMC2 and AVBMC3 simulations. It is also interesting to note that both AVBMC2 and AVBMC3 simulations finish earlier than the simulations employing the AVBMC1 algorithm or the Metropolis scheme (by about 30-40 CPU seconds) despite the increase of the computational overhead mentioned above. The reason behind this is that in both AVBMC2 and AVBMC3 simulations, most of the in f out (and also in f in for AVBMC3) intrabox swap moves are immediately rejected due to the lack of neighbors (notice that at this state condition, a majority of IAM molecules are monomers). The second set of simulations were carried out at more severe conditions with β ranging from 9.6 up to 16.8 and A equal to 3.43 × 10-5. Here, only the three versions of the AVBMC algorithm were used. This set of state conditions was also used by Wierzchowski and Kofke,8 and they found that the unbiased Metropolis algorithm failed for most cases. In Figure 2, the energies of the IAM systems are plotted as function of β. As evidently shown in this plot, all three AVBMC versions are able to yield average energies that are very close to the analytical solution, but a small deviation was observed for both AVBMC1 and AVBMC2 at β ) 16.8, where most of the molecules are associated (viz. the energy per molecule is close to - at this condition). The corresponding heat capacities are plotted versus β in Figure 3. The heat capacity can be calculated either from finite differences between simulations at adjacent temperatures or using the fluctuation formula.9,12 As the calculations of fluctuation properties are very slowly converging and require a sufficient exploration of phase space,2,12 we followed the lead of Visco and Kofke7 and used the fluctuation formula to calculate the heat capacity to demonstrate the sampling efficiency of the AVBMC algorithms. In contrast with the energy plot in Figure 2, both AVBMC1 and AVBMC2 start to deviate from the analytical solution at about β ≈ 13 (beyond the heat capacity maximum, a transition point from the monomerdominated region to the aggregate-dominated region). The AVBMC3 works exceptionally well, showing only a small divergence at β ) 16.8. This is comparable to the results obtained with the UB algorithm by Wierzchowski and Kofke.8 This points to the importance of more frequent sampling of aggregate hoppings through the particular intrabox in f in swap move between aggregates, in particular when aggregates (bonded configurations) are dominant. Although a majority of intrabox out f in swap moves in the AVBMC1 and AVBMC2 simulations for higher values of β were attempted hops between aggregates (because most molecules are associated), these moves were rarely accepted due to the huge entropic penalty involved
11280 J. Phys. Chem. B, Vol. 105, No. 45, 2001
Figure 2. Plot of the energy per particle as function of reduced inverse temperature for 256 IAM molecules at FVeff ) 3.43 × 10-5. Results are shown for the simulations using the AVBMC1 (diamonds), the AVBMC2 (squares), and the AVBMC3 algorithm (circles). The analytical solution by van Roij10 is depicted as solid line.
Figure 3. Plot of the heat capacity per particle as function of reduced inverse temperature for 256 IAM molecules at FVeff ) 3.43 × 10-5. Symbols as in Figure 2.
in the acceptance rule in the absence of the balance of the energetic factor.
Chen and Siepmann
Figure 4. Plot of the internal energy as function of temperature for the superheated hydrogen fluoride vapor at 15.5 kPa. Results are shown for the simulations using the AVBMC1 (diamonds), the AVBMC2 (squares), and the AVBMC3 algorithm (circles). Lines are drawn only as guides for the eye.
B. Superheated Hydrogen Fluoride Vapor. The superheated vapor phase of HF was previously used as a test case by Visco and Kofke7 for their MASA algorithm and also used by us for the original AVBMC algorithm.9 Here we repeated the calculations for the 15.5 kPa isobar using the OPLS HF model17 for the AVBMC2 and AVBMC3 algorithms. It should be mentioned that the bonded Vin region was bounded by an inner shell at 2.0 Å and outer shell at 4.0 Å based on the F-F separation, which differ from 5.0 Å used in the previous AVBMC1 simulations. Pbias was set to 0.5. As in the earlier AVBMC1 simulations,9 the configurational-bias Monte Carlo technique19-21 with multiple insertions22,23 of the first bead (nfirst choice ) 10) and multiple orientations (nchoice ) 10) was also used to enhance the acceptance rates for the AVBMC2 and AVBMC3 swap moves. The energies of the superheated HF vapor are plotted as a function of temperature in Figure 4. All three versions of the AVBMC algorithm yielded answers that agree with each other within the uncertainties. It is evident that at the lowest temperature (T ) 240 K) aggregates (bonded configurations) are the dominant species (with a very large negative molar energy of about -20 kJ/mol), whereas at the highest temperature (T ) 310 K), the majority of the molecules are found as monomers (nonbonded configurations with a molar energy close to zero). There is a sharp turning point at about 260 K, which denotes the transition from the aggregate-dominated region to the monomer-dominated region. This turning point is also reflected by the appearance of a heat capacity maximum at this temperature (see Figure 5). The heat capacity was again computed from the fluctuation formula. For this case, all of the three versions of the AVBMC algorithm yielded satisfactory results, which is not surprising since Wierzchowski and Kofke8 discussed that the HF vapor at this condition is comparable to an IAM model at a state condition with β less than 8.
Improving the Efficiency of the AVBMC Algorithm
Figure 5. Plot of the heat capacity as function of temperature for the superheated hydrogen fluoride vapor at 15.5 kPa. Symbols as in Figure 4.
The populations of aggregates (not shown here, but see ref 9) were analyzed using a simple distance-based criterion (with a F-F distance of less than 3.3 Å) for all eight temperatures used in this study. For an aggregate with a certain size, say n-mer, its population is temperature dependent. From a “chemical” point of view, by taking the n-mer as a new chemical species, one can estimate its enthalpy of formation, ∆H, from the temperature dependence of the population equilibrium constants via a van’t Hoff plot.24 In addition, the enthalpy of formation can be also calculated directly from our simulations using the relation ∆H ) ∆U + p∆V with ∆U ) Un - nU1 (where Un and U1 denote the energies of an n-mer and a monomer, respectively). It should be noted that these energies also contain contributions from the environment, but these are very small for a dilute gas phase. The term p∆V is approximated by (n - 1)RT (∆V could be computed directly from a population-pressure dependence of n-mers by applying another fundamental thermodynamic equation,24 but it would require several simulations at various pressures along an isotherm.). As shown in Figure 6, the two methods yield enthalpies of aggregate formation for the dimer to the tetramer (because of lower populations, the results were not estimated for the pentamer and above) that are in excellent agreement. Thus, both “chemical” and “physical” views of these aggregates are consistent.1 In addition, the enthalpies of formation of both trimer and tetramer are noticeably lower for the 15.5 kPa isobar than for the 96.1 kPa isobar. Both trimer and tetramer can exist in various architectures (e.g., linear and cyclic architectures for the trimer and, in addition, branched structures for the tetramer). These architectures share different energetic and entropic characteristics (cyclic structures are energetically favorable, whereas branched structures are entropically favorable). Because the enthalpies of formation were determined at a lower temperature range (240-310 K) for the 15.5 kPa isobar compared to the 96.1 kPa isobar (270-340 K), the isomers which are energetically more
J. Phys. Chem. B, Vol. 105, No. 45, 2001 11281
Figure 6. Plot of the enthalpy of aggregate formation as function of aggregate size for the superheated hydrogen fluoride vapor at 15.5 kPa (dashed line and squares) and 96.1 kPa (solid line and circles). Lines and symbols represent the results estimated from the macroscopic route (the van’t Hoff equation) and from the microscopic route (the energy analysis), respectively.
favorable have a higher population for the 15.5 kPa isobar which lowers the overall enthalpies of formation. Finally, it should be noted that both curves are not linear, which points to the cooperativity of hydrogen-bonded aggregate formation. In particular, the gains in the enthalpies of formation from dimer to trimer and from trimer to tetramer (being 36 and 42 kJ/mol, respectively) are significantly larger than from monomer to dimer (21 kJ/mol). 5. Conclusions The underlying idea of the AVBMC algorithm was reanalyzed. It was found that the advantages of the AVBMC algorithm over the conventional Metropolis scheme for simulations of associating fluids stem from the fact that the AVBMC algorithm distinguishes various microphase regions that can coexist in microheterogeneous, self-assembling systems. The AVBMC algorithm is able to balance the different energetic and entropic factors exhibited by these microphase regions (i.e., bonded and nonbonded regions) by (i) attempting more frequently moves that lead to low-entropic but high-energetic regions and (ii) enhancing the acceptance rates for moves that lead to high-entropic but low-energetic regions. Therefore, it allows more efficient sampling of all of the relevant microphase regions. From this analysis, two improved versions, called AVBMC2 and AVBMC3 algorithms, were presented. These new versions are as simple and general as the original AVBMC algorithm but are more efficient and robust. The AVBMC2 algorithm improves the efficiency by limiting the intrabox swap moves to either in f out moves or out f in moves. The AVBMC3 algorithm further improves on the AVBMC2 version by allowing more frequent sampling of particle hopping between
11282 J. Phys. Chem. B, Vol. 105, No. 45, 2001 aggregates through the addition of the in f in move between two different aggregates. There is a small increase of the computational overhead for both AVBMC2 and AVBMC3, that is, the need to update a list of the neighbors of each molecule. The AVBMC3 algorithm excels for super-strongly associating fluids where most of the molecules are found in aggregates of different sizes and where the equilibrium between cluster sizes needs to be sampled. However, the AVBMC2 algorithm performs well for systems which have a sufficient fraction of monomeric species or which show a strong preference for small aggregates of a specific size (e.g., dimers of carboxylic acids). The three versions of the AVBMC algorithm were applied to two cases: the IAM and the vapor phase of superheated hydrogen fluoride. It is demonstrated from the tests on the IAM that the AVBMC2 and AVBMC3 algorithms can be orders of magnitude more efficient than the AVBMC1 algorithm, whereas the AVBMC3 algorithm becomes superior when aggregates of varying size dominate. For the superheated HF vapor phase, in general, good agreement was found between the simulations using the different versions of the AVBMC algorithm. It was also found that the enthalpies of aggregate formation estimated from a macroscopic route (the van’t Hoff equation) and from a microscopic energy analysis match very well. In addition, the temperature dependence of the enthalpies for trimer and tetramer formation could be explained by the fact that both trimer and tetramer are mixtures of several isomers whose populations are temperature dependent. It was found that the enthalpies of formation decreased by larger amounts from the dimer to the trimer and from the trimer to the tetramer than from the monomer to the dimer, which emphasizes that cooperativity effects play an important role for these hydrogen-bonded aggregates. Acknowledgment. We are grateful to David Kofke for providing us with a copy of the UB manuscript prior to its publication and for stimulating discussions on the original
Chen and Siepmann AVBMC algorithm. Financial support from the National Science Foundation (CTS-9813601), a Sloan Research Fellowship, and a Doctoral Dissertation Fellowship (B.C.) is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute. References and Notes (1) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria; 2nd edition; Prentice Hall: New York, 1986. (2) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (3) Tsangaris, D. M.; de Pablo, J. J. J. Chem. Phys. 1994, 101, 1477. (4) Busch, N. A.; Wertheim, M. S.; Chiew, Y. C.; Yarmush, M. L. J. Chem. Phys. 1994, 101, 3147. (5) Busch, N. A.; Wertheim, M. S.; Yarmush, M. L. J. Chem. Phys. 1996, 104, 3962. (6) Visco, D. P., Jr.; Kofke, D. A. J. Chem. Phys. 1998, 109, 4015. (7) Visco, D. P., Jr.; Kofke, D. A. J. Chem. Phys. 1999, 110, 5493. (8) Wierzchowski, S.; Kofke, D. A. J. Chem. Phys. 2001, 114, 8752. (9) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8725. (10) van Roij, R. Ph.D. Thesis, University of Utrecht, Utrecht, The Netherlands, 1996. (11) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (12) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: Oxford, U.K., 1987. (13) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (14) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527. (15) Smit, B.; de Smedt, P.; Frenkel, D. Mol. Phys. 1989, 68, 931. (16) Lopez, J. N. C.; Tildesley, D. J. Mol. Phys. 1997, 92, 187. (17) Cournoyer, M. E.; Jorgensen, W. L. Mol. Phys. 1984, 51, 119. (18) McDonald, I. R. Mol. Phys. 1972, 23, 41. (19) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59. (20) Vlugt, T. J. H.; Martin, M. G.; Smit, B.; Siepmann, J. I.; Krishna, R. Mol. Phys. 1998, 94, 727. (21) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508. (22) Esselink, K.; Loyens, L. D. J. C.; Smit, B. Phys. ReV. E 1995, 51, 1560. (23) Mackie, A. D.; Tavitian, B.; Boutin, A.; Fuchs, A. H. Mol. Sim. 1997, 19, 1. (24) Castellan, G. W. Physical Chemistry; 3rd edition; Benjamin/ Cummings: Menlo Park, CA, 1983.
