1.4 Modifying the Data Space

It is very common that we discover at a later stage that we need to revise our computational protocol. In this case we need to carefully update existing job state points and the associated data.

Let’s assume that we realize that the ideal gas law is not sufficiently exact for our needs, so we’re going to use the van der Waals equation (vdW) for a more exact estimation of the volume for each state point

\(\left(p + \frac{N^2 a}{V^2} \right) (V - Nb) = N kT\), where

\(a\) and \(b\) are two additional parameters. For \(a=b=0\) this equation is identical to the ideal gas equation.

We start by implementing a function to calculate the volume for a given statepoint.

import numpy as np


def V_vdW(p, kT, N, a=0, b=0):
    """Solve the van der Waals equation for V."""
    coeffs = [p, -(kT * N + p * N * b), a * N**2, -a * N**3 * b]
    V = sorted(np.roots(coeffs))
    return np.real(V).tolist()

You will notice that this equation is a cubic polynomial and therefore has 3 possible solutions instead of only one!

print(V_vdW(1.0, 1.0, 1000))

[0.0, 0.0, 1000.0]

That is because the vdW system has a critical point and up to three possible solutions. These solutions correspond to a liquid, a gaseous and a meta-stable phase.

We want to make the old data compatible with the new protocol, which requires two modifications of the existing data space:

1. We need to add parameters \(a\) and \(b\) to each statepoint and set them to zero.

2. The former value V needs to be relabeled V_gas and we add a zero-value for V_liq.

We previously learned that we can use the Job.sp attribute interface to access individual state point values. We can use the same interface to modify existing state point parameters.

import signac

project = signac.get_project("projects/tutorial")

for job in project:
    if "a" not in job.sp:
        job.sp.a = 0
    if "b" not in job.sp:
        job.sp.b = 0

Please checkout the section on State Point Modifications in the reference documentation for a detailed description on how to modify state points.

Next, we need to update the existing volume data. We check whether the job document has a V value and replace it with V_liq and V_gas. The V.txt files will be rewritten to contain two comma-separated values.

for job in project:
    if "V" in job.document:
        job.document["V_liq"] = 0
        job.document["V_gas"] = job.document.pop("V")
        with open(job.fn("V.txt"), "w") as file:
            file.write("{},{}\n".format(0, job.document["V_gas"]))

Let’s verify our modifications!

for job in project:
    print(job.statepoint(), job.document)

{'p': 3.4000000000000004, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 294.1176470588235}
{'p': 7.800000000000001, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 128.2051282051282}
{'p': 1.0, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 1000.0}
{'p': 2.3000000000000003, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 434.78260869565213}
{'p': 8.9, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 112.35955056179775}
{'p': 5.05, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 198.01980198019803}
{'p': 2.575, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 388.34951456310677}
{'p': 4.5, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 222.22222222222223}
{'p': 7.525, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 132.89036544850498}
{'p': 0.1, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 10000.0}
{'p': 5.6, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 178.57142857142858}
{'p': 10.0, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 100.0}
{'p': 1.2000000000000002, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 833.3333333333333}
{'p': 6.7, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 149.2537313432836}

Next, we add a few state points with known parameters.

vdW = {
    # Source: https://en.wikipedia.org/wiki/Van_der_Waals_constants_(data_page)
    "ideal gas": {"a": 0, "b": 0},
    "argon": {"a": 1.355, "b": 0.03201},
    "water": {"a": 5.536, "b": 0.03049},
}


def calc_volume(job):
    V = V_vdW(**job.statepoint())
    job.document["V_liq"] = min(V)
    job.document["V_gas"] = max(V)
    with open(job.fn("V.txt"), "w") as file:
        file.write(f"{min(V)},{max(V)}\n")


for fluid in vdW:
    for p in np.linspace(0.1, 10.0, 10):
        sp = {"N": 1000, "p": float(p), "kT": 1.0}
        sp.update(vdW[fluid])
        job = project.open_job(sp)
        job.document["fluid"] = fluid
        calc_volume(job)

The fluid label is stored in the job document as a hint, which parameters were used, however they are deliberately not part of the state point, since our calculation is only based on the parameters N, kT, p, a, and b. In general, all state point variables should be independent of each other.

Let’s inspect the results:

ps = {job.statepoint()["p"] for job in project}
for fluid in sorted(vdW):
    print(fluid)
    for p in sorted(ps):
        jobs = project.find_jobs({"p": p, "doc.fluid": fluid})
        for job in jobs:
            print(
                round(p, 2),
                round(job.document["V_liq"], 4),
                round(job.document["V_gas"], 2),
            )
    print()

argon
0.1 32.8041 8430.94
1.2 32.8034 416.27
2.3 32.8027 216.99
3.4 32.8019 146.66
4.5 32.8012 110.72
5.6 32.8005 88.89
6.7 32.7998 74.23
7.8 32.799 63.71
8.9 32.7983 55.79
10.0 32.7976 49.61

ideal gas
0.1 0.0 10000.0
1.2 0.0 833.33
2.3 0.0 434.78
3.4 0.0 294.12
4.5 0.0 222.22
5.6 0.0 178.57
6.7 0.0 149.25
7.8 0.0 128.21
8.9 0.0 112.36
10.0 0.0 100.0

water
0.1 30.6598 4999.92
1.2 30.6598 416.58
2.3 30.6597 217.31
3.4 30.6597 146.97
4.5 30.6597 111.03
5.6 30.6596 89.2
6.7 30.6596 74.54
7.8 30.6596 64.02
8.9 30.6595 56.1
10.0 30.6595 49.92

We observe that the liquid phase is almost incompressible, while the gas phase is strongly pressure dependent.

The next section of the first chapter demonstrates how to use signac-flow for building a workflow.