Skip to content

Instantly share code, notes, and snippets.

@3ki5tj

3ki5tj/ljeos.html

Created November 18, 2014 21:11
Show Gist options
  • Save 3ki5tj/25d3cf6a116f126a61ba to your computer and use it in GitHub Desktop.
Save 3ki5tj/25d3cf6a116f126a61ba to your computer and use it in GitHub Desktop.
Download ZIP
Lennard-Jones equations of state (needs dygraph-combined.js for graphics)
Raw
ljeos.html
<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html;charset=utf-8">
<title>Lennard-Jones equation of state</title>
<script type="text/javascript" src="dygraph-combined.js"></script>
<script type="text/javascript">
//<![CDATA[
/* compute reference thermal dynamics variables
using the modified Benedic-Webb-Rubin (MBWR) equation of states
return
U: the average potential energy
P: pressure
Ar: Helmholtz free energy (potential part)
mu: Gibbs free energy (potential part) */
function lj_eos3dMBWR(rho, T, x)
{
var gam = x[33], U, P, Ar, mu;
var rhop, rho2 = rho*rho, Pa = 0., Pb = 0., i;
var a = [
x[1]*T + x[2]*Math.sqrt(T) + x[3] + x[4]/T + x[5]/(T*T),
x[6]*T + x[7] + x[8]/T + x[9]/(T*T),
x[10]*T + x[11] + x[12]/T,
x[13],
x[14]/T + x[15]/(T*T),
x[16]/T,
x[17]/T + x[18]/(T*T),
x[19]/(T*T)
];
var b = [
(x[20] + x[21]/T)/(T*T),
(x[22] + x[23]/(T*T))/(T*T),
(x[24] + x[25]/T)/(T*T),
(x[26] + x[27]/(T*T))/(T*T),
(x[28] + x[29]/T)/(T*T),
(x[30] + x[31]/T + x[32]/(T*T))/(T*T)
];
var c = [
x[2]*Math.sqrt(T)/2 + x[3] + 2*x[4]/T + 3*x[5]/(T*T),
x[7] + 2*x[8]/T + 3*x[9]/(T*T),
x[11] + 2*x[12]/T,
x[13],
2*x[14]/T + 3*x[15]/(T*T),
2*x[16]/T,
2*x[17]/T + 3*x[18]/(T*T),
3*x[19]/(T*T)
];
var d = [
(3*x[20] + 4*x[21]/T)/(T*T),
(3*x[22] + 5*x[23]/(T*T))/(T*T),
(3*x[24] + 4*x[25]/T)/(T*T),
(3*x[26] + 5*x[27]/(T*T))/(T*T),
(3*x[28] + 4*x[29]/T)/(T*T),
(3*x[30] + 4*x[31]/T + 5*x[32]/(T*T))/(T*T)
];
var F = Math.exp(-gam*rho*rho);
var G = [0, 0, 0, 0, 0, 0];
G[0] = (1 - F)/(2*gam);
for (rhop = 1, i = 1; i < 6; i++) {
rhop *= rho*rho;
G[i] = -(F*rhop - 2*i*G[i-1])/(2*gam);
}
Ar = 0;
Pa = Pb = 0;
for (U = 0, i = 7; i >= 0; i--) {
U = rho * (c[i]/(i+1) + U);
Ar = rho * (a[i]/(i+1) + Ar);
Pa = rho * (a[i] + Pa);
}
for (i = 5; i >= 0; i--) {
U += d[i]*G[i];
Ar += b[i]*G[i];
Pb = rho2*(b[i] + Pb);
}
P = rho*(T + Pa + F*Pb);
mu = Ar + P/rho - T;
return [U, P, Ar, mu];
}
/* Reference:
* J. Karl Johnson, John A. Zollweg, and Keith E. Gubbins
* The Lennard-Jones equation of states revisited,
* Molecular Physics (1993) Vol. 78, No 3, 591-618 */
var ljparam_JZG1993 = [0,
0.8623085097507421,
2.976218765822098,
-8.402230115796038,
0.1054136629203555,
-0.8564583828174598,
1.582759470107601,
0.7639421948305453,
1.753173414312048,
2.798291772190376e+03,
-4.8394220260857657e-02,
0.9963265197721935,
-3.698000291272493e+01,
2.084012299434647e+01,
8.305402124717285e+01,
-9.574799715203068e+02,
-1.477746229234994e+02,
6.398607852471505e+01,
1.603993673294834e+01,
6.805916615864377e+01,
-2.791293578795945e+03,
-6.245128304568454,
-8.116836104958410e+03,
1.488735559561229e+01,
-1.059346754655084e+04,
-1.131607632802822e+02,
-8.867771540418822e+03,
-3.986982844450543e+01,
-4.689270299917261e+03,
2.593535277438717e+02,
-2.694523589434903e+03,
-7.218487631550215e+02,
1.721802063863269e+02,
3.0 // gamma
];
/* Reference:
* Jiri Kolafa and Ivo Nezbeda
* The Lennard-Jones fluid: An accurate analytic
* and theoretically-based equation of state
* Fluid Phase Equilibria (1994) Vol. 100, 1-34
* TABLE 5
* regressed from data with T <= 6
* */
var ljparam_KN1994 = [0,
0.86230851,
2.97621877,
-8.40223012,
0.10541366,
-0.85645838,
1.39945300,
-0.20682219,
2.66555449,
1205.90355811,
0.24414200,
6.17927577,
-41.33848427,
15.14482295,
88.90243729,
-2425.74868591,
-148.52651854,
68.73779789,
2698.26346845,
-1216.87158315,
-1199.67930914,
-7.28265251,
-4942.58001124,
24.87520514,
-6246.96241113,
-235.12327760,
-7241.61133138,
-111.27706706,
-2800.52326352,
1109.71518240,
1455.47321956,
-2577.25311109,
476.67051504,
4.52000000 // gamma
];
function lj_eos3d_bAhs(eta)
{
var e1 = 1 - eta;
return Math.log(e1)*5/3 + eta*(34 - 33*eta + 4*eta*eta)/(6*e1*e1);
}
function lj_eos3d_zhs(eta)
{
var e1 = 1 - eta;
return (1 + eta*(1 + eta*(1 - eta*(1 + eta)*2/3)))/(e1*e1*e1);
}
/* equivalent hard-sphere diameter of the hybrid Barker-Henderson theory
* defined in Eq. (17)
* d = Int {0 to 2^(1/6)} (1 - exp(-beta*u)) dr
* Parameters are given by Table 2 with the functional form given by (29)
*
* To verify it WolframAlpha or Mathematica
* Integrate[(1-Exp[-(4/x^12-4/x^6+1)/T]), {x, 0, 2.0^(1.0/6)}]
* */
function lj_eos3d_dhBH(T)
{
var Ci = [1.080142248, -0.076383859, 0.011117524];
var C1 = 0.000693129, Cln = -0.063920968;
var x = Math.sqrt(T);
return [
Ci[2]/T + Ci[1]/x + Ci[0] + C1*x + Cln*Math.log(T),
Ci[2] + 0.5*Ci[1]*x - 0.5*C1*x*T - T*Cln]
}
/* The residual second virial coefficient B2_LJ - B2_hs
* from the hybrid Barker-Henderson theory
* Parameters are given by Table 2 with the functional form given by (29) */
function lj_eos3d_dB2hBH(T)
{
/* from Table 2 */
var Ci = [
0.02459877,
0,
-7.02181962,
2.90616279,
-4.13749995,
0.87361369,
0.43102052,
-0.58544978];
var f = 0, dfdbeta = 0, invx = 1/Math.sqrt(T), i;
for ( i = 7; i > 0; i-- ) {
f = (f + Ci[i]) * invx;
dfdbeta = (dfdbeta - Ci[i]*i/2) * invx;
}
dfdbeta *= -T;
f += Ci[0];
return [f, dfdbeta];
}
function lj_eos3dPVEhBH(rho, T)
{
/* Table 3 */
var Cij = [
[0, 0, /* i = 0 */
2.01546797,
-28.17881636,
28.28313847,
-10.42402873,
0],
[0, 0, /* i = -1 */
-19.58371655,
75.62340289,
-120.70586598,
93.92740328,
-27.37737354],
[0, 0, /* i = -2 */
29.34470520,
-112.35356937,
170.64908980,
-123.06669187,
34.42288969],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, /* i = -4 */
-13.37031968,
65.38059570,
-115.09233113,
88.91973082,
-25.62099890]];
var gam = 1.92907278;
var i, j, tmp;
tmp = lj_eos3d_dhBH(T);
var dia = tmp[0];
var ddia = tmp[1];
tmp = lj_eos3d_dB2hBH(T);
var dB2 = tmp[0];
var ddB2 = tmp[1];
var grho2 = gam*rho*rho;
var xpngrho2 = Math.exp(-grho2);
var eta = Math.PI * rho * dia*dia*dia/6;
var Ahs = T * lj_eos3d_bAhs(eta);
var AhBH = xpngrho2 * rho * T * dB2;
var ACij = 0, zCij = 0, UCij = 0;
for ( i = 0; i <= 4; i++ ) {
var xpT = Math.pow(T, -i*.5);
var xprho = rho;
for ( j = 2; j <= 6; j++ ) {
xprho *= rho;
var xp = xpT * xprho;
ACij += Cij[i][j] * xp;
zCij += j * Cij[i][j] * xp / T;
UCij += (1 + i*.5) * Cij[i][j] * xp;
}
}
var A = Ahs + AhBH + ACij;
var zhs = lj_eos3d_zhs(eta);
var z = zhs + rho * (1 - 2*grho2) * xpngrho2 * dB2 + zCij;
var P = rho * T * z;
var U = 3*(zhs - 1)*ddia/dia + rho * xpngrho2 * ddB2 + UCij;
var mu = P/rho + A - T;
return [U, P, A, mu];
}
function setnum(who, num)
{
var x = num.toPrecision(6);
if ( x < 0 ) x = "&minus;" + (-x);
document.getElementById(who).innerHTML = x;
}
function mkplot()
{
var rho, drho = 0.01;
var arr1, arr2, arr3;
var datU, datP, datF, datG, dat = "density,"
+"JZG1993,"
+"KN1994<sub>MBWR</sub>,"
+"KN1994<sub>PVE-hBH</sub>\n";
var temp = parseFloat( document.getElementById("temp").value );
if ( isNaN(temp) ) temp = 1;
var rhomax = parseFloat( document.getElementById("rho").value );
if ( isNaN(rhomax) ) {
rhomax = 1;
} else {
rhomax *= 1.2;
}
datU = datP = datF = datG = dat;
for ( var i = 1; i*drho <= rhomax; i++ ) {
rho = i * drho;
arr1 = lj_eos3dMBWR(rho, temp, ljparam_JZG1993);
arr2 = lj_eos3dMBWR(rho, temp, ljparam_KN1994);
arr3 = lj_eos3dPVEhBH(rho, temp);
datU += "" + rho + "," + arr1[0] + "," + arr2[0] + "," + arr3[0] + "\n";
datP += "" + rho + "," + arr1[1] + "," + arr2[1] + "," + arr3[1] + "\n";
datF += "" + rho + "," + arr1[2] + "," + arr2[2] + "," + arr3[2] + "\n";
datG += "" + rho + "," + arr1[3] + "," + arr2[3] + "," + arr3[3] + "\n";
}
var optionsU = {
showRoller: true,
ylabel: 'Internal energy per particle, <i>U</i>/<i>N</i>'
};
var optionsP = {
showRoller: true,
ylabel: 'Pressure, <i>P</i>'
};
var optionsF = {
showRoller: true,
xlabel: 'Density, <i>&rho;</i>',
ylabel: 'Helmholtz free energy per particle, <i>F</i><sup>ex</sup>/<i>N</i>'
};
var optionsG = {
showRoller: true,
xlabel: 'Density, <i>&rho;</i>',
ylabel: 'Excess chemical potential, <i>&mu;</i><sup>ex</sup>'
};
var gU = new Dygraph(document.getElementById("U_plot"), datU, optionsU);
var gP = new Dygraph(document.getElementById("P_plot"), datP, optionsP);
var gF = new Dygraph(document.getElementById("F_plot"), datF, optionsF);
var gG = new Dygraph(document.getElementById("G_plot"), datG, optionsG);
}
function getEOS()
{
var inpT = parseFloat( document.getElementById("temp").value );
var inprho = parseFloat( document.getElementById("rho").value );
if ( isNaN(inpT) ) inpT = 1;
if ( isNaN(inprho) ) inprho = 0.7;
var arr1 = lj_eos3dMBWR(inprho, inpT, ljparam_JZG1993);
setnum("U1", arr1[0]);
setnum("P1", arr1[1]);
setnum("Fex1", arr1[2]);
setnum("muex1", arr1[3]);
var arr2 = lj_eos3dMBWR(inprho, inpT, ljparam_KN1994);
setnum("U2", arr2[0]);
setnum("P2", arr2[1]);
setnum("Fex2", arr2[2]);
setnum("muex2", arr2[3]);
var arr3 = lj_eos3dPVEhBH(inprho, inpT);
setnum("U3", arr3[0]);
setnum("P3", arr3[1]);
setnum("Fex3", arr3[2]);
setnum("muex3", arr3[3]);
mkplot();
}
//]]>
</script>
</head>
<body>
<h2>Lennard-Jones equation of states</h2>
<form id="my_form">
<i>T</i>: <input type="text" id="temp" value="1" size="10" />
<i>&rho;</i>: <input type="text" id="rho" value="0.7" size="10" />
<input type="button" value="Compute" size="10" onClick="getEOS()" />
<br />
<table style="text-align: center; margin: 2px; border: 2px solid">
<tr>
<th></th>
<th style="text-align:left">Johnson, Zollweg, and Gubbins (1993)</th>
<th style="text-align:left">Kolafa and Nezbeda, MBWR (1994)</th>
<th style="text-align:left">Kolafa and Nezbeda, PVE-hBH (1994)</th>
</tr>
<tr>
<td style="text-align:left">
Internal energy per particle, <i>U/N</i>:
</td>
<td id="U1"></td>
<td id="U2"></td>
<td id="U3"></td>
</tr>
<tr>
<td style="text-align:left">
Pressure, <i>P</i>:
</td>
<td id="P1"></td>
<td id="P2"></td>
<td id="P3"></td>
</tr>
<tr>
<td style="text-align:left">
Excess free energy per particle, <i>F/N</i>:
</td>
<td id="Fex1"></td>
<td id="Fex2"></td>
<td id="Fex3"></td>
</tr>
<tr>
<td style="text-align:left">
Excess chemical potential, <i>&mu;</i>:
</td>
<td id="muex1"></td>
<td id="muex2"></td>
<td id="muex3"></td>
</tr>
</table>
<center>
<input type="button" value="Plot!"
onClick="mkplot()"
style="padding: 5px 20px; width: 400px; margin:auto; text-align:center;">
</center>
<table>
<tr>
<td> <div id="U_plot"></div> </td>
<td> <div id="P_plot"></div> </td>
</tr>
<tr>
<td> <div id="F_plot"></div> </td>
<td> <div id="G_plot"></div> </td>
</tr>
</table>
</form>
</body>
</html>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment