We propose a new generalized-ensemble algorithm, which we refer to as the multibaricmultithermal Monte Carlo method. The multibaric-multithermal Monte Carlo simulations perform random walks widely both in volume space and in potential energy space. From only one simulation run, one can calculate isobaric-isothermal-ensemble averages at any pressure and any temperature. We test the effectiveness of this algorithm by applying it to the Lennard-Jones 12-6 potential system with 500 particles. It is found that a single simulation of the new method indeed gives accurate average quantities in isobaric-isothermal ensemble for a wide range of pressure and temperature.PACS numbers: 64.70. Fx, 02.70.Ns, 47.55.Dz Monte Carlo (MC) algorithm is one of the most widely used methods of computational physics. In order to realize desired statistical ensembles, corresponding MC techniques have been proposed [1,2,3,4,5]. The first MC simulation was performed in the canonical ensemble in 1953 [1]. This method is called the Metropolis algorithm and widely used. The canonical probability distribution P N V T (E) for potential energy E is given by the product of the density of states n(E) and the Boltzmann weight factor e −β0E :where β 0 is the inverse of the product of the Boltzmann constant k B and temperature T 0 at which simulations are performed. Since n(E) is a rapidly increasing function and the Boltzmann factor decreases exponentially, P N V T (E) is a bell-shaped distribution. The isobaric-isothermal MC simulation [2] is also extensively used. This is because most experiments are carried out under the constant pressure and constant temperature conditions. Both potential energy E and volume V fluctuate in this ensemble. The distribution P N P T (E, V ) for E and V is given byHere, the density of states n(E, V ) is given as a function of both E and V , and H is the "enthalpy":where P 0 is the pressure at which simulations are performed. This ensemble has bell-shaped distributions in both E and V . Besides the above physical ensembles, it is now almost a default to simulate in artificial, generalized ensembles so that the multiple-minima problem, or the broken ergodicity problem, in complex systems can be overcome (for a recent review, see Ref.[6]). The multicanonical algorithm [7,8] is one of the most well known such methods in generalized ensemble. In multicanonical ensemble, a non-Boltzmann weight factor W mc (E) is used. This multicanonical weight factor is characterized by a flat probability distribution P mc (E):and thus a free random walk is realized in the potential energy space. This enables the simulation to escape from any local-minimum-energy state and to sample the configurational space more widely than the conventional canonical MC algorithm. Another advantage is that one can obtain various canonical ensemble averages at any temperature from a * Electronic address: hokumura@ims.ac.jp † Electronic address: okamotoy@ims.ac.jp