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",
                                        macrostate_header="N",
                                        ln_prob_header="lnPI")

# 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
    plt.plot(srsw.ln_prob())

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)
minima=srsw_rw.minimums()
assert len(minima) == 1 and minima.values[0] == 152

if show_plot:
    plt.plot(srsw_rw.ln_prob())
    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)')
    plt.xlabel('N')
    plt.ylabel(r'$\ln\Pi$')

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()
mins=srsw.minimums()
assert len(mins)==1
mn=mins.values[0]

if show_plot:
    plt.plot(srsw.ln_prob())
    plt.gca().axvline(mn, color='black')
    plt.text(20, -40, 'saturated vapor')
    plt.text(225, -40, 'saturated liquid')
    plt.xlabel('N')
    plt.ylabel(r'$\ln\Pi$')

Compare chemical potential and equilibrium properties with published results.

https://www.nist.gov/mml/csd/chemical-informatics-research-group/sat-tmmc-liquid-vapor-coexistence-properties-long-range

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]
beta=1/1.2
pressure=betaPV/beta/volume
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!