{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7a703c5f",
   "metadata": {},
   "source": [
    "# Notebook05: POPC/cholesterol mixture\n",
    "\n",
    "In this notebook, we will show you how to reconstruct hydrogens, calculate the order parameters and produce output trajectories on a POPC/cholesterol (66:34) mixture. Again, this example is based on the Berger united-atom force field.\n",
    "\n",
    "Before going on, we advise you to get started with [Notebook01](Notebook_01_buildH_calc_OP.ipynb) if you are not familiar with basic features of **buildH** and Jupyter notebooks (e.g. do you know what `!` means in a Jupyter cell?). \n",
    "\n",
    "## Checking buildH activation\n",
    "\n",
    "As explained in [Notebook01](Notebook_01_buildH_calc_OP.ipynb), you should have activated **buildH** before launching this notebook, thus you should obtain the following when invoking `buildH` in a Unix terminal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9fe5a816",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: buildH [-h] [-v] -c COORD [-t TRAJ] -l LIPID\n",
      "              [-lt LIPID_TOPOLOGY [LIPID_TOPOLOGY ...]] -d DEFOP\n",
      "              [-opx OPDBXTC] [-o OUT] [-b BEGIN] [-e END] [-igch3]\n",
      "buildH: error: the following arguments are required: -c/--coord, -l/--lipid, -d/--defop\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "!buildH"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "808ca76e",
   "metadata": {
    "tags": []
   },
   "source": [
    "If you see this, you can go on. If not, please go to `Notebook01` and see the explanations there on how to activate **buildH**.\n",
    "\n",
    "## Downloading the example files\n",
    "\n",
    "The example shown in this Notebook are taken from a POPC/cholesterol (66:34) trajectory (84 POPC molecules and 44 cholesterol molecules) taken from the [NMRlipidsIVb project](http://nmrlipids.blogspot.com/). The initial source files can be found on [Zenodo](https://zenodo.org/record/4643985). For the purpose of learning, we simplified the trajectory and extracted only 25 frames (0-25 ns, one frame / ns). We also made the molecules whole with the tool [trjconv](https://manual.gromacs.org/documentation/current/onlinehelp/gmx-trjconv.html) from GROMACS with the flag `-pc mol`. The starting gro file and the trajectory of 25 frames can be directly downloaded from the **buildH** repo in the directory https://github.com/patrickfuchs/buildH/tree/master/docs/Berger_POPCCHOL_test_case. Here we'll download them with `wget` (recall to download all files in raw format):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0dd9c85a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2021-07-07 15:17:17--  https://raw.githubusercontent.com/patrickfuchs/buildH/master/docs/Berger_POPCCHOL_test_case/endCONF_OK.gro\n",
      "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
      "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 1795883 (1.7M) [text/plain]\n",
      "Saving to: ‘endCONF_OK.gro’\n",
      "\n",
      "endCONF_OK.gro      100%[===================>]   1.71M  --.-KB/s    in 0.08s   \n",
      "\n",
      "2021-07-07 15:17:17 (22.4 MB/s) - ‘endCONF_OK.gro’ saved [1795883/1795883]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://raw.githubusercontent.com/patrickfuchs/buildH/master/docs/Berger_POPCCHOL_test_case/endCONF_OK.gro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "28b93568",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2021-07-07 15:17:20--  https://github.com/patrickfuchs/buildH/raw/master/docs/Berger_POPCCHOL_test_case/popcCHOL34molPER0-25ns_OK_dt1000.xtc\n",
      "Resolving github.com (github.com)... 140.82.121.3\n",
      "Connecting to github.com (github.com)|140.82.121.3|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 Found\n",
      "Location: https://raw.githubusercontent.com/patrickfuchs/buildH/master/docs/Berger_POPCCHOL_test_case/popcCHOL34molPER0-25ns_OK_dt1000.xtc [following]\n",
      "--2021-07-07 15:17:20--  https://raw.githubusercontent.com/patrickfuchs/buildH/master/docs/Berger_POPCCHOL_test_case/popcCHOL34molPER0-25ns_OK_dt1000.xtc\n",
      "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.111.133, 185.199.108.133, 185.199.109.133, ...\n",
      "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.111.133|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 2472992 (2.4M) [application/octet-stream]\n",
      "Saving to: ‘popcCHOL34molPER0-25ns_OK_dt1000.xtc’\n",
      "\n",
      "popcCHOL34molPER0-2 100%[===================>]   2.36M  --.-KB/s    in 0.09s   \n",
      "\n",
      "2021-07-07 15:17:21 (25.8 MB/s) - ‘popcCHOL34molPER0-25ns_OK_dt1000.xtc’ saved [2472992/2472992]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://github.com/patrickfuchs/buildH/raw/master/docs/Berger_POPCCHOL_test_case/popcCHOL34molPER0-25ns_OK_dt1000.xtc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "87e085c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data\n",
      "endCONF_OK.gro\n",
      "img\n",
      "Notebook_01_buildH_calc_OP.ipynb\n",
      "Notebook_02_buildH_calc_OP_outputwH.ipynb\n",
      "Notebook_03_mixture_POPC_POPE.ipynb\n",
      "Notebook_04_library.ipynb\n",
      "Notebook_05_mixture_POPC_cholesterol.ipynb\n",
      "popcCHOL34molPER0-25ns_OK_dt1000.xtc\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3919d24",
   "metadata": {},
   "source": [
    "We also need the two def files corresponding to POPC and cholesterol. They can be downloaded from the **buildH** repo in the directory https://github.com/patrickfuchs/buildH/tree/master/def_files (again remember to use raw files):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d4f09307",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2021-07-07 15:17:27--  https://raw.githubusercontent.com/patrickfuchs/buildH/master/def_files/Berger_POPC.def\n",
      "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.110.133, 185.199.111.133, 185.199.108.133, ...\n",
      "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.110.133|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 2135 (2.1K) [text/plain]\n",
      "Saving to: ‘Berger_POPC.def’\n",
      "\n",
      "Berger_POPC.def     100%[===================>]   2.08K  --.-KB/s    in 0s      \n",
      "\n",
      "2021-07-07 15:17:27 (8.49 MB/s) - ‘Berger_POPC.def’ saved [2135/2135]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://raw.githubusercontent.com/patrickfuchs/buildH/master/def_files/Berger_POPC.def"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f8fe0ee4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2021-07-07 15:17:30--  https://raw.githubusercontent.com/patrickfuchs/buildH/master/def_files/Berger_CHOL.def\n",
      "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
      "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 1275 (1.2K) [text/plain]\n",
      "Saving to: ‘Berger_CHOL.def’\n",
      "\n",
      "Berger_CHOL.def     100%[===================>]   1.25K  --.-KB/s    in 0s      \n",
      "\n",
      "2021-07-07 15:17:31 (5.90 MB/s) - ‘Berger_CHOL.def’ saved [1275/1275]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://raw.githubusercontent.com/patrickfuchs/buildH/master/def_files/Berger_CHOL.def"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cb8a936b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Berger_CHOL.def\t\t\t  Notebook_02_buildH_calc_OP_outputwH.ipynb\n",
      "Berger_POPC.def\t\t\t  Notebook_03_mixture_POPC_POPE.ipynb\n",
      "data\t\t\t\t  Notebook_04_library.ipynb\n",
      "endCONF_OK.gro\t\t\t  Notebook_05_mixture_POPC_cholesterol.ipynb\n",
      "img\t\t\t\t  popcCHOL34molPER0-25ns_OK_dt1000.xtc\n",
      "Notebook_01_buildH_calc_OP.ipynb\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91b6ec3d",
   "metadata": {},
   "source": [
    "Recall, the def files present on the **buildH** `def_files` directory contain all the possible (apolar) H that can be reconstructed on each lipids. The polar hydrogen of cholesterol is already present, thus it does not need to be reconstructed. **buildH** will let it as is, like all other non apolar H atoms.\n",
    "\n",
    "## Dealing with POPC\n",
    "\n",
    "When we have a mixture of lipids, **buildH** cannot reconstruct hydrogens on all lipids at the same time. So we start with H reconstruction and OP calculation on POPC (we'll focus on cholesterol later). Since Berger POPC is a supported lipid, we can launch easily **buildH** with the classical flags: `-c` for reading the initial pdb or gro file, `-t` for reading the trajectory, `-l` for specifying the type of lipid, `-d` for specifying the def file, `-o` for the output file with order parameters, and last `-opx` for the output trajectory with hydrogens (recall not to put any extension, they'll be added automatically by **buildH**). Importantly, when we want to output a trajectory with newly reconstructed hydrogens (option `-opx`), *all possible H to reconstruct need to be present in the def file*. This is fine since this is the case for all lipids in the [def file directory](https://github.com/patrickfuchs/buildH/tree/master/def_files) of the  **buildH** repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "79174caf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Constructing the system...\n",
      "System has 26026 atoms\n",
      "Writing new pdb with hydrogens.\n",
      "Writing trajectory with hydrogens in xtc file.\n",
      "Dealing with frame 0 at 30000.0 ps.\n",
      "Dealing with frame 1 at 31000.0 ps.\n",
      "Dealing with frame 2 at 32000.0 ps.\n",
      "Dealing with frame 3 at 33000.0 ps.\n",
      "Dealing with frame 4 at 34000.0 ps.\n",
      "Dealing with frame 5 at 35000.0 ps.\n",
      "Dealing with frame 6 at 36000.0 ps.\n",
      "Dealing with frame 7 at 37000.0 ps.\n",
      "Dealing with frame 8 at 38000.0 ps.\n",
      "Dealing with frame 9 at 39000.0 ps.\n",
      "Dealing with frame 10 at 40000.0 ps.\n",
      "Dealing with frame 11 at 41000.0 ps.\n",
      "Dealing with frame 12 at 42000.0 ps.\n",
      "Dealing with frame 13 at 43000.0 ps.\n",
      "Dealing with frame 14 at 44000.0 ps.\n",
      "Dealing with frame 15 at 45000.0 ps.\n",
      "Dealing with frame 16 at 46000.0 ps.\n",
      "Dealing with frame 17 at 47000.0 ps.\n",
      "Dealing with frame 18 at 48000.0 ps.\n",
      "Dealing with frame 19 at 49000.0 ps.\n",
      "Dealing with frame 20 at 50000.0 ps.\n",
      "Dealing with frame 21 at 51000.0 ps.\n",
      "Dealing with frame 22 at 52000.0 ps.\n",
      "Dealing with frame 23 at 53000.0 ps.\n",
      "Dealing with frame 24 at 54000.0 ps.\n",
      "Dealing with frame 25 at 55000.0 ps.\n",
      "Results written to OP_POPC.out\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "!buildH -c endCONF_OK.gro -t popcCHOL34molPER0-25ns_OK_dt1000.xtc -l Berger_PLA -d Berger_POPC.def -o OP_POPC.out -opx traj_POPC_wH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c459ebcd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Berger_CHOL.def\n",
      "Berger_POPC.def\n",
      "data\n",
      "endCONF_OK.gro\n",
      "img\n",
      "Notebook_01_buildH_calc_OP.ipynb\n",
      "Notebook_02_buildH_calc_OP_outputwH.ipynb\n",
      "Notebook_03_mixture_POPC_POPE.ipynb\n",
      "Notebook_04_library.ipynb\n",
      "Notebook_05_mixture_POPC_cholesterol.ipynb\n",
      "OP_POPC.out\n",
      "popcCHOL34molPER0-25ns_OK_dt1000.xtc\n",
      "traj_POPC_wH.pdb\n",
      "traj_POPC_wH.xtc\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f11f0273",
   "metadata": {},
   "source": [
    "The order parameters of POPC lipids have been written to the file `OP_POPC.dat`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "69396299",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# OP_name            resname atom1 atom2  OP_mean OP_stddev OP_stem\n",
      "#--------------------------------------------------------------------\n",
      "gamma1_1             PLA     C1    H11    0.03609  0.16545  0.01805\n",
      "gamma1_2             PLA     C1    H12   -0.01327  0.09031  0.00985\n",
      "gamma1_3             PLA     C1    H13   -0.01714  0.11162  0.01218\n",
      "gamma2_1             PLA     C2    H21    0.03573  0.16452  0.01795\n",
      "gamma2_2             PLA     C2    H22   -0.01077  0.10390  0.01134\n",
      "gamma2_3             PLA     C2    H23    0.00015  0.09757  0.01065\n",
      "gamma3_1             PLA     C3    H31    0.03542  0.16105  0.01757\n",
      "gamma3_2             PLA     C3    H32   -0.01826  0.09667  0.01055\n"
     ]
    }
   ],
   "source": [
    "!head OP_POPC.out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9194ef3",
   "metadata": {},
   "source": [
    "A pdb `traj_POPC_wH.pdb` and an xtc `traj_POPC_wH.xtc` file have been produced with hydrogens on POPC only. So let's check whether the POPC have hydrogens and not cholesterol (recall, here POPC residue name is PLA):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f70f9cc9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATOM      1  C1   PLA    1      10.670  42.050  27.850  1.00  0.00             C\n",
      "ATOM      2  H11  PLA    1      10.727  41.123  27.280  1.00  0.00             H\n",
      "ATOM      3  H12  PLA    1       9.962  41.929  28.670  1.00  0.00             H\n",
      "ATOM      4  H13  PLA    1      11.654  42.290  28.252  1.00  0.00             H\n",
      "ATOM      5  C2   PLA    1       8.900  42.780  26.440  1.00  0.00             C\n",
      "ATOM      6  H21  PLA    1       8.966  41.826  25.916  1.00  0.00             H\n",
      "ATOM      7  H22  PLA    1       8.564  43.552  25.748  1.00  0.00             H\n",
      "ATOM      8  H23  PLA    1       8.189  42.694  27.262  1.00  0.00             H\n",
      "ATOM      9  C3   PLA    1      11.120  43.130  25.800  1.00  0.00             C\n",
      "ATOM     10  H31  PLA    1      11.123  42.136  25.353  1.00  0.00             H\n",
      "ATOM     11  H32  PLA    1      12.130  43.389  26.117  1.00  0.00             H\n",
      "ATOM     12  H33  PLA    1      10.771  43.857  25.067  1.00  0.00             H\n",
      "ATOM     13  N4   PLA    1      10.220  43.140  26.970  1.00  0.00             N\n",
      "ATOM     14  C5   PLA    1      10.120  44.420  27.670  1.00  0.00             C\n",
      "ATOM     15  H51  PLA    1       9.433  44.285  28.506  1.00  0.00             H\n",
      "ATOM     16  H52  PLA    1       9.704  45.146  26.971  1.00  0.00             H\n",
      "ATOM     17  C6   PLA    1      11.420  45.010  28.230  1.00  0.00             C\n",
      "ATOM     18  H61  PLA    1      12.062  44.193  28.559  1.00  0.00             H\n",
      "ATOM     19  H62  PLA    1      11.177  45.647  29.081  1.00  0.00             H\n",
      "ATOM     20  O7   PLA    1      12.130  45.790  27.260  1.00  0.00             O\n",
      "grep: write error: Broken pipe\n"
     ]
    }
   ],
   "source": [
    "!grep PLA traj_POPC_wH.pdb | head -20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "049bbfae",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATOM   2681  C1  CHOL   21       8.750  35.790  36.370  1.00  0.00             C\n",
      "ATOM   2682  C2  CHOL   21       9.800  36.660  35.670  1.00  0.00             C\n",
      "ATOM   2683  C3  CHOL   21       9.370  37.990  35.030  1.00  0.00             C\n",
      "ATOM   2684  C4  CHOL   21       8.600  37.850  33.710  1.00  0.00             C\n",
      "ATOM   2685  C5  CHOL   21       9.340  36.970  32.700  1.00  0.00             C\n",
      "ATOM   2686  O6  CHOL   21       8.500  36.770  31.560  1.00  0.00             O\n",
      "ATOM   2687  H   CHOL   21       8.170  37.670  31.260  1.00  0.00             H\n",
      "ATOM   2688  C8  CHOL   21       9.890  35.630  33.200  1.00  0.00             C\n",
      "ATOM   2689  C9  CHOL   21      10.480  35.810  34.600  1.00  0.00             C\n",
      "ATOM   2690  C10 CHOL   21      11.690  35.180  34.860  1.00  0.00             C\n",
      "ATOM   2691  C11 CHOL   21      12.450  35.230  36.190  1.00  0.00             C\n",
      "ATOM   2692  C12 CHOL   21      11.660  35.960  37.280  1.00  0.00             C\n",
      "ATOM   2693  C13 CHOL   21      10.840  37.110  36.690  1.00  0.00             C\n",
      "ATOM   2694  C14 CHOL   21      10.220  37.900  37.850  1.00  0.00             C\n",
      "ATOM   2695  C15 CHOL   21      11.100  38.390  39.000  1.00  0.00             C\n",
      "ATOM   2696  C16 CHOL   21      11.860  37.210  39.620  1.00  0.00             C\n",
      "ATOM   2697  C17 CHOL   21      10.960  36.250  40.390  1.00  0.00             C\n",
      "ATOM   2698  C18 CHOL   21      12.510  36.490  38.440  1.00  0.00             C\n",
      "ATOM   2699  C19 CHOL   21      13.680  35.670  38.990  1.00  0.00             C\n",
      "ATOM   2700  C20 CHOL   21      14.070  36.390  40.280  1.00  0.00             C\n",
      "grep: write error: Broken pipe\n"
     ]
    }
   ],
   "source": [
    "!grep CHOL traj_POPC_wH.pdb | head -20"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "046dfef0",
   "metadata": {},
   "source": [
    "Yeah that is fine! POPC has all hydrogens reconstructed, but not cholesterol (we can only see the polar hydrogen on the OH group which is natively present in a united-atom structure)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b52fadb7",
   "metadata": {},
   "source": [
    "## Dealing with cholesterol\n",
    "\n",
    "Now that POPC have all hydrogens, we can restart from that pdb file to reconstruct the hydrogens of cholesterol. Thus we can use the flags `-c` with the pdb file we just generated `traj_POPC_wH.pdb` (Hs on POPC but cholesterol). The `-t` will be used with the same trajectory with H on POPC (but not POPE!) `traj_POPC_wH.xtc`. We output order parameters for POPE with `-o`. Last, `-opx` will produce a file with hydrogens on POPE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "713f658d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Constructing the system...\n",
      "System has 32914 atoms\n",
      "Writing new pdb with hydrogens.\n",
      "Writing trajectory with hydrogens in xtc file.\n",
      "Dealing with frame 0 at 0.0 ps.\n",
      "Dealing with frame 1 at 0.0 ps.\n",
      "Dealing with frame 2 at 0.0 ps.\n",
      "Dealing with frame 3 at 0.0 ps.\n",
      "Dealing with frame 4 at 0.0 ps.\n",
      "Dealing with frame 5 at 0.0 ps.\n",
      "Dealing with frame 6 at 0.0 ps.\n",
      "Dealing with frame 7 at 0.0 ps.\n",
      "Dealing with frame 8 at 0.0 ps.\n",
      "Dealing with frame 9 at 0.0 ps.\n",
      "Dealing with frame 10 at 0.0 ps.\n",
      "Dealing with frame 11 at 0.0 ps.\n",
      "Dealing with frame 12 at 0.0 ps.\n",
      "Dealing with frame 13 at 0.0 ps.\n",
      "Dealing with frame 14 at 0.0 ps.\n",
      "Dealing with frame 15 at 0.0 ps.\n",
      "Dealing with frame 16 at 0.0 ps.\n",
      "Dealing with frame 17 at 0.0 ps.\n",
      "Dealing with frame 18 at 0.0 ps.\n",
      "Dealing with frame 19 at 0.0 ps.\n",
      "Dealing with frame 20 at 0.0 ps.\n",
      "Dealing with frame 21 at 0.0 ps.\n",
      "Dealing with frame 22 at 0.0 ps.\n",
      "Dealing with frame 23 at 0.0 ps.\n",
      "Dealing with frame 24 at 0.0 ps.\n",
      "Dealing with frame 25 at 0.0 ps.\n",
      "Results written to OP_CHOL.out\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "!buildH -c traj_POPC_wH.pdb -t traj_POPC_wH.xtc -l Berger_CHOL -d Berger_CHOL.def -o OP_CHOL.out -opx traj_POPC_CHOL_wH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c4f03dc7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Berger_CHOL.def\n",
      "Berger_POPC.def\n",
      "data\n",
      "endCONF_OK.gro\n",
      "img\n",
      "Notebook_01_buildH_calc_OP.ipynb\n",
      "Notebook_02_buildH_calc_OP_outputwH.ipynb\n",
      "Notebook_03_mixture_POPC_POPE.ipynb\n",
      "Notebook_04_library.ipynb\n",
      "Notebook_05_mixture_POPC_cholesterol.ipynb\n",
      "OP_CHOL.out\n",
      "OP_POPC.out\n",
      "popcCHOL34molPER0-25ns_OK_dt1000.xtc\n",
      "traj_POPC_CHOL_wH.pdb\n",
      "traj_POPC_CHOL_wH.xtc\n",
      "traj_POPC_wH.pdb\n",
      "traj_POPC_wH.xtc\n"
     ]
    }
   ],
   "source": [
    "!ls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c606845",
   "metadata": {},
   "source": [
    "Great, we can now look at the order parameters of cholesterol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "837ea0fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# OP_name            resname atom1 atom2  OP_mean OP_stddev OP_stem\n",
      "#--------------------------------------------------------------------\n",
      "cholesterol1a        CHOL    C1    H11   -0.13027  0.25110  0.03785\n",
      "cholesterol1b        CHOL    C1    H12   -0.25855  0.18820  0.02837\n",
      "cholesterol1c        CHOL    C1    H13    0.77977  0.24458  0.03687\n",
      "cholesterol3a        CHOL    C3    H31   -0.27535  0.24273  0.03659\n",
      "cholesterol3b        CHOL    C3    H32   -0.28172  0.23349  0.03520\n",
      "cholesterol4a        CHOL    C4    H41   -0.27938  0.21471  0.03237\n",
      "cholesterol4b        CHOL    C4    H42   -0.30714  0.14779  0.02228\n",
      "cholesterol5a        CHOL    C5    H51   -0.30473  0.20913  0.03153\n"
     ]
    }
   ],
   "source": [
    "!head OP_CHOL.out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b7d076c",
   "metadata": {},
   "source": [
    "And check whether POPC **and** POPE have all hydrogens reconstructed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "73fc7da5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATOM   2681  C1  CHOL   21       8.750  35.790  36.370  1.00  0.00             C\n",
      "ATOM   2682  H11 CHOL   21       8.249  36.375  37.141  1.00  0.00             H\n",
      "ATOM   2683  H12 CHOL   21       9.238  34.929  36.827  1.00  0.00             H\n",
      "ATOM   2684  H13 CHOL   21       8.017  35.447  35.640  1.00  0.00             H\n",
      "ATOM   2685  C2  CHOL   21       9.800  36.660  35.670  1.00  0.00             C\n",
      "ATOM   2686  C3  CHOL   21       9.370  37.990  35.030  1.00  0.00             C\n",
      "ATOM   2687  H31 CHOL   21      10.267  38.580  34.839  1.00  0.00             H\n",
      "ATOM   2688  H32 CHOL   21       8.733  38.517  35.740  1.00  0.00             H\n",
      "ATOM   2689  C4  CHOL   21       8.600  37.850  33.710  1.00  0.00             C\n",
      "ATOM   2690  H41 CHOL   21       7.627  37.405  33.917  1.00  0.00             H\n",
      "ATOM   2691  H42 CHOL   21       8.462  38.841  33.277  1.00  0.00             H\n",
      "ATOM   2692  C5  CHOL   21       9.340  36.970  32.700  1.00  0.00             C\n",
      "ATOM   2693  H51 CHOL   21      10.245  37.528  32.460  1.00  0.00             H\n",
      "ATOM   2694  O6  CHOL   21       8.500  36.770  31.560  1.00  0.00             O\n",
      "ATOM   2695  H   CHOL   21       8.170  37.670  31.260  1.00  0.00             H\n",
      "ATOM   2696  C8  CHOL   21       9.890  35.630  33.200  1.00  0.00             C\n",
      "ATOM   2697  H81 CHOL   21      10.667  35.280  32.521  1.00  0.00             H\n",
      "ATOM   2698  H82 CHOL   21       9.084  34.898  33.237  1.00  0.00             H\n",
      "ATOM   2699  C9  CHOL   21      10.480  35.810  34.600  1.00  0.00             C\n",
      "ATOM   2700  C10 CHOL   21      11.690  35.180  34.860  1.00  0.00             C\n",
      "grep: write error: Broken pipe\n"
     ]
    }
   ],
   "source": [
    "!grep CHOL traj_POPC_CHOL_wH.pdb | head -20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d0ac2551",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATOM      1  C1   PLA    1      10.670  42.050  27.850  1.00  0.00             C\n",
      "ATOM      2  H11  PLA    1      10.730  41.120  27.280  1.00  0.00             H\n",
      "ATOM      3  H12  PLA    1       9.960  41.930  28.670  1.00  0.00             H\n",
      "ATOM      4  H13  PLA    1      11.650  42.290  28.250  1.00  0.00             H\n",
      "ATOM      5  C2   PLA    1       8.900  42.780  26.440  1.00  0.00             C\n",
      "ATOM      6  H21  PLA    1       8.970  41.830  25.920  1.00  0.00             H\n",
      "ATOM      7  H22  PLA    1       8.560  43.550  25.750  1.00  0.00             H\n",
      "ATOM      8  H23  PLA    1       8.190  42.690  27.260  1.00  0.00             H\n",
      "ATOM      9  C3   PLA    1      11.120  43.130  25.800  1.00  0.00             C\n",
      "ATOM     10  H31  PLA    1      11.120  42.140  25.350  1.00  0.00             H\n",
      "ATOM     11  H32  PLA    1      12.130  43.390  26.120  1.00  0.00             H\n",
      "ATOM     12  H33  PLA    1      10.770  43.860  25.070  1.00  0.00             H\n",
      "ATOM     13  N4   PLA    1      10.220  43.140  26.970  1.00  0.00             N\n",
      "ATOM     14  C5   PLA    1      10.120  44.420  27.670  1.00  0.00             C\n",
      "ATOM     15  H51  PLA    1       9.430  44.290  28.510  1.00  0.00             H\n",
      "ATOM     16  H52  PLA    1       9.700  45.150  26.970  1.00  0.00             H\n",
      "ATOM     17  C6   PLA    1      11.420  45.010  28.230  1.00  0.00             C\n",
      "ATOM     18  H61  PLA    1      12.060  44.190  28.560  1.00  0.00             H\n",
      "ATOM     19  H62  PLA    1      11.180  45.650  29.080  1.00  0.00             H\n",
      "ATOM     20  O7   PLA    1      12.130  45.790  27.260  1.00  0.00             O\n",
      "grep: write error: Broken pipe\n"
     ]
    }
   ],
   "source": [
    "!grep PLA traj_POPC_CHOL_wH.pdb | head -20"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16a5a4c4",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "In this notebook, we showed you how to use **buildH** on a mixture of POPC/cholesterol lipids. **buildH** has to be launched for each lipid separately. If a trajectory with all hydrogens on all lipids is wanted, one has to produce intermediate trajectories at each step. In the figure below is shown a snapshot with all the reconstructed hydrogens on cholesterol molecules (POPC and water are not shown for clarity, image rendered with [VMD](https://www.ks.uiuc.edu/Research/vmd/)).\n",
    "\n",
    "![H reconstruction on a POPC/POPE mixture with buildH](img/cholesterols.png)"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}