J. Phys. Chem. 1996, 100, 1153-1158
1153
Free Energy Calculations by Expanded Ensemble Method for Lattice and Continuous Polymers P. N. Vorontsov-Velyaminov,* A. V. Broukhno, T. V. Kuznetsova, and A. P. Lyubartsev† Faculty of Physics, St. Petersburg State UniVersity, 198904, St. Petersburg, Russia ReceiVed: May 9, 1995X
Expanded ensemble Monte Carlo method suggested earlier in J. Chem. Phys. (1992, 96, 1776) was applied to calculate free energy (entropy) of a lattice polymer (self-avoiding random walk on a simple cubic lattice) with respect to the standard reference systemsa pure random walk with no restrictions (in the latter case the number of conformations is known exactly). An effective scheme for the initial choice of balancing factors is proposed with the following iteration procedure for their optimization. Calculations were performed for chain lengths N ) 5, 8, 25, 50, and 100. Comparison of simulation results with the exact data available for short chains (N ) 5, 8) indicate a high accuracy of the simulation method. For the continuous polymer model (several bead and spring chains in the periodic box of a fixed volume) the expanded ensemble Monte Carlo method was used to calculate free energy with respect to the ideal gas as a reference system. Configurations were changed with the aid of the constant temperature molecular dynamics procedure, while changes in expansion parameter were made within the usual Monte Carlo scheme. For several cases (a single chain of N ) 32 monomers in the periodic box, two chains, and four chains), the balancing factors were optimized and free energies were calculated in a wide range of temperatures.
Introduction Knowledge of free energy is crucial for understanding such molecular phenomena as phase transitions, solvation, and conformational transitions in macromolecules, chemical equilibria, etc. Monte Carlo (MC) or molecular dynamics (MD) simulations are straightforward in calculating such characteristics of condensed molecular systems as internal energy, pressure, and radial distribution functions for which there exists an appropriate microscopic analoguesan estimator. The Helmholtz free energy and other quantities which include entropy cannot practically be obtained by simple (direct) averaging of an appropriate estimator. For calculation of free energy within MC or MD simulations a number of methods were worked out. They are thermodynamical integration,1-3 particle insertion,4,5 multistage sampling,6 different versions of umbrella sampling,7,8 perturbation method,2,9,10 acceptance ratio method,11,12 and some others (see also reviews13,14). However, most of them require a large number of computer runs, face difficulties, and even fail for systems with a strong coupling parameter or at high densities. Recently a new approach, the expanded ensemble (EE) method, was suggested15 (see also preceding papers16). Similar types of approaches have also been reported in other works: semigrand-canonical ensemble for chemical potential calculation,17 simulated tempering,18 and force balance method.19 The expanded ensemble method allows us to calculate absolute free energies with a high precision for arbitrary systems using a rather simple algorithm within the MC15 or MD20 scheme. The method was tested and proved its effectiveness for different systems: primitive electrolyte model,15,21 Lennard-Jones fluids and water,20 quantum Heisenberg model,21,22 or ordered polyelectrolyte system.23 In the present work the expanded ensemble method is used for calculation of the free energy (entropy) of a polymer model, † Present address: Division of Physical Chemistry, Arrhenius Laboratory, University of Stockholm, S-10691, Stockholm, Sweden. X Abstract published in AdVance ACS Abstracts, December 1, 1995.
0022-3654/96/20100-1153$12.00/0
both lattice and continuum. In the lattice case a polymer is presented as a self-avoiding random walk on a lattice. The partition function Z in this case is the total number of configurations, and the entropy is
S ) -(βF) ) ln(Z)
(1)
There are no clear ways to calculate Z for long enough polymers. For short chains one can count all configurations exactly. For longer chains one can randomly construct fully penetrable (Gaussian or ideal) chains and then count the ratio of nonoverlapping chains, but this ratio rapidly tends to zero with an increase in the chain length, thus making this way also ineffective. The idea is to construct “an umbrella” of semipenetrable chains, which fills the gap between the self-avoiding and the ideal chain, and use the expanded ensemble to calculate the ratio of partition functions or the entropy difference. This idea was first suggested in ref 15. Recently this approach was successfully applied for calculation of the chemical potential in a dense polymer system.24 The existence of exact solution for short lattice polymers makes this system very attractive for methodological study of the expanded ensemble method. In the following sections we give a short description of this method for entropy calculation of lattice polymers together with a scheme for optimization of algorithm parameters. Then we present results of calculations of the lattice polymer model and comparison with exact results for a short chain. In the last section we apply the expanded ensemble method to the continuous polymer model. Expanded Ensemble Method for a Lattice Polymer The general formulation of the expanded ensemble method has been presented elsewhere.15 Here we describe its application to the polymer model. Consider a “penetrable” lattice polymer with total potential energy
U(qi) ) N(qi) © 1996 American Chemical Society
(2)
1154 J. Phys. Chem., Vol. 100, No. 4, 1996
Vorontsov-Velyaminov et al.
where N(qi) is the number of overlaps in the given conformation (the energy of a single overlap is taken equal to unity). The reduced Hamiltonian is
H(qi) ) βU(qi) where β is the reciprocal temperature. In the case β ) 0 we have an ideal polymersa random walksfor which the configurational partition function and free energy is known exactly; for sufficiently large β (β . 1, or T f 0) we have a self-avoiding random walk, i.e., a polymer with excluded volume. Now we consider a set of canonical ensembles with reduced reciprocal temperatures βm ) 1/kTm:
0 ) β0 < β1 < β2 < ... < βM . 1 The canonical partition function (here and below we always mean the configurational partition function) for each m is
Zm ) ∑ exp(-Hm(qi))
(3)
i
where the sum is taken over all possible conformations {qi}. Further we construct a partition function for the expanded ensemble as
Z)
∑
Zm exp(ηm)
(4)
0emeM
where ηm are so-called balancing factors to be determined later. The MC procedure includes two types of steps:
βmFm ) -L ln(6) - ln(Pm/P0) + ηm
(7)
that for a sufficiently large value of βm is practically equal to entropy S for a self-avoiding lattice polymer. Hereafter “free energy” will mean excess free energy, i.e., the free energy difference between the given system and the ideal one, as well as “entropy” will mean excess entropy. Optimal Choice of Balancing Factors In principle the result does not depend on balancing factors ηm, but they strongly (exponentially) affect probabilities Pm. The probabilities can be estimated in a MC procedure only if they are not vanisingly small; moreover it is evident that for minimizing the statistical error it is desirable to have the Pm as close to each other as possible. Unfortunately the optimal values of ηm, yielding equal probabilities for all subensembles, are unknown a priori. It is easy to show that these optimal values are
ηm ) βmFm
(8)
that are just the quantities which we want to calculate from simulations. Nevertheless it is easy to organize a simple iterative procedure for obtaining suitable values of ηm. This procedure was described in detail in previous works.15,20 The weak point of this iterative procedure is the initial choice of the trial balancing factors. Here we suggest an algorithm for obtaining the initial values of ηm. For each subensemble we made a separate short MC run and roughly estimated the internal energy U(βm) ) 〈N(q)〉m. Then for ηm the following expressions were used:
(1) change of conformation at fixed m with transition probabilities
ηm(0) - ηm-1(0) ) (βm - βm-1)(〈N(q)〉m + 〈N(q)〉m-1)/2 (9)
w(βm ) constant, qifqj) ) min{1, exp[-βm(N(qj) - N(qi))]}
with η0(0) ) 0 (the index (0) means the zeroth approximation in the following iterational procedure). This formula is based on the Gibbs-Helmholtz equation for internal and free energies:
(2) change of “temperature” at a fixed conformation with probabilities of transitions w(qi ) constant, βmfβm(1) ) min{1, exp[-N(qi)(βm(1 - βm) + ηm(1 - ηm]}
U(β) ) ∂(βF(β))/∂β
(10)
Really (10) means that for U(β) being a smooth function of β and any pair β2 > β1 we can write
β2F2 - β1F1 ) ∫β U(β) dβ ) U(β*)(β2 - β1) β2 1
In the course of a MC run we calculate nmsthe number of times when the system is in the mth subensemble. Then the probability of occurrence in this subensemble can be estimated as Pm ) nm/n (n is the total number of MC steps). On the other hand, it is clear from (4) that
Pm ) Zm exp(ηm)/Z
(5)
and hence
Pm/Pk ) (Zm/Zk) exp(ηm - ηk) ) exp[-(βmFm - βkFk)] exp(ηm - ηk) (6) So for the free energy difference between two arbitrary states (subensembles) we obtain
βmFm - βkFk ) - ln(Pm/Pk) + ηm - ηk As far as, for the ideal system (m ) 0), the partition function is known (Z ) 6L for a 3-dimensional simple cubic lattice where L + 1 is the number of monomers), and η0 can be chosen equal to zero, we can calculate
where β1 < β* < β2. For small (β1 - β2), the dependence of U can be considered linear, so
β2F2 - β1F1 ) (β1 - β2)[U(β1) + U(β2)]/2 Taking into account (8) and (2), we get the initial estimation of balancing factors (9) for each pair of neighboring subensembles. Actually (9) means that we use the thermodynamic integration for the initial rough estimation of the free energy and hence of the initial (and more or less close to optimal) values of the balancing factors. At the second stage we use the expanded ensemble to finally optimize ηm and calculate the free energies more precisely. At the second stage, optimization of ηm was achieved by means of runs in expanded canonical ensemble with application of the iterative procedure:
ηm(1) - ηk(1) ) ηm(0) - ηk(0) - ln(Pm(1)/Pk(1)) As far as for any iteration we could consider η0(i) ) 0, the
Free Energy Calculations for Selected Polymers
J. Phys. Chem., Vol. 100, No. 4, 1996 1155
TABLE 1: Example of a Procedure for the Optimal Choice of Parameters ηm for L ) 25a iteration separate runs βm 0.0 0.15 0.45 0.9 1.5 2.25 3.15 4.2 5.4 6.75 8.25 statistical error
〈Nm〉(0) 7.53 5.27 3.62 2.13 1.08 0.4 0.3 0.11 0. 0. 0.
first (0)
ηm
0.0 0.96 2.29 3.58 4.55 5.10 5.42 5.63 5.7 5.7 5.7
(1)
Pm
0.0766 0.0807 0.0807 0.0902 0.0885 0.0872 0.0843 0.0931 0.100 0.106 0.113
second (1)
ηm
0.0 0.934 2.266 3.504 4.476 5.039 5.371 5.536 5.566 5.54 5.505 0.05
(2)
Pm
0.0978 0.0923 0.0934 0.0835 0.0811 0.0834 0.0939 0.1067 0.975 0.088 0.0824
final results (2)
ηm
0.0 0.99 2.33 3.587 4.519 5.067 5.35 5.448 5.513 5.543 5.569 0.05
Pm
βmFm
〈Nm〉
βmUm
0.093 0.0935 0.0927 0.0939 0.0948 0.0923 0.0915 0.0863 0.0865 0.0871 0.0885
0.0 0.985 2.333 3.577 4.50 5.075 5.366 5.523 5.585 5.608 5.618 0.005
7.315 5.424 3.504 2.054 1.104 0.511 0.193 0.064 0.017 0.0033 0.0044
0.0 0.814 1.577 1.85 1.656 1.149 0.609 0.269 0.092 0.0224 0.0005
a 〈N 〉(0) are the average number of overlaps (internal energy) obtained in preliminary separate short runs (only 1000 MC steps in each subensemble). m ηm(0) are calculated using (9). Then two iterations in the expanded ensemble (20 000 MC steps each) follow which produce final optimal values of ηm(2). The final run was of 250 000 MC steps.
simplest way was to compute the difference between mth and ideal systems:
ηm(i) ) ηm(i-1) - ln(Pm(i)/P0(i))
(11)
Each iteration step makes values of Pm closer and closer to each other. This is illustrated in Table 1 for the case of a polymer chain of the length L ) 25 (details of the simulations are given below). When the difference in Pm is about several percent, these optimized ηm are used to perform the long final MC run in the expanded ensemble, yielding very accurate values of the free energy. Computational Details The maximum reduced temperature βM was chosen so as to make transitional probability w(βm ) constant, N(q)fN(q)+1) sufficiently small in the Mth subensemble. For example, βM ) 8 gives w ) exp(-8) ) 0.0003, that makes such an subensemble practically coincident with the ensemble for self-avoiding random walks. Values βM ) 8.25 and 9.09 were used in different calculations. The number of subensembles M ) 11 was chosen to provide high enough probabilities of transitions between neighboring subensembles. For the differences (βm - βm-1) we used the linear dependence on m. The change of conformation was carried out by choosing a monomer at random and replacing a small neighboring piece of the lattice polymer (three to five monomers) by a new piece of the same length. The remaining “tail” was added to the new conformation by a parallel shift (without changing its conformation). The procedure was tested by comparison of MC results obtained in a single subensemble at β ) 0 with corresponding exact results which can be obtained for a short chain. For example we obtained perfect coincidence of distribution over the number of overlaps N for the polymer of length L ) 8. For the given model this distribution is in fact an energy distribution. Results and Discussion
Figure 1. Specific excess free energy and specific internal energy for lattice polymer of length L ) 25 (dot lines), L ) 50 (dashed lines), and L ) 100 (solid lines) via β obtained from expanded ensemble simulation. Marked points correspond to inverse temperatures βm in the expanded ensemble (Table 1).
run length was only 20 000 MC steps, while final runs consisted of 250 000 MC steps. Final results for chains of length L ) 25, 50, and 100 are presented in Figure 1. It is noted that the differences between the corresponding specific values are small. The value of βU/L passes through a maximum as a result of superposition of two opposite factorssgrowth of β and the fall of the average number of overlaps for growing β. At high β, the latter tendency prevails and the product tends to zero. So at β > 5 free energy βF/L practically equals -S/L. For a short chain the total number of conformations is, relatively, not very great so that thermodynamics of the system can be calculated exactly. Table 2 gives the total exact number of configurations nc for each number of overlaps N. So one can calculate exactly the free and internal energies for each βm:
βmFm ) -ln[∑nc(N) exp(-βmN)]
(12)
N
All simulations in the expanded ensemble were started with the described-above procedure for the optimization of balancing factors ηm. The results of optimization for L ) 25 are shown in Table 1; in other cases we similarly obtained values of balancing factors providing nearly equal probability distribution over subensembles after two to four iterations. Each preliminary
Um )
∑N Nnc(N) exp(-βmN) ∑N nc(N) exp(-βmN)
(13)
1156 J. Phys. Chem., Vol. 100, No. 4, 1996
Vorontsov-Velyaminov et al.
TABLE 2: Results of Simulation in the Expanded Ensemble Compared with Exact Results for L ) 8a MC result
exact result
βm
-βmFm
Um
-βmFm
Um
0.0 0.1 0.45 0.9 1.5 2.25 3.15 4.2 5.4 6.75 8.25
0 0.228 0.572 0.898 1.153 1.311 1.409 1.441 1.456 1.468 1.468
1.627 1.317 0.909 0.552 0.291 0.133 0.051 0.019 0.005 0.0011 0.0005
0 0.2244 0.5260 0.8944 1.1525 1.3175 1.4052 1.4444 1.4591 1.4638 1.4650
1.6601 1.3488 0.9363 0.5756 0.3128 0.1478 0.0601 0.0211 0.0063 0.0016 0.0003
a The length of the MC run was 100 000 MC steps. Exact data were obtained using (12) and (13), where the number of configurations for each number of overlaps nc(N) ) 64 661, 90 806, 63 225, 25 894, 21 832, 5733, 3535, 2836, 549, 574, 110, 120, 40, 20, 0, 0, 1) for N ) 0, ..., 16.
Figure 2. 1/L extrapolation of specific entropy for a self-avoiding walk on the simple cubic lattice.
TABLE 3: L Dependence of Excess Entropy -∆S ) βMFM - β0F0 (Where β0F0 ) -N ln(6) Is Entropy of the Ideal Chain) and Specific Excess Entropy -S/L ) λ for Lattice Polymers MC result
exact result
L
-∆S
λ
-∆S
λ
5 8 25 50 100
0.814 1.468 5.619 11.628 24.077
0.163 0.184 0.225 0.234 0.241
0.789 1.465
0.158 0.183
In Table 2 a comparison between exact results and that of the expanded ensemble simulations is shown. For the free energy one can see perfect agreement between themsdifferences are in the limit of the statistical error (0.005). Disagreements for internal energy are larger, but this has an explanation: 100 000 MC steps in the expanded ensemble correspond to about 10 000 MC steps in each subensemble, and a statistical error of 1-2% is natural for such simulations. Final values of excess free energies and corresponding specific entropies are shown in Table 3. These data make it possible to extrapolate results to large L when entropy should become an extensive quantity (SL/L should not depend on L). Figure 2 presents such an extrapolation, which gives an asymptotic value ∆SL = 0.2485. Existence of such a limit corresponds to the earlier observation in ref 25 that WSAW/W0 ) exp(-λL) where WSAW is a number of self-avoiding conformations, while W0 is the total number of free random walks for a chain of L monomers. In our terms ∆S ) SSAW - S0 ) ln(WSAW/W0) ) -λL so ∆SL/L ) -λ is the specific excess entropy. In ref 25 for the tetrahedral lattice was obtained λ ) 0.039. If we recalculate our value (0.2485) to another reference systemsrandom walk with forbidden back steps with W0 ) 6 × 5L-1 it gives the value 0.063 which is of the same order as that of ref 25 (unfortunately we have not performed our simulations directly for the tetrahedral lattice). Another important observation is presented in Figure 3. Straight dashed lines with different slopes from the origin are dependences of reduced Hamiltonian Hm ) βmU on reduced energy (the number of overlaps) N for various reciprocal temperatures βm (different m). Another series of straight solid lines are corresponding dependences of “shifted Hamiltonian” Hm - ηm obtained by parallel shift down of the initial lines with optimal values of balancing factors ηm for polymer length L ) 25 (Table 1). The shifted lines intersect with each other at different points, while the initial ones all pass through the
Figure 3. Reduced and shifted reduced Hamiltonians for different subensembles. Dashed lines are βmHm for different m (reduced Hamiltonians), solid lines are βmHm - ηm (reduced shifted Hamiltonians) for optimal values of balancing factors for L ) 25. Points of intersection of solid lines corresponding to mth and (m + 1)th subensembles are about in the middle between the values of internal energies Um and Um+1.
origin. Mean energies (mean number of overlaps) for different m are also presented in the U axis. It is clearly seen that mean energies and the points of the intersections of corresponding lines follow each other; e.g., the intersection point of the shifted lines with m ) 1 and 2 lies between 〈N(q)〉1 and 〈N(q)〉2. MC transitions between neigboring subensembles are especially easy in the intersection points and in their close vicinities. It makes it possible to follow all the paths between m ) 0 and 10 which might be impossible for a nonoptimal choice of ηm. Also important is the behavior of the mean values of energy Um ) 〈N(q)〉m if we shift them according to Um f Um - ηm/βm. Simple calculation shows that the range of shifted energies Um - ηm/βm is about 2, while that for nonshifted values, as one may see in Figure 3 it is about 7.5. So the distribution interval becomes almost 4 times more narrow, and energy distributions should strongly overlap, especially for neighboring subensembles. This feature was discussed earlier in ref 15. Continuous Polymer Model The expanded ensemble technique described above has been applied also to simulate polymer melt: a periodic cell filled with one or several chain molecules. We used the algorithm of the work in ref 20, combining molecular dynamics steps facilitating a more effective global update, with MC steps for transitions between the subensembles. The computations were performed for a single polymer chain of 32 monomers connected
Free Energy Calculations for Selected Polymers
J. Phys. Chem., Vol. 100, No. 4, 1996 1157
TABLE 4: Results of the Simulation for the Continuum Polymer Modela polymer with LJ interaction
ideal polymer
LJ system
λ ) Tmin/T
βF
βUbond
βULJ
S
βF
βUbond
S0
0.0 0.001 0.002 0.003 0.005 0.007 0.01 0.02 0.04 0.07 0.1 0.15 0.2 0.3 0.4 0.6 0.8 1.0
0 0.458 0.706 0.863 1.059 1.188 1.327 1.622 1.955 2.22 2.349 2.411 2.359 2.067 1.63 0.529 -0.737 -2.096
0 0.011 0.020 0.030 0.049 0.067 0.092 0.160 0.254 0.332 0.378 0.418 0.422 0.415 0.403 0.402 0.407 0.417
0.0 0.31 0.36 0.36 0.33 0.32 0.31 0.30 0.24 0.09 -0.09 -0.43 -0.79 -.155 -2.35 -4.00 -5.69 -7.38
0.0 -0.15 -0.34 -0.50 -0.72 -0.087 -1.02 -1.33 -1.71 -2.13 -2.44 -2.84 -3.15 -3.62 -3.98 -4.53 -4.95 -5.29
0.0 0.011 0.022
0.0 0.011 0.021
0.0 0.00 -0.001
0.053
0.051
-0.002
0.103 0.194 0.345 0.533 0.669 0.85 0.988 1.192 1.339 1.547 1.692 1.803
0.096 0.171 0.265 0.362 0.442 0.465 0.495 0.510 0.512 0.507 0.501 0.495
-0.007 -0.017 -0.08 -0.171 -0.247 -0.385 -0.493 -0.682 -0.827 -1.04 -1.191 -1.308
βF
βULJ
S - S0 - SLJ
0.0 0.45 0.68 0.84 1.00
0.0 0.31 0.36 0.36 0.34
0.0 -0.01 -0.02 -0.06
1.23
0.31
-0.09
1.64
0.25
-0.31
1.72
-0.07
-0.41
1.46 1.04 0.46 -0.78 -2.13 -3.56
-0.78 -1.51 -2.31 -3.88 -5.47 -7.25
-0.42 -0.39 -0.39 -0.40 -0.41 -0.39
a Compared are the results of simulations for the continuum polymer model (LJ and bond interactions) with ideal polymer (bond potential, no LJ) and the pure Lennard-Jones system (the latter data are taken from ref 20). In the last column the values of the entropy difference (S - S0 SLJ) are shown. All values are given per atom.
by a “soft” harmonic potential. Excluded volume effects were represented by the modified Lennard-Jones (LJ) potential (used in ref 20):
|
a b + U(r) ) r6 r12 A - Br
r > rlin r < rlin
(14)
where A and B were chosen to provide the smoothness of the resulting potential curve: A ) 7a/rlin6 + 13b/rlin12 and B ) 6a/ rlin7 + 12b/rlin13. For comparison of the results with simulations of the work of ref 20, performed for argon atoms, beads with the same mass and Lennard-Jones parameters where chosen for polymer simulation (σ ) (-b/a)1/6 ) 3.41 Å, ) a2/4b ) 119.6 K). Linearization of the potential at small distances r < rlin ) 0.72σ was made to provide easy transition from the LJ potential to the ideal gas. For minimal temperature in our expanded ensemble this potential is equivalent to the true Lennard-Jones potential since particles never come to each other closer than rlin. Computations were performed at relative density F* ) Nσ3/V ) 0.846. For imitation of covalent bonds a harmonic interaction of adjacent monomers was used:
Un,n+1 ) κ(r-r0)2/2
(15)
where equilibrium bond length was chosen as r0 ) 1.25σ, the force constant κ was chosen from κr02 ) 200kTmin, Tmin ) 86 K being the minimal temperature in the expanded ensemble (it corresponds to the reduced temperature T* ) Tmin/ ) 0.719). This condition for κ provides average bond length fluctuations of about 0.1r0. The selected potential parameters prevent the chain links from slipping through each other and allow use of the integration MD step ∆t ) 1 fs. The MD procedure used was a version of the conventional leapfrog velocity algorithm.26 The temperature constraint was provided by application of Berendsen’s “weak-coupling” scheme.27 The latter was chosen because it does not cause temperature oscillations intrinsic for the Nose scheme.28 When implementing MC steps with changing subensembles, we left the kinetic part of the Hamiltonian constant and multiplied the potential part by λ, changing from 0 to 1 (as it was done in ref 20). The case λ ) 1 corresponds to our original
system at T ) Tmin, while λ ) 0 corresponds to the ideal gas case; thus free energy differences computed via molecular modeling correspond to the excess free energies. While calculating averages over the configurations with fixed m, we obtain equilibrium results corresponding to a temperature Tm ) Tmin/λm. Simulational results presented in Table 4 show that our scheme combining MD and MC steps enables computation of free energy differences between polymers at different temperatures. In the course of relatively short runs (20-50 ps) we optimized balancing factors so that all subensembles had nearly equal probabilities; then we used a long run (2 ns) for the final averaging. We carried out also a simulation of an “ideal” polymer with a switched off Lennard-Jones interaction. The ideal continuum polymer is a chain of noninteracting beads connected by harmonic potential (15); thus it is an analogue of a random walk on the lattice. These data are also presented in Table 4. Free energy differences between polymer with LJ interaction and the ideal one show the effects of interaction between monomers (excluded volume effects). In Table 4 we also display data from ref 20 for a pure Lennard-Jones system with the same LJ parameters. One can do several observations. The LJ parts of the internal energy are practically equal through the whole range of λ. This means that the local structure of such a polymer melt is nearly the same as that for the Lennard-Jones system and may be a consequence so that we used the equilibrium distance r0 close to the minimum of the LJ potential (14). The bond energies in the interacting polymer are a little bit less than in the case of the ideal polymer; this means that the interaction between beads suppresses slightly bond fluctuations. Another interesting observation is the behavior of the entropies. The difference S - S0 is an entropy difference due to Lennard-Jones interaction between monomers (excluded volume effect) and in this sense an analogue of the entropy difference which we calculated for the lattice polymers. This difference is determined mainly by the Lennard-Jones interactions between beads. If we further subtract from this difference the entropy of a pure Lennard-Jones system, we obtain quantities listed in the last column of Table 4. One can see that the quantity (S - S0 - SLJ) practically does not depend on β in the
1158 J. Phys. Chem., Vol. 100, No. 4, 1996 range of β from 0.1 to 1 and is equal to about -0.4. This range of β corresponds to conditions when the polymer structure becomes well defined (βUbond reaches saturation). Conclusion In this work we presented calculations of free energy by the expanded ensemble method for two simple polymer models. The reference system in both cases has an exactly defined partition function. The procedures used are simple and provide accurate results. This is supported by independent calculations for polymers by Wilding and Muller.24 The suggested method can be easily applied to more complicated systems of polymer molecules such as closed chains, brushes, stars, polymer nets, copolymers, biopolymers, and polyelectrolytes. Other systems can be used as a reference. Acknowledgment. P.N.V.-V. expresses gratitude to Professors M. A. Coplan and H. L. Friedman for their interest, their patience in our long discussions, and their determinant help and continuous encouragement along the course of this contribution. References and Notes (1) Strastma, T. P.; Berendsen, H. J. C.; Postma, J. P. M. J. Chem. Phys. 1986, 85, 6720. (2) Pearlman, D. A. J. Phys. Chem. 1994, 98, 1487; J. Comput. Chem. 1994, 15, 105. (3) Mezei, M. Mol. Simul. 1993, 10, 225. (4) Widom, B. J. Chem. Phys. 1963, 39, 2804. (5) Adams, D. J. Mol. Phys. 1974, 28, 1241. (6) Valleau, J. P.; Card, D. N. J. Chem. Phys. 1972, 57, 5457. (7) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. Valleau, J. P. J. Comput. Phys. 1991, 96, 193. (8) Mezei, M. J. Comput. Phys. 1987, 68, 237.
Vorontsov-Velyaminov et al. (9) Pearlman, D. E.; Kollman, P. A. J. Chem. Phys. 1989, 91, 7831. Kollman, P. A. Chem. ReV. 1993, 93, 2395. (10) Reynolds, C. A.; King, P. M.; Richard, W. G. Mol. Phys. 1992, 76, 251. (11) Bennet, C. H. J. Comput. Phys. 1976, 22, 245. (12) Rittger, E. Mol. Phys. 1993, 79, 1073. (13) Frenkel, D. In Molecular dynamics simulation of statisticalmechanical systems; Ciccotti, G., Hoover, W. G., Eds.; North-Holland: Amsterdam, 1986; p. 151. Frenkel, D. Computer simulation in materials science; NATO ASI Series Applied Science; Dordrecht: Boston/London, 1991; Vol. 205, p 85. (14) Beveridge, D. L.; DiCapua, F. M. Annu. ReV. Biophys. Chem. 1989, 18, 431. (15) Lyubartsev, A. P.; Martsinovskii, A. A.; Shevkunov, S. V.; Vorontsov-Velyaminov, P. N. J. Chem. Phys. 1992, 96, 1776. (16) Shevkunov, S. V.; Martsinovskii, A. A.; Vorontsov-Velyaminov, P. N. High Temp. Phys. (USSR) 1988, 26, 246. Shevkunov, S. V.; Vorontsov-Velyaminov, P. N.; Martsinovskii, A. A. Mol. Simul. 1991, 5, 119. (17) Nezbeda, I.; Kolafa, J. A. Mol. Simul. 1991, 5, 391. (18) Marinari, E.; Parisi, G. Europhys. Lett. 1992, 19, 451. (19) Attard, P. J. Chem. Phys. 1993, 98, 2225. (20) Lyubartsev, A. P.; Laaksonen, A.; Vorontsov-Velyaminov, P. N. Mol. Phys. 1994, 82, 455. (21) Lyubartsev, A. P.; Martsinovskii, A. A.; Vorontsov-Velyaminov, P. N.; Kuznetsova, T. V. Russ. J. Phys. Chem. 1993, 66, 230. (22) Kuznetsova, T. V.; Voronsov-Velyaminov, P. N. J. Phys.: Condens. Matter 1993, 5, 717. (23) Lyubartsev, A. P.; Nordenskio¨ld, L. J. Phys. Chem. 1995, 99, 10373. (24) Wilding, N. B.; Muller, M. J. Chem. Phys. 1994, 101, 4324. (25) Wall, F. T.; Rubin, R. J.; Isaakson, L. M. J. Chem. Phys. 1957, 27, 186. (26) Allen, M. P.; Tildesley, D. J. The computer simulation of liquids; Clarendon: Oxford, U.K., 1987. (27) Berendsen, H. J. C. In Computer simulation in material science; Meyer, M., Pontikis, V., Eds.; NATO ASI Series Applied Science: Dordrecht: Boston/London, 1991; Vol. 205, p 139. (28) Nose, S. Molec. Phys. 1984, 52, 255.
JP951285I
