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\).

import numpy as np
import pandas as pd
import feasst as fst
import pyfeasst

# load the SRSW data from file
dfsrsw = pd.read_csv(fst.install_dir() + "/plugin/flat_histogram/test/data/stat120.csv")
ln_prob = fst.LnProbability(fst.DoubleVector(dfsrsw["lnPI"]))
beta = 1./1.2

# plot macrostate and energy
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
ax = dfsrsw.plot('N', 'lnPI')

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

assert(len(ln_prob.minima()) == 0)

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

gce = fst.GrandCanonicalEnsemble(
    fst.Histogram(fst.args({"width": "1", "max": str(dfsrsw["N"].iloc[-1])})),
    ln_prob,
    -2.902929) # original conjugate, beta_mu = lnz
gce.reweight(-0.1)
plt.plot(gce.ln_prob().values())
minima = gce.ln_prob().minima()
assert(len(minima) == 1 and minima[0] == 152)
plt.gca().axvline(minima[0], color='black')
plt.text(20, -40, 'low density phase')
plt.text(200, -40, 'high density phase\n(more probable)')

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.

gce = pyfeasst.find_equilibrium(gce)

plt.plot(gce.ln_prob().values())
plt.gca().axvline(gce.ln_prob().minima()[0], color='black')
plt.text(20, -40, 'saturated vapor')
plt.text(225, -40, 'saturated liquid')

Text(225, -40, 'saturated liquid')

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(abs(-3.03 - gce.beta_mu()) < 0.01)
mins = gce.ln_prob().minima()
assert(len(mins) == 1 and mins[0] == 166)

# pressure
volume = 8**3
assert(abs(0.07722559602858005 - gce.betaPV(0)/beta/volume) < 1e-8)

# pressure of both phases are equivalent
assert(abs(gce.betaPV(0) - gce.betaPV(1)) < 1e-5)

# density of vapor
assert(abs(0.1003523021979188 - gce.average_macrostate(0)/volume) < 1e-8)

# density of liquid
assert(abs(0.5631867734067794 - gce.average_macrostate(1)/volume) < 1e-8)

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