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

The python package pyfeasst (https://pypi.org/project/pyfeasst/) is used for analysis.

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

[1]:
import numpy as np
import pandas as pd
from pyfeasst import macrostate_distribution

srsw = macrostate_distribution.read_csv(file_name="../test/data/stat120.csv",
                                        macrostate_name="N",
                                        ln_prob_name="lnPI")
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
plt.plot(srsw.ln_prob())
[1]:
[<matplotlib.lines.Line2D at 0x7f880265b700>]
../../../_images/plugin_flat_histogram_tutorial_tutorial_03_reweight_phase_separation_1_1.svg

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.

[2]:
assert len(srsw.minimums())==0

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

[3]:
srsw_rw = srsw.reweight(delta_beta_mu=-0.1)
plt.plot(srsw_rw.ln_prob())
minima=srsw_rw.minimums()
assert len(minima) == 1 and minima.values[0] == 152
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$')
[3]:
Text(0, 0.5, '$\\ln\\Pi$')
../../../_images/plugin_flat_histogram_tutorial_tutorial_03_reweight_phase_separation_5_1.svg

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

[4]:
delta_beta_mu_equilibrium = srsw.equilibrium()
mins=srsw.minimums()
assert len(mins)==1
mn=mins.values[0]

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$')
[4]:
Text(0, 0.5, '$\\ln\\Pi$')
../../../_images/plugin_flat_histogram_tutorial_tutorial_03_reweight_phase_separation_7_1.svg

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

[5]:
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.563188088494581 - liquid_density) < 1e-8
pressure: 0.07722557238264459
vapor_density: 0.10035100340092687
liquid_density: 0.5631880884945807

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