- Email: [email protected]
Molecular simulations in the multibaric-multithermal ensembles
Computer Physics Communications 169 (2005) 317–321 www.elsevier.com/locate/cpc
Molecular simulations in the multibaric-multithermal ensembles Hisashi Okumura a,b,∗ , Yuko Okamoto a,b a Department of Theoretical Studies, Institute for Molecular Science, Okazaki, Aichi 444-8585, Japan b Departments of Functional and Structural Molecular Science, The Graduate University for Advanced Studies,
Okazaki, Aichi 444-8585, Japan Available online 20 April 2005
Abstract We present new generalized-ensemble molecular dynamics and Monte Carlo simulation methods, which we refer to as the multibaric-multithermal algorithms. The multibaric-multithermal simulations perform random walks widely both in the potential-energy space and in the volume space. Thus, these algorithms allow the simulations to escape from any localminimum-energy state. In order to test the effectiveness of these algorithms, we investigate the liquid–gas phase transition of a Lennard–Jones 12–6 potential system by the multibaric-multithermal Monte Carlo simulation. 2005 Elsevier B.V. All rights reserved. Keywords: Multibaric-multithermal ensemble; Monte Carlo simulation; Molecular dynamics simulation; Lennard–Jones
1. Introduction Molecular dynamics (MD) algorithm and Monte Carlo (MC) algorithm are two of the most widely used methods of computational physics. In order to realize desired statistical ensembles, such as the canonical and isobaric-isothermal ensembles, corresponding MD and MC techniques have been proposed [1–7]. Besides the above physical ensembles, it is also common to simulate in artificial, generalized ensembles so that the multiple-minima problem in complex systems can be overcome. The multicanonical algorithm [8– 11] is one of the most well-known such methods in * Corresponding author.
E-mail address: [email protected] (H. Okumura). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2005.03.071
generalized ensemble. In the multicanonical ensemble, a non-Boltzmann weight factor is used so that a free one-dimensional random walk is realized in the potential-energy space. Thus, one can obtain various canonical-ensemble averages at any temperature from a single simulation run by the reweighting techniques [12]. However, it is impossible to change volume V and to specify pressure P as in experimental environments, because the simulation is performed in a fixed volume. In order to overcome these difficulties, we recently proposed new algorithms, the multibaric-multithermal ensemble MC [13,14] and MD [15] algorithms, in which one can obtain various isobaric-isothermal ensembles from only one simulation run. This is an extension of the multicanonical algorithm to realize
318
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
two-dimensional random walks both in the potentialenergy space and in the volume space. In this article, we describe our algorithms in this ensemble and apply these algorithms to investigate the liquid–gas phase transition of a Lennard–Jones 12–6 system [16].
Hmbt =
2. Methods
i
In the multibaric-multithermal ensemble [13,14], every state is sampled with a weight factor Wmbt (E, V ) ≡ exp −β0 Hmbt (E, V ) (1) so that a uniform distribution for the potential energy E and volume V may be obtained: Pmbt (E, V ) = n(E, V )Wmbt (E, V ) = constant.
we employ the reweighting techniques [12] for the results of the production run. In order to obtain the multibaric-multithermal ensemble by an MD simulation, the Hamiltonian Hmbt is given by
(2)
Here, Wmbt (E, V ) and Hmbt are referred to as the multibaric-multithermal weight factor and the multibaric-multithermal enthalpy, respectively, and β0 is the inverse of the product of the Boltzmann constant kB and absolute temperature T0 at which simulations are performed. In order to perform a multibaric-multithermal MC simulation, we perform Metropolis sampling on the scaled coordinates q i = L−1 r i (r i are the real coordinates) and the volume V (here, the particles √ are placed in a cubic box of a side of size L ≡ 3 V ). The trial moves of q i and V are generated in the same way as in the isobaric-isothermal MC simulation [3,4]. The multibaric-multithermal enthalpy ≡ is changed from Hmbt (E(q (N ) , V ), V ) to Hmbt (N ) , V ), V ) by these trial moves. The trial Hmbt (E(q moves will now be accepted with the probability w(o → n) = min 1, exp −β0 Hmbt − Hmbt − N kB T0 ln(V /V ) . (3) Instead of volume-space sampling in Eq. (3), the density-space sampling is also possible by w(o → n) = min 1, exp − Hmbt − Hmbt + (N + 2)kB T0 log(ρ /ρ) /kB T0 , (4) with the trail moves of number density to ρ from ρ, where ρ ≡ N/V . The multibaric-multithermal probability distribution is obtained by this scheme. In order to calculate the isobaric-isothermal-ensemble average,
π i 2 PV 2 1/3 + E(V q) + + P0 V 2M 2mi V 2/3 s 2
Ps2 + gkB T0 log s + δH E(V 1/3 q), V . 2Q (5) This Hamiltonian is based on the Nosé–Andersen algorithm [5–7] for the isobaric-isothermal MD simulation. The variables π i , PV , and Ps are the conjugate momenta for q i , V , and Nosé’s additional degree of freedom s, respectively. The real time interval dt is associated with the virtual time interval dt by the relation dt = dt /s. The constant mi is the mass of particle i and M and Q are the artificial “mass” associated with V and s, respectively. The constant g corresponds to the number of degrees of freedom. In the case of an atomic system, g is equal to 3N (N is the number of particles) when the time development is performed in the real time t and g is equal to 3N + 1 when the time development is performed in the virtual time t . The last term δH (E, V ) is characteristic of the multibaricmultithermal MD simulation. This is the difference between multibaric-multithermal enthalpy Hmbt and “isobaric-isothermal enthalpy” H ≡ E + P0 V : +
Hmbt (E, V ) = H + δH (E, V ).
(6)
The equations of motion for q i , V , and s are obtained from the Hamiltonian in Eq. (5). The equations of motion in the real time development are ∂ δH ∂E −1 s˙ 2V˙ 1+ + q˙i , (7) − q¨i = ∂E ∂r i s 3V mi V 1/3
s2 1 mi V 2/3 q˙i 2 V¨ = M 3V i ∂ δH ∂E − 1+ ri · ∂E ∂r i i s2 s˙ ∂δH − (8) P0 + + V˙ , M ∂V s
s˙ 2 s 2/3 2 mi V q˙i − 3N kB T0 + . s¨ = (9) Q s i
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
319
Fig. 1. The probability distributions PN P T (E ∗ /N, ρ ∗ ; T0∗ , P0∗ ) from the isobaric-isothermal simulations. (a) The initial condition is the liquid phase. (b) The initial condition is the gas phase.
Performing MD simulations by these equations of motion, the multibaric-multithermal ensemble distribution Pmbt (E, V ) is obtained. The multibaric-multithermal weight factor is, however, not a priori known and has to be determined by the iterations of short MC or MD simulations [13–18]. After an optimal weight factor Wmbt (E, V ) is determined, a long production simulation is performed for data collection.
3. Liquid–gas phase transition of a Lennard–Jones 12–6 potential system In order to examine the effectiveness of our methods, we investigate the liquid–gas phase transition by the multibaric-multithermal MC simulation. We considered a Lennard–Jones 12–6 potential system of 500 particles (N = 500) in a cubic unit cell with periodic boundary conditions. The length and the energy are scaled in units of the Lennard–Jones diameter σ and the minimum value of the potential , respectively. We use an asterisk (∗) for the reduced quantities. We performed the multibaric-multithermal simulation at T0∗ = 1.1 and at P0∗ = 0.0455. These values of temperature and pressure are respectively lower than the critical temperature Tc∗ and the critical pressure Pc∗ [19–24]. Recent reliable data are Tc∗ = 1.3207(4) and Pc∗ = 0.1288(5) [19]. The pressure value of P0∗ = 0.0455 is the equilibrium vapor-pressure at T0∗ = 1.1. The cutoff radius rc∗ was taken to be L∗ /2. A cutoff correction was added for the pressure and the potential energy. In one MC sweep we made the trial moves of all particle coordinates and the number density (N + 1 trial moves altogether). For each trial move
the Metropolis evaluation of Eq. (4) was made. In order to obtain the weight factor Wmbt (E, ρ), a generalized version of the energy landscape paving method [16,25] was employed here. The production run was performed for 1 × 107 MC sweeps. In order to compare the multibaric-multithermal MC method with the conventional one, we also performed the isobaric-isothermal MC simulations. Two simulations are carried out from different initial conditions: One is the liquid state of ρ ∗ = 0.659 and the other is the gas state of ρ ∗ = 0.060. Fig. 1 shows the probability distribution PN P T (E ∗ /N, ρ ∗ ) from the isobaric-isothermal simulation. The simulation started from the liquid state has a single peak in the liquid phase (Fig. 1(a)) and the simulation started from the gas state has a single peak in the gas phase (Fig. 1(b)). These results mean that the isobaric-isothermal simulation samples only one phase. On the other hand, Fig. 2(a) shows that the multibaric-multithermal simulation has a broad and flat distribution Pmbt (E ∗ /N, ρ ∗ ) both in the potentialenergy space and in the density space, that is, it covers the wide ranges of E ∗ /N and ρ ∗ . We employed the reweighting techniques for the calculation of PN P T (E ∗ /N, ρ ∗ ). In order to satisfy the phase-coexistence conditions, the peak volumes of PN P T (E, ρ; T , P ) should be equal between the liquid phase and the gas phase. Fig. 2(b) shows the probability distribution PN P T (E ∗ /N, ρ ∗ ) reweighted at T ∗ = T0∗ = 1.1. Both liquid phase and gas phase coexist corresponding to the isobaric-isothermal peak in the liquid phase (Fig. 1(a)) and that in the gas phase (Fig. 1(b)), respectively.
320
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
Fig. 2. (a) The probability distribution Pmbt (E ∗ /N, ρ ∗ ) from the multibaric-multithermal simulation and (b) the reweighted distribution PN P T (E ∗ /N, ρ ∗ ; T ∗ , P ∗ ) to satisfy the coexistence conditions at T ∗ = T0∗ = 1.1.
Fig. 3. The phase coexistence curve in T ∗ –ρ ∗ plane. •: The multibaric-multithermal simulation data, ◦: the critical point estimated from the multibaric-multithermal simulation data of •, 2: the N V T plus test particle method [19], : the N P T plus test particle method by Okumura and Yonezawa [20], and : the N P T plus test particle method by Lotfi et al. [21].
multithermal simulation performed random walks in wide ranges of the potential energy and the number density, thus, it sampled both liquid phase and gas phase. On the other hand, the conventional isobaricisothermal simulation sampled only one phase whether the liquid phase or the gas phase depending on the initial condition. We obtained the critical constants as Tc∗ = 1.309 ± 0.013, ρc∗ = 0.317 ± 0.004, and Pc∗ = 0.124 ± 0.007. Our coexistence data agreed well with those determined previously by other methods. The multibaric-multithermal algorithms can be utilized for studying not only the liquid–gas phase transition but also more pressure-induced phenomena such as the pressure denaturation of proteins.
References Fig. 3 shows the liquid–gas coexistence curve in T ∗ –ρ ∗ plane obtained by the reweighting techniques. The critical constants are estimated as Tc∗ = 1.309 ± 0.013, ρc∗ = 0.317 ± 0.04, and Pc∗ = 0.124 ± 0.007. Fig. 3 also shows the coexistence data by the N V T plus test particle method [19] and the N P T plus test particle method [20,21]. Our coexistence data agreed well with those determined previously by other methods [19–24].
4. Conclusions We presented new MC and MD algorithms in the multibaric-multithermal ensemble. We successfully investigated the liquid–gas phase transition of a Lennard–Jones 12–6 potential system by the multibaric-multithermal MC simulation. The multibaric-
[1] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987, p. 110. [2] D. Frenkel, B. Smit, Understanding Molecular Simulation From Algorithms to Applications, Academic Press, San Diego, 2002, p. 111. [3] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [4] I.R. McDonald, Mol. Phys. 23 (1972) 41. [5] S. Nosé, Mol. Phys. 52 (1984) 255. [6] S. Nosé, J. Chem. Phys. 81 (1984) 511. [7] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [8] B.A. Berg, T. Neuhaus, Phys. Lett. B 267 (1991) 249. [9] B.A. Berg, T. Neuhaus, Phys. Rev. Lett. 68 (1992) 9. [10] U.H.E. Hansmann, Y. Okamoto, F. Eisenmenger, Chem. Phys. Lett. 259 (1996) 321. [11] N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B 101 (1997) 817. [12] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635; Phys. Rev. Lett. 63 (1989) 1658(E). [13] H. Okumura, Y. Okamoto, Chem. Phys. Lett. 383 (2004) 391.
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
[14] [15] [16] [17] [18]
H. Okumura, Y. Okamoto, Phys. Rev. E 70 (2004) 026702. H. Okumura, Y. Okamoto, Chem. Phys. Lett. 391 (2004) 248. H. Okumura, Y. Okamoto, J. Phys. Soc. Jpn. 73 (2004) 3304. B.A. Berg, T. Celik, Phys. Rev. Lett. 69 (1992) 2292. Y. Okamoto, U.H.E. Hansmann, J. Phys. Chem. 99 (1995) 11276. [19] H. Okumura, F. Yonezawa, J. Phys. Soc. Jpn. 70 (2001) 1990. [20] H. Okumura, F. Yonezawa, J. Chem. Phys. 113 (2000) 9162.
321
[21] A. Lotfi, J. Vrabec, J. Fischer, Mol. Phys. 76 (1992) 1319. [22] J.M. Caillol, J. Chem. Phys. 109 (1998) 4885. [23] J.J. Potoff, A.Z. Panagiotopoulos, J. Chem. Phys. 109 (1998) 10914. [24] A.Z. Panagiotopoulos, J. Phys.: Condens. Matter 12 (2000) R25. [25] U.H.E. Hansmann, L.T. Wille, Phys. Rev. Lett. 88 (2002) 068105.
Molecular simulations in the multibaric-multithermal ensembles Hisashi Okumura a,b,∗ , Yuko Okamoto a,b a Department of Theoretical Studies, Institute for Molecular Science, Okazaki, Aichi 444-8585, Japan b Departments of Functional and Structural Molecular Science, The Graduate University for Advanced Studies,
Okazaki, Aichi 444-8585, Japan Available online 20 April 2005
Abstract We present new generalized-ensemble molecular dynamics and Monte Carlo simulation methods, which we refer to as the multibaric-multithermal algorithms. The multibaric-multithermal simulations perform random walks widely both in the potential-energy space and in the volume space. Thus, these algorithms allow the simulations to escape from any localminimum-energy state. In order to test the effectiveness of these algorithms, we investigate the liquid–gas phase transition of a Lennard–Jones 12–6 potential system by the multibaric-multithermal Monte Carlo simulation. 2005 Elsevier B.V. All rights reserved. Keywords: Multibaric-multithermal ensemble; Monte Carlo simulation; Molecular dynamics simulation; Lennard–Jones
1. Introduction Molecular dynamics (MD) algorithm and Monte Carlo (MC) algorithm are two of the most widely used methods of computational physics. In order to realize desired statistical ensembles, such as the canonical and isobaric-isothermal ensembles, corresponding MD and MC techniques have been proposed [1–7]. Besides the above physical ensembles, it is also common to simulate in artificial, generalized ensembles so that the multiple-minima problem in complex systems can be overcome. The multicanonical algorithm [8– 11] is one of the most well-known such methods in * Corresponding author.
E-mail address: [email protected] (H. Okumura). 0010-4655/$ – see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2005.03.071
generalized ensemble. In the multicanonical ensemble, a non-Boltzmann weight factor is used so that a free one-dimensional random walk is realized in the potential-energy space. Thus, one can obtain various canonical-ensemble averages at any temperature from a single simulation run by the reweighting techniques [12]. However, it is impossible to change volume V and to specify pressure P as in experimental environments, because the simulation is performed in a fixed volume. In order to overcome these difficulties, we recently proposed new algorithms, the multibaric-multithermal ensemble MC [13,14] and MD [15] algorithms, in which one can obtain various isobaric-isothermal ensembles from only one simulation run. This is an extension of the multicanonical algorithm to realize
318
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
two-dimensional random walks both in the potentialenergy space and in the volume space. In this article, we describe our algorithms in this ensemble and apply these algorithms to investigate the liquid–gas phase transition of a Lennard–Jones 12–6 system [16].
Hmbt =
2. Methods
i
In the multibaric-multithermal ensemble [13,14], every state is sampled with a weight factor Wmbt (E, V ) ≡ exp −β0 Hmbt (E, V ) (1) so that a uniform distribution for the potential energy E and volume V may be obtained: Pmbt (E, V ) = n(E, V )Wmbt (E, V ) = constant.
we employ the reweighting techniques [12] for the results of the production run. In order to obtain the multibaric-multithermal ensemble by an MD simulation, the Hamiltonian Hmbt is given by
(2)
Here, Wmbt (E, V ) and Hmbt are referred to as the multibaric-multithermal weight factor and the multibaric-multithermal enthalpy, respectively, and β0 is the inverse of the product of the Boltzmann constant kB and absolute temperature T0 at which simulations are performed. In order to perform a multibaric-multithermal MC simulation, we perform Metropolis sampling on the scaled coordinates q i = L−1 r i (r i are the real coordinates) and the volume V (here, the particles √ are placed in a cubic box of a side of size L ≡ 3 V ). The trial moves of q i and V are generated in the same way as in the isobaric-isothermal MC simulation [3,4]. The multibaric-multithermal enthalpy ≡ is changed from Hmbt (E(q (N ) , V ), V ) to Hmbt (N ) , V ), V ) by these trial moves. The trial Hmbt (E(q moves will now be accepted with the probability w(o → n) = min 1, exp −β0 Hmbt − Hmbt − N kB T0 ln(V /V ) . (3) Instead of volume-space sampling in Eq. (3), the density-space sampling is also possible by w(o → n) = min 1, exp − Hmbt − Hmbt + (N + 2)kB T0 log(ρ /ρ) /kB T0 , (4) with the trail moves of number density to ρ from ρ, where ρ ≡ N/V . The multibaric-multithermal probability distribution is obtained by this scheme. In order to calculate the isobaric-isothermal-ensemble average,
π i 2 PV 2 1/3 + E(V q) + + P0 V 2M 2mi V 2/3 s 2
Ps2 + gkB T0 log s + δH E(V 1/3 q), V . 2Q (5) This Hamiltonian is based on the Nosé–Andersen algorithm [5–7] for the isobaric-isothermal MD simulation. The variables π i , PV , and Ps are the conjugate momenta for q i , V , and Nosé’s additional degree of freedom s, respectively. The real time interval dt is associated with the virtual time interval dt by the relation dt = dt /s. The constant mi is the mass of particle i and M and Q are the artificial “mass” associated with V and s, respectively. The constant g corresponds to the number of degrees of freedom. In the case of an atomic system, g is equal to 3N (N is the number of particles) when the time development is performed in the real time t and g is equal to 3N + 1 when the time development is performed in the virtual time t . The last term δH (E, V ) is characteristic of the multibaricmultithermal MD simulation. This is the difference between multibaric-multithermal enthalpy Hmbt and “isobaric-isothermal enthalpy” H ≡ E + P0 V : +
Hmbt (E, V ) = H + δH (E, V ).
(6)
The equations of motion for q i , V , and s are obtained from the Hamiltonian in Eq. (5). The equations of motion in the real time development are ∂ δH ∂E −1 s˙ 2V˙ 1+ + q˙i , (7) − q¨i = ∂E ∂r i s 3V mi V 1/3
s2 1 mi V 2/3 q˙i 2 V¨ = M 3V i ∂ δH ∂E − 1+ ri · ∂E ∂r i i s2 s˙ ∂δH − (8) P0 + + V˙ , M ∂V s
s˙ 2 s 2/3 2 mi V q˙i − 3N kB T0 + . s¨ = (9) Q s i
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
319
Fig. 1. The probability distributions PN P T (E ∗ /N, ρ ∗ ; T0∗ , P0∗ ) from the isobaric-isothermal simulations. (a) The initial condition is the liquid phase. (b) The initial condition is the gas phase.
Performing MD simulations by these equations of motion, the multibaric-multithermal ensemble distribution Pmbt (E, V ) is obtained. The multibaric-multithermal weight factor is, however, not a priori known and has to be determined by the iterations of short MC or MD simulations [13–18]. After an optimal weight factor Wmbt (E, V ) is determined, a long production simulation is performed for data collection.
3. Liquid–gas phase transition of a Lennard–Jones 12–6 potential system In order to examine the effectiveness of our methods, we investigate the liquid–gas phase transition by the multibaric-multithermal MC simulation. We considered a Lennard–Jones 12–6 potential system of 500 particles (N = 500) in a cubic unit cell with periodic boundary conditions. The length and the energy are scaled in units of the Lennard–Jones diameter σ and the minimum value of the potential , respectively. We use an asterisk (∗) for the reduced quantities. We performed the multibaric-multithermal simulation at T0∗ = 1.1 and at P0∗ = 0.0455. These values of temperature and pressure are respectively lower than the critical temperature Tc∗ and the critical pressure Pc∗ [19–24]. Recent reliable data are Tc∗ = 1.3207(4) and Pc∗ = 0.1288(5) [19]. The pressure value of P0∗ = 0.0455 is the equilibrium vapor-pressure at T0∗ = 1.1. The cutoff radius rc∗ was taken to be L∗ /2. A cutoff correction was added for the pressure and the potential energy. In one MC sweep we made the trial moves of all particle coordinates and the number density (N + 1 trial moves altogether). For each trial move
the Metropolis evaluation of Eq. (4) was made. In order to obtain the weight factor Wmbt (E, ρ), a generalized version of the energy landscape paving method [16,25] was employed here. The production run was performed for 1 × 107 MC sweeps. In order to compare the multibaric-multithermal MC method with the conventional one, we also performed the isobaric-isothermal MC simulations. Two simulations are carried out from different initial conditions: One is the liquid state of ρ ∗ = 0.659 and the other is the gas state of ρ ∗ = 0.060. Fig. 1 shows the probability distribution PN P T (E ∗ /N, ρ ∗ ) from the isobaric-isothermal simulation. The simulation started from the liquid state has a single peak in the liquid phase (Fig. 1(a)) and the simulation started from the gas state has a single peak in the gas phase (Fig. 1(b)). These results mean that the isobaric-isothermal simulation samples only one phase. On the other hand, Fig. 2(a) shows that the multibaric-multithermal simulation has a broad and flat distribution Pmbt (E ∗ /N, ρ ∗ ) both in the potentialenergy space and in the density space, that is, it covers the wide ranges of E ∗ /N and ρ ∗ . We employed the reweighting techniques for the calculation of PN P T (E ∗ /N, ρ ∗ ). In order to satisfy the phase-coexistence conditions, the peak volumes of PN P T (E, ρ; T , P ) should be equal between the liquid phase and the gas phase. Fig. 2(b) shows the probability distribution PN P T (E ∗ /N, ρ ∗ ) reweighted at T ∗ = T0∗ = 1.1. Both liquid phase and gas phase coexist corresponding to the isobaric-isothermal peak in the liquid phase (Fig. 1(a)) and that in the gas phase (Fig. 1(b)), respectively.
320
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
Fig. 2. (a) The probability distribution Pmbt (E ∗ /N, ρ ∗ ) from the multibaric-multithermal simulation and (b) the reweighted distribution PN P T (E ∗ /N, ρ ∗ ; T ∗ , P ∗ ) to satisfy the coexistence conditions at T ∗ = T0∗ = 1.1.
Fig. 3. The phase coexistence curve in T ∗ –ρ ∗ plane. •: The multibaric-multithermal simulation data, ◦: the critical point estimated from the multibaric-multithermal simulation data of •, 2: the N V T plus test particle method [19], : the N P T plus test particle method by Okumura and Yonezawa [20], and : the N P T plus test particle method by Lotfi et al. [21].
multithermal simulation performed random walks in wide ranges of the potential energy and the number density, thus, it sampled both liquid phase and gas phase. On the other hand, the conventional isobaricisothermal simulation sampled only one phase whether the liquid phase or the gas phase depending on the initial condition. We obtained the critical constants as Tc∗ = 1.309 ± 0.013, ρc∗ = 0.317 ± 0.004, and Pc∗ = 0.124 ± 0.007. Our coexistence data agreed well with those determined previously by other methods. The multibaric-multithermal algorithms can be utilized for studying not only the liquid–gas phase transition but also more pressure-induced phenomena such as the pressure denaturation of proteins.
References Fig. 3 shows the liquid–gas coexistence curve in T ∗ –ρ ∗ plane obtained by the reweighting techniques. The critical constants are estimated as Tc∗ = 1.309 ± 0.013, ρc∗ = 0.317 ± 0.04, and Pc∗ = 0.124 ± 0.007. Fig. 3 also shows the coexistence data by the N V T plus test particle method [19] and the N P T plus test particle method [20,21]. Our coexistence data agreed well with those determined previously by other methods [19–24].
4. Conclusions We presented new MC and MD algorithms in the multibaric-multithermal ensemble. We successfully investigated the liquid–gas phase transition of a Lennard–Jones 12–6 potential system by the multibaric-multithermal MC simulation. The multibaric-
[1] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987, p. 110. [2] D. Frenkel, B. Smit, Understanding Molecular Simulation From Algorithms to Applications, Academic Press, San Diego, 2002, p. 111. [3] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [4] I.R. McDonald, Mol. Phys. 23 (1972) 41. [5] S. Nosé, Mol. Phys. 52 (1984) 255. [6] S. Nosé, J. Chem. Phys. 81 (1984) 511. [7] H.C. Andersen, J. Chem. Phys. 72 (1980) 2384. [8] B.A. Berg, T. Neuhaus, Phys. Lett. B 267 (1991) 249. [9] B.A. Berg, T. Neuhaus, Phys. Rev. Lett. 68 (1992) 9. [10] U.H.E. Hansmann, Y. Okamoto, F. Eisenmenger, Chem. Phys. Lett. 259 (1996) 321. [11] N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B 101 (1997) 817. [12] A.M. Ferrenberg, R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635; Phys. Rev. Lett. 63 (1989) 1658(E). [13] H. Okumura, Y. Okamoto, Chem. Phys. Lett. 383 (2004) 391.
H. Okumura, Y. Okamoto / Computer Physics Communications 169 (2005) 317–321
[14] [15] [16] [17] [18]
H. Okumura, Y. Okamoto, Phys. Rev. E 70 (2004) 026702. H. Okumura, Y. Okamoto, Chem. Phys. Lett. 391 (2004) 248. H. Okumura, Y. Okamoto, J. Phys. Soc. Jpn. 73 (2004) 3304. B.A. Berg, T. Celik, Phys. Rev. Lett. 69 (1992) 2292. Y. Okamoto, U.H.E. Hansmann, J. Phys. Chem. 99 (1995) 11276. [19] H. Okumura, F. Yonezawa, J. Phys. Soc. Jpn. 70 (2001) 1990. [20] H. Okumura, F. Yonezawa, J. Chem. Phys. 113 (2000) 9162.
321
[21] A. Lotfi, J. Vrabec, J. Fischer, Mol. Phys. 76 (1992) 1319. [22] J.M. Caillol, J. Chem. Phys. 109 (1998) 4885. [23] J.J. Potoff, A.Z. Panagiotopoulos, J. Chem. Phys. 109 (1998) 10914. [24] A.Z. Panagiotopoulos, J. Phys.: Condens. Matter 12 (2000) R25. [25] U.H.E. Hansmann, L.T. Wille, Phys. Rev. Lett. 88 (2002) 068105.