References: “Introduction to Statistical Thermodynamics” (c) 1960 by Terrell L. Hill 8.14.2: Properties of Nobel Gases – Chemistry LibreTexts
First, we approximate the hard sphere second virial coefficient for the Noble Gases using Monte Carlo integration. Analytic integration yields:
B2(T) = 2 * pi * a ^ 3 / 3
Where a is the diameter of the atom. The volume utilized was V = X ^ 3 where X = 25.0 * 100.0 * a. Then we compute four cluster integrals.
// HardSphereZ3.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Sunday, November 1, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// all potential energy is for a hard sphere
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
// second virial coefficient
static double B2(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
if (r12 >= a)
I += 1.0;
}
return V * I / N;
}
// three particle cluster integral
static double Z30(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r13 = std::sqrt(
std::pow(x1 - x3, 2.0) +
std::pow(y1 - y3, 2.0) +
std::pow(z1 - z3, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r13 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
// four particle cluster integral
static double Z40(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double x4 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double y4 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double z4 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double r34 = std::sqrt(
std::pow(x3 - x4, 2.0) +
std::pow(y3 - y4, 2.0) +
std::pow(z3 - z4, 2.0));
double r41 = std::sqrt(
std::pow(x4 - x1, 2.0) +
std::pow(y4 - y1, 2.0) +
std::pow(z4 - z1, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (r34 >= a)
u += 1.0;
if (r41 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
static double Z41(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double x4 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double y4 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double z4 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double r34 = std::sqrt(
std::pow(x3 - x4, 2.0) +
std::pow(y3 - y4, 2.0) +
std::pow(z3 - z4, 2.0));
double r41 = std::sqrt(
std::pow(x4 - x1, 2.0) +
std::pow(y4 - y1, 2.0) +
std::pow(z4 - z1, 2.0));
double r24 = std::sqrt(
std::pow(x2 - x4, 2.0) +
std::pow(y2 - y4, 2.0) +
std::pow(z2 - z4, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (r34 >= a)
u += 1.0;
if (r41 >= a)
u += 1.0;
if (r24 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
static double Z42(double X, double a, int N)
{
double I = 0.0, V = X * X * X;
for (int n = 0; n < N; n++)
{
double x1 = (X * rand()) / RAND_MAX;
double x2 = (X * rand()) / RAND_MAX;
double x3 = (X * rand()) / RAND_MAX;
double x4 = (X * rand()) / RAND_MAX;
double y1 = (X * rand()) / RAND_MAX;
double y2 = (X * rand()) / RAND_MAX;
double y3 = (X * rand()) / RAND_MAX;
double y4 = (X * rand()) / RAND_MAX;
double z1 = (X * rand()) / RAND_MAX;
double z2 = (X * rand()) / RAND_MAX;
double z3 = (X * rand()) / RAND_MAX;
double z4 = (X * rand()) / RAND_MAX;
double r12 = std::sqrt(
std::pow(x1 - x2, 2.0) +
std::pow(y1 - y2, 2.0) +
std::pow(z1 - z2, 2.0));
double r23 = std::sqrt(
std::pow(x2 - x3, 2.0) +
std::pow(y2 - y3, 2.0) +
std::pow(z2 - z3, 2.0));
double r34 = std::sqrt(
std::pow(x3 - x4, 2.0) +
std::pow(y3 - y4, 2.0) +
std::pow(z3 - z4, 2.0));
double r41 = std::sqrt(
std::pow(x4 - x1, 2.0) +
std::pow(y4 - y1, 2.0) +
std::pow(z4 - z1, 2.0));
double r24 = std::sqrt(
std::pow(x2 - x4, 2.0) +
std::pow(y2 - y4, 2.0) +
std::pow(z2 - z4, 2.0));
double r31 = std::sqrt(
std::pow(x3 - x1, 2.0) +
std::pow(y3 - y1, 2.0) +
std::pow(z3 - z1, 2.0));
double u = 0.0;
if (r12 >= a)
u += 1.0;
if (r23 >= a)
u += 1.0;
if (r34 >= a)
u += 1.0;
if (r41 >= a)
u += 1.0;
if (r24 >= a)
u += 1.0;
if (r31 >= a)
u += 1.0;
if (u > 0.0)
I += u;
}
return V * I / N;
}
int main()
{
double pi = 4.0 * atan(1.0);
// radii of Noble Gases in picometers (10^-12 m)
const double Radii[6] = { 31, 38, 71, 88, 108, 120 };
// spherically symmetric atom's name
const std::string Names[6] = {
"Helium", "Neon", "Argon", "Krypton", "Xenon", "Radon" };
std::cout << "Hard Sphere Second Virial Coefficients for Noble Gases";
std::cout << std::endl;
std::cout << "Monte Carlo Integration is Used";
std::cout << std::endl;
std::cout << "B2(T)" << std::endl;
std::cout << "Atom\tExperiment\tTheory" << std::endl;
srand(1);
for (int i = 0; i < 6; i++)
{
double a = 2.0 * Radii[i] * 1.0e-12;
double X = 25.0 * 100.0 * a;
double b2 = B2(X, a, 1000000);
double bt = 2.0 * pi * std::pow(a, 3.0) / 3.0;
std::cout << Names[i] << "\t";
std::cout << std::scientific;
std::cout << std::setprecision(6);
std::cout << std::right;
std::cout << b2 << '\t';
std::cout << bt << std::endl;
}
std::cout << "Hard Sphere Configuration Integrals";
std::cout << std::endl;
std::cout << "Atom\tZ30\t\tZ40\t\tZ41\t\tZ42";
std::cout << std::endl;
for (int i = 0; i < 6; i++)
{
double a = 2.0 * Radii[i] * 1.0e-12;
double X = 25.0 * 100.0 * a;
double z30 = Z30(X, a, 1000000);
double z40 = Z40(X, a, 1000000);
double z41 = Z41(X, a, 1000000);
double z42 = Z42(X, a, 1000000);
std::cout << Names[i] << '\t';
std::cout << std::scientific;
std::cout << std::setprecision(6);
std::cout << std::right;
std::cout << z30 << '\t';
std::cout << z40 << '\t';
std::cout << z41 << '\t';
std::cout << z42 << '\t';
std::cout << std::endl;
}
return 0;
}
