Analysis of a two-phase (phase separated) simulation

In this previous example, we analyzed a simulation conducted above the critical temperature. In this tutorial, we now consider an LJ simulation with potentially two phases at \(T^*=1.2\).

This tutorial uses data from the SRSW. For an example of running and analyzing a phase separated system, see tutorial 10 in the flat_histogram plugin.

import numpy as np
import pandas as pd
from pyfeasst import macrostate_distribution

srsw = macrostate_distribution.MacrostateDistribution(file_name="../test/data/stat120.csv",

# Set show_plot to True to see the plots. Plots take up too much space in the git repo.
show_plot = False
if show_plot:
    %matplotlib inline
    import matplotlib.pyplot as plt

Note that this macrostate distribution may contain a second maximum depending on the value of the chemical potential. But at the chemical potential shown above, the low density phase is unstable and there is no minimum.

assert len(srsw.minimums()) == 0

However, a minima appears after reweighting to a lower chemical potential.

srsw_rw = srsw.reweight(delta_beta_mu=-0.1)
assert len(minima) == 1 and minima.values[0] == 152

if show_plot:
    plt.gca().axvline(minima.values[0], color='black')
    plt.text(20, -40, 'low density phase')
    plt.text(200, -40, 'high density phase\n(more probable)')

Find the chemical potential at phase equilibrium subject to the constraint that the probability of observing the two phases are equal.

delta_beta_mu_equilibrium = srsw.equilibrium()
assert len(mins)==1

if show_plot:
    plt.gca().axvline(mn, color='black')
    plt.text(20, -40, 'saturated vapor')
    plt.text(225, -40, 'saturated liquid')

Compare chemical potential and equilibrium properties with published results.

assert len(mins) == 1 and mins.values[0] == 166
assert abs(-0.128071 - delta_beta_mu_equilibrium) < 0.01

volume = 8**3
vapor, liquid = srsw.split()
betaPV = -vapor.ln_prob()[0]
print('pressure:', pressure)
assert abs(0.07722557238264459 - pressure) < 1e-8

vapor_density = vapor.average_macrostate()/volume
print('vapor_density:', vapor_density)
assert abs(0.10035100340092683 - vapor_density) < 1e-8

liquid_density = liquid.average_macrostate()/volume
print('liquid_density:', liquid_density)
assert abs(0.5631867324202517 - liquid_density) < 1e-8
pressure: 0.07722557238264459
vapor_density: 0.10035100340092695
liquid_density: 0.563186732420252

Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!