notebooks/signac_104_Modifying_the_Data_Space.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.4 Modifying the Data Space\n",
    "\n",
    "It is very common that we discover at a later stage that we need to revise our computational protocol.\n",
    "In this case we need to carefully update existing job state points and the associated data.\n",
    "\n",
    "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)](https://en.wikipedia.org/wiki/Van_der_Waals_equation) for a more exact estimation of the volume for each state point\n",
    "\n",
    "$\\left(p + \\frac{N^2 a}{V^2} \\right) (V - Nb) = N kT$, where\n",
    "\n",
    "$a$ and $b$ are two additional parameters.\n",
    "For $a=b=0$ this equation is identical to the ideal gas equation.\n",
    "\n",
    "We start by implementing a function to calculate the volume for a given statepoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def V_vdW(p, kT, N, a=0, b=0):\n",
    "    \"\"\"Solve the van der Waals equation for V.\"\"\"\n",
    "    coeffs = [p, -(kT * N + p * N * b), a * N**2, -a * N**3 * b]\n",
    "    V = sorted(np.roots(coeffs))\n",
    "    return np.real(V).tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will notice that this equation is a cubic polynomial and therefore has [3 possible solutions](http://docs.scipy.org/doc/numpy/reference/generated/numpy.roots.html) instead of only one!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 0.0, 1000.0]\n"
     ]
    }
   ],
   "source": [
    "print(V_vdW(1.0, 1.0, 1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is because the vdW system has a critical point and up to three possible solutions.\n",
    "These solutions correspond to a liquid, a gaseous and a meta-stable phase.\n",
    "\n",
    "We want to make the old data compatible with the new protocol, which requires two modifications of the existing data space:\n",
    "\n",
    "1. We need to add parameters $a$ and $b$ to each statepoint and set them to zero.\n",
    "2. The former value `V` needs to be relabeled `V_gas` and we add a zero-value for `V_liq`.\n",
    "\n",
    "We previously learned that we can use the `Job.sp` attribute interface to access individual state point values.\n",
    "We can use the same interface to modify existing state point parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import signac\n",
    "\n",
    "project = signac.get_project(\"projects/tutorial\")\n",
    "\n",
    "for job in project:\n",
    "    if \"a\" not in job.sp:\n",
    "        job.sp.a = 0\n",
    "    if \"b\" not in job.sp:\n",
    "        job.sp.b = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please checkout the section on [State Point Modifications](https://signac.readthedocs.io/en/latest/projects.html#state-point-modifications) in the reference documentation for a detailed description on how to modify state points.\n",
    "\n",
    "Next, we need to update the existing volume data.\n",
    "We check whether the *job document* has a `V` value and replace it with `V_liq` and `V_gas`.\n",
    "The `V.txt` files will be rewritten to contain two comma-separated values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "for job in project:\n",
    "    if \"V\" in job.document:\n",
    "        job.document[\"V_liq\"] = 0\n",
    "        job.document[\"V_gas\"] = job.document.pop(\"V\")\n",
    "        with open(job.fn(\"V.txt\"), \"w\") as file:\n",
    "            file.write(\"{},{}\\n\".format(0, job.document[\"V_gas\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's verify our modifications!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'p': 3.4000000000000004, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 294.1176470588235}\n",
      "{'p': 7.800000000000001, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 128.2051282051282}\n",
      "{'p': 1.0, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 1000.0}\n",
      "{'p': 2.3000000000000003, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 434.78260869565213}\n",
      "{'p': 8.9, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 112.35955056179775}\n",
      "{'p': 5.05, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 198.01980198019803}\n",
      "{'p': 2.575, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 388.34951456310677}\n",
      "{'p': 4.5, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 222.22222222222223}\n",
      "{'p': 7.525, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 132.89036544850498}\n",
      "{'p': 0.1, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 10000.0}\n",
      "{'p': 5.6, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 178.57142857142858}\n",
      "{'p': 10.0, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 100.0}\n",
      "{'p': 1.2000000000000002, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 833.3333333333333}\n",
      "{'p': 6.7, 'kT': 1.0, 'N': 1000, 'a': 0, 'b': 0} {'V_liq': 0, 'V_gas': 149.2537313432836}\n"
     ]
    }
   ],
   "source": [
    "for job in project:\n",
    "    print(job.statepoint(), job.document)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we add a few state points with known parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "vdW = {\n",
    "    # Source: https://en.wikipedia.org/wiki/Van_der_Waals_constants_(data_page)\n",
    "    \"ideal gas\": {\"a\": 0, \"b\": 0},\n",
    "    \"argon\": {\"a\": 1.355, \"b\": 0.03201},\n",
    "    \"water\": {\"a\": 5.536, \"b\": 0.03049},\n",
    "}\n",
    "\n",
    "\n",
    "def calc_volume(job):\n",
    "    V = V_vdW(**job.statepoint())\n",
    "    job.document[\"V_liq\"] = min(V)\n",
    "    job.document[\"V_gas\"] = max(V)\n",
    "    with open(job.fn(\"V.txt\"), \"w\") as file:\n",
    "        file.write(f\"{min(V)},{max(V)}\\n\")\n",
    "\n",
    "\n",
    "for fluid in vdW:\n",
    "    for p in np.linspace(0.1, 10.0, 10):\n",
    "        sp = {\"N\": 1000, \"p\": float(p), \"kT\": 1.0}\n",
    "        sp.update(vdW[fluid])\n",
    "        job = project.open_job(sp)\n",
    "        job.document[\"fluid\"] = fluid\n",
    "        calc_volume(job)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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*.\n",
    "In general, all state point variables should be independent of each other.\n",
    "\n",
    "Let's inspect the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "argon\n",
      "0.1 32.8041 8430.94\n",
      "1.2 32.8034 416.27\n",
      "2.3 32.8027 216.99\n",
      "3.4 32.8019 146.66\n",
      "4.5 32.8012 110.72\n",
      "5.6 32.8005 88.89\n",
      "6.7 32.7998 74.23\n",
      "7.8 32.799 63.71\n",
      "8.9 32.7983 55.79\n",
      "10.0 32.7976 49.61\n",
      "\n",
      "ideal gas\n",
      "0.1 0.0 10000.0\n",
      "1.2 0.0 833.33\n",
      "2.3 0.0 434.78\n",
      "3.4 0.0 294.12\n",
      "4.5 0.0 222.22\n",
      "5.6 0.0 178.57\n",
      "6.7 0.0 149.25\n",
      "7.8 0.0 128.21\n",
      "8.9 0.0 112.36\n",
      "10.0 0.0 100.0\n",
      "\n",
      "water\n",
      "0.1 30.6598 4999.92\n",
      "1.2 30.6598 416.58\n",
      "2.3 30.6597 217.31\n",
      "3.4 30.6597 146.97\n",
      "4.5 30.6597 111.03\n",
      "5.6 30.6596 89.2\n",
      "6.7 30.6596 74.54\n",
      "7.8 30.6596 64.02\n",
      "8.9 30.6595 56.1\n",
      "10.0 30.6595 49.92\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ps = {job.statepoint()[\"p\"] for job in project}\n",
    "for fluid in sorted(vdW):\n",
    "    print(fluid)\n",
    "    for p in sorted(ps):\n",
    "        jobs = project.find_jobs({\"p\": p, \"doc.fluid\": fluid})\n",
    "        for job in jobs:\n",
    "            print(\n",
    "                round(p, 2),\n",
    "                round(job.document[\"V_liq\"], 4),\n",
    "                round(job.document[\"V_gas\"], 2),\n",
    "            )\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the liquid phase is almost incompressible, while the gas phase is strongly pressure dependent.\n",
    "\n",
    "The [next section](signac_105_signac-flow_Ideal_Gas_Example.ipynb) of the first chapter demonstrates how to use **signac-flow** for building a workflow."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}