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.
In [1]:
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!
In [2]:
print(V_vdW(1.0, 1.0, 1000))
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:
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.
In [3]:
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.
In [4]:
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!
In [5]:
for job in project:
print(job.statepoint(), job.document)
Next, we add a few state points with known parameters.
In [6]:
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('{},{}\n'.format(min(V), max(V)))
for fluid in vdW:
for p in np.linspace(0.1, 10.0, 10):
sp = {'N': 1000, 'p': 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:
In [7]:
ps = set((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_filter={'fluid': fluid})
for job in jobs:
print(round(p, 2), round(job.document['V_liq'], 4), round(job.document['V_gas'], 2))
print()
We observe that the liquid phase is almost incompressible, while the gas phase is strongly pressure dependent.
The final section of the first chapter demonstrates how to interact with a signac on the command line.