Metropolis Monte Carlo Phase Transition for small systems Metropolis Monte Carlo simulation for the argon liquid fails to reproduce correctly the pressure-temperature isotherm from the point of liquid-vapour coexistence:

Why such a phenomenon appears ? The reason invoked from the book is that small systems (here 200 particles are simulated) do not form different phases "due to the relative important free energy cost associated with liquid-vapour interphase". 
However, as the system obeys periodic boundary conditions, the system simulated is supposed to be infinite..
Thank you for reading!
Reference: Frenkel and Smit, Understanding molecular simulation from Algorithms to Applications
After the comment of LonelyProf: here are some interesting surface shapes 

 A: The book is correct. 
On a pressure-density isotherm, in the constant-volume ensemble (i.e. $NVT$ fixed, measuring $P$ as a function of $\rho=N/V$) there will be a range of densities near the liquid-vapour (first-order) phase transition which corresponds to the two-phase region. In the thermodynamic limit this gives a horizontal line on the plot, at the coexistence pressure $P_c$, stretching between the equilibrium coexisting densities $\rho_{\ell}$ (liquid) and $\rho_v$ (vapour). As the overall density is changed, all that happens is that the proportion of the two phases changes. There is an interface between the two phases, but it doesn't contribute anything to the thermodynamics in this limit.
For a small system, as studied in computer simulation, the interface does contribute. Even in a periodic system, this is the case, because the system of two phases (with an interface) is duplicated in every periodic box. All possible configurations of liquid, vapour, and interface, give a contribution, but the ones with smallest interfacial area (and hence, smallest interfacial energy) contribute most: in practice this means the ones with a spherical droplet of one phase in the other, or with a planar slab of each phase, stretching all the way across the box. It may be shown (by calculating the free energy of the system as a function of density from a simple model, including an approximate estimate of the interfacial contribution) that this gives an apparent pressure-density curve of the kind shown in your diagram: the shape of an "S" on its side, sometimes called a "small system loop". 
If you simulate larger and larger systems, the loop becomes flatter and flatter, and you eventually reach the horizontal straight line expected in the thermodynamic limit.
This shape is often confused with the van der Waals loop which is part of the van der Waals equation of state. But that loop arises because of the approximate (and unphysical) nature of the equation of state. It is not the same thing.
If you were to simulate the same system in the constant pressure ensemble ($NPT$ specified, measuring $\rho$ as a function of $P$) you would not see the same loop. Instead, you would get a probability distribution of densities, which has two peaks in the region of the phase transition, one around $\rho_{\ell}$ and the other around $\rho_v$. The relative weight of the two peaks would reflect the amount of each phase present. The minimum between the peaks would correspond to the less likely configurations, which have an interface present between the two phases. There would be an extremely small range of pressures over which the average density changed over from $\rho_v$ to $\rho_{\ell}$; in other words, in this ensemble, the measured curve would look more like a slightly smoothed out version of the straight line seen in the thermodynamic limit.
This relation between the observed behaviour in the two ensembles goes back 50 years, to a paper by W. W. Wood, J Chem Phys, 48, 415 (1968). I'll include a bit of the derivation here. The probability distribution for volume in the constant-$NPT$ ensemble at pressure $P_c$ is
$$
\mathcal{P}(V) \propto Q(N,V,T) \exp(-P_cV/k_BT)
$$
where $Q$ is the canonical (constant-$NVT$) ensemble partition function. This follows from the basic expression of the constant-$NPT$ probability density in terms of particle coordinates and volume. Taking logs, and remembering that $Q$ is related to the Helmholtz free energy by $F(N,V,T)=-k_BT\ln Q(N,V,T)$,
$$
\ln \mathcal{P}(V) = \text{constant} - \frac{F(N,V,T)}{k_BT} - \frac{P_c}{k_BT}V
$$
and differentiating with respect to volume
$$
k_BT\frac{\partial \ln \mathcal{P}(V)}{\partial V} = P(N,V,T) - P_c .
$$
Here $P(N,V,T)$ is the pressure that would be measured in a constant-$NVT$ simulation, in the vicinity of $P_c$. It follows that, if $\ln\mathcal{P}(V)$ (measured at constant pressure) is a function that has two maxima separated by a minimum, then the curve $P(N,V,T)$ measured as a function of $V$ in a set of constant-$NVT$ simulations will cross the line $P=P_c$ at three points, giving a horizontal-S loop form. The same applies when we do the plot as a function of $\rho$ rather than $V$.
